Changeset 68

Timestamp:

08/30/11 16:13:04 (13 years ago)

Author:

berrada

Message:

commit final avec filtre

Location:

branches/branche-mb

Files:

• data_in/file_rest/GYRE_00000400_restart_CONTROL_T.nc (modified) (previous)
• scripts/sinobad.i (modified) (2 diffs)
• scripts/sinobad_2D.i (modified) (4 diffs)
• src/cost.dat (modified) (1 diff)
• src/dta_lfexp.h (modified) (4 diffs)
• src/dta_zfimp.h (added)
• src/dta_zfimp_t.h (added)
• src/sinobad.d (modified) (5 diffs)
• src/sinobad.h (modified) (4 diffs)
• src/ta_c.h (modified) (1 diff)

branches/branche-mb/scripts/sinobad.i

@@ -18 +18 @@
 xwriteout 0 true ../data_out/exp_T.nc
 #setstate nu 1
-#init_nu
+init_nu
+xset_fdt_fco th 100. 
+xset_fdt_fco tz 100.
+xset_fdt_fco kh 5.e+02
+xset_fdt_fco kz 1.2e-05
+
 print_time OFF
 set_modeltime 0
@@ -36 +41 @@
 savestate ub  1   ijk   50    A       1       ../data_out/ub_obs_48.dat   
 savestate vb  1   ijk   50    A       1       ../data_out/vb_obs_48.dat   
-#GOTO TESTAD 
+GOTO TESTAD 
 
 loadobs sshb  1   ij   50    A       1       ../data_out/sshb_obs_48.dat D   

branches/branche-mb/scripts/sinobad_2D.i

@@ -2 +2 @@
 #               fichier de D'INSTRUCTIONS de référence
 #Julien Brajard
+# modified by mohamed Berrada (filtre August 2011)
 #==============================================================================
 #script faisant l'assimilation de référence
@@ -14 +15 @@
 load_eb ta_c   // pour tenir en compte ta_c sinon sera 0 cause filtre
 true_target_in_tab ta_c
+
+xset_fdt_fco th 1000. 
+xset_fdt_fco tz 10.
+xset_fdt_fco kh 5.e+05
+xset_fdt_fco kz 1.2e-03
+
 xinitnc ../data_out/exp_T_2D.nc
 
@@ -40 +47 @@
 setstate nu 0
 outoebx nu 1 0
-set_bcoef nu 1.e-05
+set_bcoef nu 1.e-8
 set_modeltime 0
 print_time OFF
@@ -69 +76 @@
 goto nomulti
 load_eb ta_c  // charger l'ebauche '
-set_bcoef nu 1.e-06
+set_bcoef nu 5.e-05
 runm
 
 load_eb ta_c  // charger l'ebauche '
-set_bcoef nu 1.e-07
+set_bcoef nu 1.e-04
 runm
+
+##SAUVEGARDE
+xwriteout 0 fin ../data_out/exp_T_2D.nc
+xwriteout 50 forwfin_48 ../data_out/exp_T_2D.nc
+
+xrst_save ../data_in/file_rest/GYRE_00000400_restart_CONTROL_T_2D.nc
+
+
 nomulti
 

branches/branche-mb/src/cost.dat

@@ -2 +2 @@
  0.0000000000000000e+00
  0.0000000000000000e+00
+ 2.1474823164582097e+02
+ 1.2274234090104164e+01
+ 2.3881198434552111e+01
+ 8.8771059162513382e+00
+ 5.5557975737505458e+00
+ 3.4803779742760317e+00
+ 1.7128978455386892e+00
+ 1.1973455752757605e+00
+ 9.8335187344394392e-01
+ 7.6117334075192200e-01
+ 5.9846985043808543e-01
+ 6.7366624706570466e-01
+ 5.5874401856066103e-01
+ 7.3595639204467600e-01
+ 5.0953722678532765e-01
+ 1.1965283691171660e+00
+ 6.5986339544255534e-01
+ 6.6013045954860117e-01
+ 6.5392286280927792e-01
+ 5.6394373996837344e-01
+ 5.7610237016601717e-01
+ 5.8032116526810040e-01
+ 5.8309360765246787e-01
+ 5.8519910883568982e-01
+ 5.8664512715255190e-01
+ 5.8761359807265934e-01
+ 5.8831541023256795e-01
+ 5.8880781610867361e-01
+ 5.8915693207370023e-01
+ 5.0893021741293576e-01
+ 5.0885469878063883e-01
+ 5.8920763747258753e-01
+ 5.0875296736771736e-01
+ 5.0874187732734499e-01
+ 5.8921509431560293e-01
+ 5.8922030000369108e-01
+ 5.8922395197080846e-01
+ 5.8922651344473387e-01

branches/branche-mb/src/dta_lfexp.h

@@ -37 +37 @@
         tim1= 0.;
 #ifdef K_NEU
-        tim1= - bim1 * x2;
+        tim1= bim1 * x2;
 #endif
       }
       else{
-        tim1= -bim1 * x3;
+        tim1= bim1 * x3;
       }
       if (Yi==NX-2){
         tip1= 0.;
 #ifdef K_NEU
-        tip1= - bi * x2;
+        tip1= bi * x2;
 #endif
       }
       else{
-        tip1= -bi * x4;
+        tip1= bi * x4;
       }
       //cdt au bord j
@@ -56 +56 @@
         tjm1= 0.;
 #ifdef K_NEU
-        tjm1= - bjm1 * x2;
+        tjm1= bjm1 * x2;
 #endif
       }
       else{
-        tjm1= -bjm1 * x5;
+        tjm1= bjm1 * x5;
       }
       if (Yj==NY-2){
         tjp1= 0.;
 #ifdef K_NEU
-        tjp1= - bj * x2;
+        tjp1= bj * x2;
 #endif
       }
       else{
-        tjp1= -bj * x6;
+        tjp1= bj * x6;
       }
 
-      tij= (1./zbtr+bi+bim1+bj+bjm1) * x2;
+      tij= (1./zbtr-bi-bim1-bj-bjm1) * x2;
       YS1 = zbtr*(tip1 + tim1 + tjp1 + tjm1 + tij);
     }
@@ -105 +105 @@
         dtim1_x2= 0.;
 #ifdef K_NEU
-        dtim1_x2= - bim1;
+        dtim1_x2= bim1;
 #endif
       }
       else{
-        dtim1_x3= -bim1;
+        dtim1_x3= bim1;
       }
       if (Yi==NX-2){
         dtip1_x2= 0.;
 #ifdef K_NEU
-        dtip1_x2= - bi;
+        dtip1_x2= bi;
 #endif
       }
       else{
-        dtip1_x4= -bi;
+        dtip1_x4= bi;
       }
       //cdt au bord j
@@ -124 +124 @@
         dtjm1_x2= 0.;
 #ifdef K_NEU
-        dtjm1_x2= - bjm1;
+        dtjm1_x2= bjm1;
 #endif
       }
       else{
-        dtjm1_x5= -bjm1;
+        dtjm1_x5= bjm1;
       }
       if (Yj==NY-2){
         dtjp1_x2= 0.;
 #ifdef K_NEU
-        dtjp1_x2= - bj;
+        dtjp1_x2= bj;
 #endif
       }
       else{
-        dtjp1_x6= -bj;
+        dtjp1_x6= bj;
       }
-      dtij_x2= (1./zbtr+bi+bim1+bj+bjm1);
+      dtij_x2= (1./zbtr-bi-bim1-bj-bjm1);
 
       YJ1I1 = 0.;

branches/branche-mb/src/sinobad.d

@@ -1 @@
-defval K_NEU
+//defval K_NEU
 #define FZIMP
 #define FILTER
-//#define FILZNOL
-#define FILLNOZ
+#define FILZNOL
+//#define FILLNOZ
 //#define SOLSORYAO
 //#define OPTIMORDER
@@ -22 +22 @@
 defval TU       2 | start time step
 #ifdef FILTER
-defval OFTL      10
-defval OFTZ      0
-defval OFT       10
+defval OFTL      0
+defval OFTZ      4
+defval OFT       4
 #endif
 #ifndef FILTER
@@ -329 +329 @@
 //#---------->dta_zfil
 #ifndef FILLNOZ
+
 #ifndef FILZNOL
 #ifndef FZIMP
@@ -367 +368 @@
 ctin     ta_c           1       from    dta_zfexp               1  i    j    k   t
 #endif
-#endif
+
+#endif
+
 #ifdef FILLNOZ
 ctin     ta_c           1       from    dta_lfexp               1  i    j    k   t
@@ -1209 +1212 @@
 #ifdef FILTER
 insert_fct  arg load_eb
-insert_fct  init_nu
+insert_fct  arg init_nu
+insert_fct  arg xset_fdt_fco
 #endif
 insert_fct arg xwriteout

branches/branche-mb/src/sinobad.h

@@ -87 +87 @@
 double pdth_fil=100.;//7200.;
 double pdtz_fil=100.;//100.;//7200.;
-double fsaht_fil= 5.e+2;//1000.;
+double fsaht_fil= 5.e+02;//1000.;
 double avt_fil=1.2e-05;//1.2e-5;
 //double lscale= 212000.; // length-scale
 double norm_fil; //coef de normalisation
 void normfil();// calcul du coef de normalisation norm_fil
+void xset_fdt_fco();// set fdt fcoef
 #endif
 
@@ -259 +260 @@
 //_____________________________________________________________________________
 #ifdef K_FILTER
+//_____________________________________________________________________________
 void normfil(){
   double norm_zfil=1.;
   double norm_lfil=1.;
-#if defined (K_FILZNOL) /* && !defined (K_FILLNOZ)*/ 
-  norm_zfil= 2.*sqrt(PI*avt_fil*pdtz_fil*2.*OFTZ);
-#ifdef K_FZIMP
+  norm_lfil= 2.*sqrt(PI*fsaht_fil*pdth_fil*2.*(OFTL-1.));// schema explicite
+  norm_zfil= 2.*sqrt(PI*avt_fil*pdtz_fil*2.*(OFTZ-1.));// schema explicite 
+#ifdef K_FZIMP  // schema implicite 
   double fact=1.; 
   int M=2*(OFTZ-1);
@@ -271 +273 @@
   norm_zfil=2.*fact*sqrt(avt_fil*pdtz_fil);
 #endif
-#elif  defined (K_FILLNOZ) /*&& !defined (K_FILZNOL)*/  
-  norm_lfil= 2.*sqrt(PI*fsaht_fil*pdth_fil*2.*OFTL);
-#endif
+
+#if defined (K_FILZNOL)
+  norm_lfil=1.;
+#endif
+#if  defined (K_FILLNOZ)
+  norm_zfil=1.;
+#endif
+
   norm_fil= norm_zfil*norm_lfil;
   printf("norm_fil = %23.16e\n",norm_fil);
 }
 //_____________________________________________________________________________
-void init_nu(){
-
-  for(int tt=0;tt<24;tt++)
+void xset_fdt_fco(int argc, char *argv[]){// set fdt fcoef
+  if (argc!=3){
+    printf("\n inappropriate arguments choice in command xset_fdt_fco\n");
+    exit(99);
+  }
+  if(!strcmp(argv[1],"th")){
+    pdth_fil=atof(argv[2]);
+  }
+  else if(!strcmp(argv[1],"tz")){
+    pdtz_fil=atof(argv[2]);
+  }
+  else if(!strcmp(argv[1],"kh")){
+    fsaht_fil=atof(argv[2]);
+  }
+  else if(!strcmp(argv[1],"kz")){
+    avt_fil=atof(argv[2]);
+  }
+  else{
+    printf("\n inappropriate first argument choice in command xset_fdt_fco\n");
+    exit(99);
+  }
+}
+//_____________________________________________________________________________
+void init_nu(int argc, char *argv[]){
+  if (argc==1){
+    for(int tt=0;tt<24;tt++)
+      for(int k=0;k<NZ;k++)
+        for(int j=0;j<NY;j++)
+          for(int i=0;i<NX;i++){
+            YS_nu(0,i,j,k)=YS_nu(0,i,j,k)+0.1*exp(-((i+1.-(3.*NX)/4.)*(i+1.-(3.*NX)/4.)+ 
+                                                    (j+1.-(3.*NY)/4.)*(j+1.-(3.*NY)/4.))/4.)
+              *(1.-exp(-double(k+1.-NZ)/NZ))*(1.-exp(-double(tt)/24.));
+          }
+  }
+  else if  (argc==5){
     for(int k=0;k<NZ;k++)
       for(int j=0;j<NY;j++)
         for(int i=0;i<NX;i++){
-          YS_nu(0,i,j,k)=YS_nu(0,i,j,k)+0.1*exp(-((i+1.-(3.*NX)/4.)*(i+1.-(3.*NX)/4.)+ 
-                                              (j+1.-(3.*NY)/4.)*(j+1.-(3.*NY)/4.))/4.)
-            *(1.-exp(-double(k+1.-NZ)/NZ))*(1.-exp(-double(tt)/24.));
-            }
+          YS_nu(0,i,j,k)=0.;
+          }
+
+    float v=atof(argv[1]);
+    int ji=atoi(argv[2])-1;
+    int jj=atoi(argv[3])-1;
+    int jk=atoi(argv[4])-1;
+    YS_nu(0,ji,jj,jk)=v;
+  }
+  else{
+    printf("\n inappropriate arguments choice in command init_nu\n");
+    exit(99);
+  }
 }
 //_____________________________________________________________________________
@@ -300 +348 @@
       }
   double a0,a1;
+  double a2,b2,c2;
   for(int j=1;j<NY-1;j++)
     for(int i=1;i<NX-1;i++){
       a1=avt_fil*pdtz_fil/fse3w(i,j,1);
-      tab_luzimp( i, j, 0 , 0 ) = 0.;// 0-->gam
-      tab_luzimp( i, j, 0 , 1 ) = fse3t(i,j,0)+a1;// 1-->bet
+      b2= 1.-a1/fse3t(i,j,0);
 #ifndef K_NEUM
       a0=avt_fil*pdtz_fil/fse3w(i,j,0);
-      tab_luzimp( i, j, 0 , 1 ) += a0;// 1-->bet
-#endif
+      b2 -= (a0/fse3t(i,j,0));// 1-->bet
+#endif
+      tab_luzimp( i, j, 0 , 0 ) = 0.;// 0-->gam
+      tab_luzimp( i, j, 0 , 1 ) = b2;// 1-->bet
       for(int k=1;k<NZ-2;k++){
         a0=avt_fil*pdtz_fil/fse3w(i,j,k);
         a1=avt_fil*pdtz_fil/fse3w(i,j,k+1);
-        tab_luzimp( i, j, k , 0 ) = -a0/tab_luzimp( i, j, k-1 , 1);
-        tab_luzimp( i, j, k , 1 ) = fse3t(i,j,k)+a1+a0 + a0 * tab_luzimp( i, j, k , 0 );
+        a2 = a0/fse3t(i,j,k);
+        b2 =  1.-a1/fse3t(i,j,k)-a0/fse3t(i,j,k);
+        c2 = a0/fse3t(i,j,k-1);
+        tab_luzimp( i, j, k , 0 ) = c2/tab_luzimp( i, j, k-1 , 1);
+        tab_luzimp( i, j, k , 1 ) = b2 - a2 * tab_luzimp( i, j, k , 0 );
       }
       a0=avt_fil*pdtz_fil/fse3w(i,j,NZ-2);
-      tab_luzimp( i, j, NZ-2 , 0 ) = -a0/tab_luzimp( i, j, NZ-3 , 1);
-      tab_luzimp( i, j, NZ-2 , 1 ) = fse3t(i,j,NZ-2)+a0 + a0 * tab_luzimp( i, j, NZ-2 , 0 );
+      a2 = a0/fse3t(i,j,NZ-2);
+      b2 =  1.-a0/fse3t(i,j,NZ-2);
 #ifndef K_NEUM
       a1=avt_fil*pdtz_fil/fse3w(i,j,NZ-1);
-      tab_luzimp( i, j, NZ-2 , 1 ) += a1;
-#endif
-    }
-  /*  for(int k=0;k<NZ;k++){ // init a 0
-    printf("gam(%d)=%f \t",k+1,tab_luzimp( 5, 5, k , 0 ));
-    printf("bet(%d)=%f \n",k+1,tab_luzimp( 5, 5, k , 1 ));
-    }*/
+      b2 -= a1/fse3t(i,j,NZ-2);
+#endif
+      c2 = a0/fse3t(i,j,NZ-3);
+     
+      tab_luzimp( i, j, NZ-2 , 0 ) = c2/tab_luzimp( i, j, NZ-3 , 1);
+      tab_luzimp( i, j, NZ-2 , 1 ) = b2 - a2 * tab_luzimp( i, j, NZ-2 , 0 );
+    }
+  for(int k=0;k<NZ;k++){ // init a 0
+    printf("gam(%d)=%f \t",k+1,tab_luzimp( 23, 15, k , 0 ));
+    printf("bet(%d)=%f \n",k+1,tab_luzimp( 23, 15, k , 1 ));
+    }
 }
 #endif

branches/branche-mb/src/ta_c.h

@@ -20 +20 @@
   ////
   //  printf("norm_fac= %24.16e\n",norm_fac);
-  if(fabs(norm_fil*x1)>5.) 
+  if(fabs(norm_fil*x1)>5) 
       printf("increm(%d,%d,%d) = %24.16e (%24.16e) \n",Yi+1,Yj+1,Yk+1,norm_fil*x1,YS1);
 }