Changeset 62


Ignore:
Timestamp:
07/29/11 16:44:58 (13 years ago)
Author:
berrada
Message:

filtre vertical-schema explicite- et normalisation

Location:
branches/branche-mb
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • branches/branche-mb/data_out/plot_results.m

    r38 r62  
    1 f=netcdf('exp_T.nc'); 
     1f=netcdf('exp_T_2D.nc'); 
    22 
    33addpath('../scripts/matlab_toolbox'); 

  • branches/branche-mb/scripts/sinobad.i

    r59 r62  
    1717 
    1818xwriteout 0 true ../data_out/exp_T.nc 
    19  
     19#setstate nu 1 
     20#init_nu 
    2021print_time OFF 
    2122set_modeltime 0 
    … …  
    3536savestate ub  1   ijk   50    A       1       ../data_out/ub_obs_48.dat    
    3637savestate vb  1   ijk   50    A       1       ../data_out/vb_obs_48.dat    
     38#GOTO TESTAD  
    3739 
    3840loadobs sshb  1   ij   50    A       1       ../data_out/sshb_obs_48.dat D    
    … …  
    4244loadobs vb  1   ijk   50    A       1       ../data_out/vb_obs_48.dat D    
    4345 
    44  
     46GOTO FINTEST  
     47#TEST DE L'ADJOINT' 
     48TESTAD 
     49testad  1.e-01 
     50stop 
     51#----------------- 
     52#TEST DE LA PROGRAMMATION DES DERIVEES 
     53TESTDF 
     54testdf 3  3  1  4    f 1.e-10  0.001  10 
     55stop 
     56#----------------- 
     57# TEST DE LINEAIRE TANGENT 
     58TESTLT 
     59 testlt 1.e+00 0.5e-01 2. 15 
     60stop 
     61FINTEST 
    4562#goto FINRUN 
    4663 
    … …  
    5067setstate nu 0 
    5168outoebx nu 1 0 
     69set_bcoef nu 1.e-16 
    5270set_modeltime 0 
    5371print_time OFF 
    … …  
    6381FINTEST 
    6482##DEBUT DE L'ASSIMILATION' 
    65  
     83print_cost ON 
    6684set_modeltime 0 
    6785!echo "RUN OPTIMIZATION WITH M1QN3 ......." 

  • branches/branche-mb/scripts/sinobad_2D.i

    r34 r62  
    1212##INITIALISATION 
    1313xistate_init 1 ../data_in/file_rest/GYRE_00000424_restart.nc 
     14load_eb ta_c   // pour tenir en compte ta_c sinon sera 0 cause filtre 
    1415true_target_in_tab ta_c 
    1516xinitnc ../data_out/exp_T_2D.nc 
    … …  
    3637##INITIALISATION DE L'ASSIMILATION' 
    3738xistate_init 1 ../data_in/file_rest/GYRE_00000400_restart.nc 
     39load_eb ta_c  // charger l'ebauche ' 
     40setstate nu 0 
     41outoebx nu 1 0 
     42set_bcoef nu 1.e-05 
    3843set_modeltime 0 
    3944print_time OFF 
    … …  
    4954FINTEST 
    5055##DEBUT DE L'ASSIMILATION' 
    51  
     56print_cost ON 
    5257set_modeltime 0 
    5358!echo "RUN OPTIMIZATION WITH M1QN3 ......." 
    … …  
    6166setm_ddf1    1 
    6267runm 
     68 
     69goto nomulti 
     70load_eb ta_c  // charger l'ebauche ' 
     71set_bcoef nu 1.e-06 
     72runm 
     73 
     74load_eb ta_c  // charger l'ebauche ' 
     75set_bcoef nu 1.e-07 
     76runm 
     77nomulti 
     78 
    6379##SAUVEGARDE 
    6480xwriteout 0 fin ../data_out/exp_T_2D.nc 

  • branches/branche-mb/src/Makefile

    r55 r62  
    1 #GENOPT=-DOPTIMORDER -DSOLSORYAO -E 
     1# GENOPT=-DOPTIMORDER -DSOLSORYAO -E 
    22#GENOPT= -E 
    33GENOPT=-DOPTIMORDER -E 

  • branches/branche-mb/src/cost.dat

    r59 r62  
    22 0.0000000000000000e+00 
    33 0.0000000000000000e+00 
     4 2.1474823164582097e+02 
     5 1.2975758308681927e+01 
     6 2.4264061059597573e+01 
     7 9.5294911320411000e+00 
     8 6.1849686337309064e+00 
     9 4.0738038971828123e+00 
     10 2.3898782803950009e+00 
     11 2.0888660877807852e+00 
     12 1.6129218002581345e+00 
     13 1.4602912585844716e+00 
     14 1.5154815138350144e+00 
     15 1.4661975187573568e+00 
     16 1.4668321473655250e+00 
     17 1.4708862822664150e+00 
     18 1.4854466172307939e+00 
     19 1.4205507107251405e+00 
     20 1.7030517837610164e+00 
     21 1.4345968631501795e+00 
     22 1.4351405506167443e+00 
     23 1.4397453657339783e+00 
     24 1.4537877548984335e+00 
     25 1.4674691588039481e+00 
     26 1.4760376865940141e+00 
     27 1.4825757968508098e+00 
     28 1.4094545811622352e+00 
     29 1.4080263806488316e+00 
     30 1.4835348660561365e+00 
     31 1.4064118013093134e+00 
     32 1.4062029823934583e+00 
     33 1.4836757396029698e+00 
     34 1.4837748681619887e+00 
     35 1.4838440857485093e+00 
     36 1.4838925284218420e+00 
     37 1.4839264406875585e+00 
     38 1.4061456937219163e+00 

  • branches/branche-mb/src/dta_lfil.h

    r59 r62  
    1818 
    1919  if(Yt==1){ 
    20     YS1=x1; 
     20    double zdep; 
     21    if (Yk==0)  
     22      zdep=fsdept( Yi , Yj , 1 ) - fsdept( Yi , Yj , 0 ); 
     23    else     
     24      zdep=fsdept( Yi , Yj , Yk ) - fsdept( Yi , Yj , Yk-1 ); 
     25    double wm12=1./(106000*sqrt(zdep)); //1/sqrt(ds1*ds2*ds3) 
     26    YS1=wm12*x1; 
    2127  } 
    2228  else{ 
    … …  
    3036  } 
    3137  // 
     38    if(Yt==1 && Yk==0 && Yi<26 && Yi>22 && Yj<17 && Yj>13)  
     39      printf("nul1(%d,%d) = %24.16e \n",Yi+1,Yj+1,YS1); 
     40    if(Yt==40 && Yk==0 && Yi<26 && Yi>22 && Yj<17 && Yj>13)  
     41      printf("nul40(%d,%d) = %24.16e \n",Yi+1,Yj+1,YS1); 
    3242} 
    3343 
    … …  
    4151 
    4252  if(Yt==1){ 
    43     YJ1I1=1.; 
     53    double zdep; 
     54    if (Yk==0)  
     55      zdep=fsdept( Yi , Yj , 1 ) - fsdept( Yi , Yj , 0 ); 
     56    else     
     57      zdep=fsdept( Yi , Yj , Yk ) - fsdept( Yi , Yj , Yk-1 ); 
     58    double wm12=1./(106000*sqrt(zdep)); //1/sqrt(ds1*ds2*ds3) 
     59    YJ1I1=wm12; 
    4460  } 
    4561  else{ 
    … …  
    5470  } 
    5571  // 
     72  if(Ycurward==BACKWARD){ 
     73    if(Yt==10 && Yk==0 && Yi<26 && Yi>22 && Yj<17 && Yj>13)  
     74      printf("YG_1L(%d,%d) = %24.16e \n",Yi+1,Yj+1,YG1); 
     75  } 
    5676} 

  • branches/branche-mb/src/dta_zfil.h

    r59 r62  
    1515  //printf("dta_zfil : %d\n",Yt); 
    1616  if(Yt==1){ 
    17     YS1=x1; 
     17    double wm12=1./sqrt(fse3t(Yi,Yj,Yk)); //1/sqrt(dz) 
     18    YS1=wm12*x1; 
    1819  } 
    1920  else{ 
    … …  
    2627    } 
    2728  } 
     29  /*    if(Yt==1 && Yk==0 && Yi<26 && Yi>22 && Yj<17 && Yj>13)  
     30printf("nuz1(%d,%d) = %24.16e \n",Yi+1,Yj+1,YS1); 
     31    if(Yt==40 && Yk==0 && Yi<26 && Yi>22 && Yj<17 && Yj>13)  
     32    printf("nuz40(%d,%d) = %24.16e \n",Yi+1,Yj+1,YS1);*/ 
    2833  // 
    2934} 
    … …  
    3641  YJ1I1=0.;     YJ1I2=0.;    YJ1I3=0.;    YJ1I4=0.; 
    3742  if(Yt==1){ 
    38     YJ1I1=1.; 
     43    double wm12=1./sqrt(fse3t(Yi,Yj,Yk)); //1/sqrt(dz) 
     44    YJ1I1=wm12; 
    3945  } 
    4046  else{ 
    … …  
    4955  // 
    5056 
     57  /* if(Ycurward==BACKWARD){ 
     58    if(Yt==1 && Yk==0 && Yi<26 && Yi>22 && Yj<17 && Yj>13)  
     59      printf("YG_1(%d,%d) = %24.16e \n",Yi+1,Yj+1,YG1); 
     60      }*/ 
    5161} 

  • branches/branche-mb/src/sinobad.d

    r60 r62  
    1 //#define FILTER 
     1#define FILTER 
     2#define FILZNOL 
     3//#define FILNOZL 
    24//#define SOLSORYAO 
    35//#define OPTIMORDER 
    … …  
    1820defval TU       2 | start time step 
    1921#ifdef FILTER 
    20 defval OFTL      40 
    21 defval OFTZ      40 
    22 defval OFT      80 
     22defval OFTL      0 
     23defval OFTZ      25 
     24defval OFT       25 
    2325#endif 
    2426#ifndef FILTER 
    … …  
    2830#ifdef FILTER 
    2931defval K_FILTER 
     32#endif 
     33#ifdef FILZNOL 
     34defval K_FILZNOL 
     35#endif 
     36#ifdef FILNOZL 
     37defval K_FILNOZL 
    3038#endif 
    3139#ifdef SOLSORYAO 
    … …  
    5361#ifdef FILTER 
    5462traj  Tcstf      M   0  1 
     63#ifndef FILNOZL 
     64#ifndef FILZNOL 
    5565traj  Tlfil      M   1  0    1  OFTL 
     66#endif 
    5667traj  Tzfil      M   1  OFTL 1  OFTZ 
     68#endif 
    5769#endif 
    5870traj  Tcst       M   0   OFT  1 1 
    … …  
    6274#ifdef FILTER 
    6375space  S3df      M    NX NY NZ    Tcstf 
     76#ifndef FILNOZL 
     77#ifndef FILZNOL 
    6478space  S3dtlf    M    NX NY NZ    Tlfil 
     79#endif 
    6580space  S3dtzf    M    NX NY NZ    Tzfil 
    6681#endif 
     82#endif 
    6783 
    6884space  S0d       M    1           Tcst 
    69 //#space  S3d    M    NX NY NZ    Tcst 
    7085space  S3d       M    NX NY NZ    Tcst 
    7186space  S3dt      M    NX NY NZ    Tsbd 
    … …  
    96111#ifdef FILTER 
    97112modul    nu                     space   S3df                    noward            output 1      target 
     113#ifndef FILNOZL 
     114#ifndef FILZNOL 
    98115modul    dta_lfil               space   S3dtlf                  input 6           output 1      tempo 
    99116modul    ztuv_lfil              space   S3dtlf                  input 3           output 2      tempo 
     117#endif 
    100118modul    dta_zfil               space   S3dtzf                  input 4           output 1      tempo 
    101119modul    zwy_zfil               space   S3dtzf                  input 2           output 1      tempo 
     120#endif 
    102121modul    ta_c                   space   S3d                     input 1           output 1       
    103122#endif 
    … …  
    285304 
    286305#ifdef FILTER 
     306#ifndef FILNOZL 
     307#ifndef FILZNOL 
    287308//#---------->dta_lfil 
    288309ctin     dta_lfil       1       from    nu                      1  i    j    k   
    … …  
    297318ctin     ztuv_lfil      2       from    dta_lfil                1  i    j    k     t-1  
    298319ctin     ztuv_lfil      3       from    dta_lfil                1  i    j+1  k     t-1  
     320#endif 
    299321 
    300322//#---------->dta_zfil 
     323#ifdef FILZNOL 
     324ctin     dta_zfil       1       from    nu                      1  i    j    k   
     325#endif 
     326#ifndef FILZNOL 
    301327ctin     dta_zfil       1       from    dta_lfil                1  i    j    k     t 
     328#endif 
    302329ctin     dta_zfil       2       from    dta_zfil                1  i    j    k     t-1  
    303330ctin     dta_zfil       3       from    zwy_zfil                1  i    j    k+1   t  
    … …  
    309336 
    310337//#---------->ta_c 
    311 ctin     ta_c           1       from    dta_zfil                1  i    j    k-1   t 
     338ctin     ta_c           1       from    dta_zfil                1  i    j    k   t 
     339#endif 
     340#ifdef FILNOZL 
     341ctin     ta_c           1       from    nu                      1  i    j    k    
     342#endif 
    312343#endif 
    313344//#ctin  ta_c                   1..NPCA from    pca_ta                  1..NPCA  1 
    314345 
    315 //exec disp_ct_in 
     346//#exec disp_ct_in 
    316347exec disp_modul 
    317348 
    … …  
    878909#ifdef FILTER 
    879910// calcul du filtre 
    880  
     911#ifndef FILNOZL 
     912#ifndef FILZNOL 
    881913order modinspace S3dtlf 
    882914      order YA3   
    … …  
    889921      forder 
    890922forder 
    891  
     923#endif 
    892924order modinspace S3dtzf 
    893925    order YA3 YA2 YA1 
    … …  
    898930    forder 
    899931forder 
     932#endif 
    900933order modinspace S3d 
    901934       order YA3 YA2 YA1 
    … …  
    11041137#ifdef FILTER 
    11051138 
     1139#ifndef FILNOZL 
     1140#ifndef FILZNOL 
    11061141order spaceintraj Tlfil 
    11071142        S3dtlf 
    11081143forder 
     1144#endif 
    11091145order spaceintraj Tzfil 
    11101146        S3dtzf 
    11111147forder 
     1148#endif 
    11121149order spaceintraj Tcst 
    11131150        S3d 
    … …  
    11311168insert_fct  arg xchangesavemode 
    11321169 
     1170#ifdef FILTER 
    11331171insert_fct  arg load_eb 
    1134  
     1172insert_fct  init_nu 
     1173#endif 
    11351174insert_fct arg xwriteout 
    11361175insert_fct arg xwritegrad 

  • branches/branche-mb/src/sinobad.h

    r59 r62  
    1313#include"../include/meshmask.h" 
    1414 
     15#ifdef K_FILTER 
    1516#define ta_eb( i, j, k )  (ta_eb[(i)*(NY*NZ)+(j)*(NZ)+(k)]) // ebauche 
    1617double ta_eb[NZ*NY*NX];   // ebauche 
    1718void load_eb (int argc, char *argv[]);  // load "ebauche" from current variable ta_c   
     19#endif 
    1820// For PPCA  
    1921/*double shfs_ta[NZ*NY*NX][NPCA]; 
    … …  
    7880 
    7981#ifdef K_FILTER 
    80 double pdth_fil=7200.; 
    81 double pdtz_fil=7200.; 
     82double pdth_fil=100.;//7200.; 
     83double pdtz_fil=100.;//7200.; 
    8284double fsaht_fil= 1.2e+2;//1000.; 
    83 double avt_fil=1.2e-1; 
     85double avt_fil=1.2e-5; 
     86double lscale= 212000.; // length-scale 
    8487#endif 
    8588 
    … …  
    106109  printf("//====================================================================//\n"); 
    107110  printf("\n"); 
     111#ifdef K_FILTER 
     112  printf("option FILTER enable\n"); 
     113#else 
     114  printf("option FILTER disable\n"); 
     115#endif 
    108116#ifdef K_SOLSORYAO 
    109117  printf("option SOLSORYAO enable\n"); 
    … …  
    236244  normsup_alla_c(); 
    237245} 
    238  
     246#ifdef K_FILTER 
     247void init_nu(){ 
     248 
     249  for(int tt=0;tt<24;tt++) 
     250    for(int k=0;k<NZ;k++) 
     251      for(int j=0;j<NY;j++) 
     252        for(int i=0;i<NX;i++){ 
     253          YS_nu(0,i,j,k)=YS_nu(0,i,j,k)+1.*exp(-((i+1.-(3.*NX)/4.)*(i+1.-(3.*NX)/4.)+  
     254                                              (j+1.-(3.*NY)/4.)*(j+1.-(3.*NY)/4.))/4.) 
     255            *(1.-exp(-double(k+1.-NZ)/NZ))*(1.-exp(-double(tt)/24.)); 
     256            } 
     257} 
     258#endif 
    239259//_____________________________________________________________________________ 
    240260#define hu( i, j )  (hu[(j)*(NX)+(i)]) 
    … …  
    540560} 
    541561//______________________________________________________________________________ 
     562#ifdef K_FILTER 
    542563void load_eb (int argc, char *argv[]){  // load "ebauche" 
    543564  // on doit presiser la variable exple : load_eb ta_c 
    … …  
    560581  } 
    561582} 
     583#endif 
    562584 
    563585 

  • branches/branche-mb/src/ta_c.h

    r59 r62  
    88forward (YREAL x1) 
    99{  
    10   YS1 = ta_eb(Yi,Yj,Yk) +x1; 
     10 
     11  //if(Yk==0 && Yi<14 && Yi>10 && Yj<14 && Yj>10) printf("nu(%d,%d) = %24.16e \n",Yi+1,Yj+1,x1); 
     12 
     13  double norm_fac; 
     14#ifdef K_FILZNOL 
     15  norm_fac= 2.*sqrt(PI*avt_fil*pdtz_fil*2.*OFTZ);///1.e-01;//1./sqrt(2.*PI)/lscale; 
     16#elif defined (K_FILNOZL)  
     17  norm_fac= 1.e+00; 
     18#else 
     19  norm_fac= 1.e+00; 
     20#endif 
     21  YS1 = ta_eb(Yi,Yj,Yk) +norm_fac*x1; 
     22  if(Yi==0 || Yi==NX-1 || Yj==0 || Yj==NY-1 || Yk==NZ-1){ 
     23    YS1=0.; 
     24  } 
    1125 
    1226 
    … …  
    1630  //YS1=sum; 
    1731  //// 
     32  /* if(fabs(x1)>0.1)  
     33     printf("increm(%d,%d,%d) = %24.16e (%24.16e) \n",Yi+1,Yj+1,Yk+1,x1,YS1);*/ 
    1834} 
    1935//=========================================================================== 
    … …  
    2238backward (YREAL x1) 
    2339{ 
    24   YJ1I1=1.;  
    25      
     40  double norm_fac; 
     41#ifdef K_FILZNOL 
     42  norm_fac= 2.*sqrt(PI*avt_fil*pdtz_fil*2.*OFTZ);///1.e-01;//1./sqrt(2.*PI)/lscale; 
     43#elif defined (K_FILNOZL)  
     44  norm_fac= 1.e+00; 
     45#else 
     46  norm_fac= 1.e+00; 
     47#endif 
     48  YJ1I1=norm_fac;  
     49  if(Yi==0 || Yi==NX-1 || Yj==0 || Yj==NY-1 || Yk==NZ-1){ 
     50    YJ1I1=0.; 
     51  } 
     52    
    2653  //for(int n=0;n<NPCA;n++)      
    2754  //  Yjac[0][n]    = stdev_ta[n]*shfs_ta(Yi,Yj,Yk,n); 
Note: See TracChangeset for help on using the changeset viewer.