Changeset 147 for altifloat


Ignore:
Timestamp:
04/30/15 11:39:49 (9 years ago)
Author:
jbrlod
Message:

add divergence term in cost function

Location:
altifloat
Files:
28 added
4 edited

Legend:

Unmodified
Added
Removed

  • altifloat/matlab_toolbox/script_perso.m

    r145 r147  
    11close all 
    22clear all 
    3 indir='../exp_forw_dan/'; 
     3indir='../exp_twin_aviso/'; 
    44%indir='../../code_leila/altifloat/obs_float/'; 
    55%indir='../obs_float/'; 
    … …  
    88robs=load([ indir 'obs_jum.dat']); 
    99rsim=load([ indir 'rfloat_total_jum1.dat']); 
    10 uvb=load([ indir 'uv_bck_dan_s.dat']); 
    11 uvt=load([ indir 'uv_bck_dan.dat']); 
     10uvb=load([ indir 'uv_bck.dat']); 
     11uvt=load([ indir 'uv_truth.dat']); 
    1212uvr=load([indir 'uv_total_jum1.dat']); 
    13 meshg=load([indir 'meshgrid_dan.dat']); 
    14 drifterdir='../altifloat/'; 
    15 fname={'a300234011043340.mat','a300234011120770.mat','a300234011249440.mat'}; 
     13meshg=load([indir 'meshgrid_aviso.dat']); 
     14 
     15%drifterdir='../altifloat/'; 
     16%fname={'a300234011043340.mat','a300234011120770.mat','a300234011249440.mat'}; 
    1617%X1=ut(:,end-3); 
    1718%Y1=ut(:,end-2); 
    18 iech=3; 
    19 jech=3; 
    20 LIM=[34.8 35.3 33.2 33.6]; 
     19iech=1; 
     20jech=1; 
     21LIM=[34.5 36 32.5 35]; 
    2122%Format : t/s, ilon, ilat, val 
    2223 
    … …  
    6667 
    6768for j=1:nfloat      
    68   data= load([drifterdir fname{j}]); 
    69    
    70  
    71    
     69   
     70   
     71 
    7272  rfloat=rsim(rsim(:,2)==j-1&rsim(:,end)>0,[end-1 end]); 
    7373  rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1); 
    7474  rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1); 
    75   Hsim=plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'.-r'); 
     75  Hsim=plot(rfloat_lon,rfloat_lat,'.-r'); 
    7676    rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]); 
    7777  rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1); 
    7878  rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1); 
    79   Hobs=plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'.-m'); 
     79  Hobs=plot(rfloat_lon,rfloat_lat,'.-m'); 
    8080 
    8181  %   plot(rfloat(:,2)+1,rfloat(:,1)+1,'.-m'); 
    … …  
    9595 
    9696%legend([Hb Hr],'background','retrieved','Location','SouthEast'); 
    97 legend([Hb Hr Ht],'background.','retrieved','truth'); 
     97legend([Hb Hr Ht],'background.','retrieved','truth','Location','SouthEast'); 
    9898 
    9999%print -dpng ../exp_forw_dan/simu_drifter_dan.png 
    100100 
    101101 
    102  
     102if (false) 
    103103delta_b=sqrt(mean((Ut(:,:,1:end-1)-Ub(:,:,1:end-1)),3).^2+mean((Vt(:,:,1:end-1)-Vb(:,:,1:end-1)),3).^2); 
    104104delta_r=sqrt(mean((Ut(:,:,1:end-1)-Ur(:,:,1:end-1)),3).^2+mean((Vt(:,:,1:end-1)-Vr(:,:,1:end-1)),3).^2); 
    … …  
    155155%print -f1 -dpng ../exp_forw_dan/twin_exp_uv_field.png 
    156156%print -f2 -dpng ../exp_forw_dan/twin_exp_errors.png 
     157end 
     158 
    157159 
    158160%% Erreur le long du trajet 
    … …  
    171173end %id 
    172174end %if false 
     175 
     176%% Divergence 
     177load ../exp_twin_aviso/Div_ana.dat 
     178load ../exp_twin_aviso/Div_fg.dat 
     179 
     180Div2_ana=reshape(Div_ana(:,end),nlat,nlon,ntime-1); 
     181Div2_fg=reshape(Div_fg(:,end),nlat,nlon,ntime-1); 
     182 
     183figure(3) 
     184subplot(3,1,1) 
     185imagesc(lon,lat,mean(Div2_fg,3),[-2 2]*1e-5); 
     186axis xy; 
     187axis(LIM); 
     188title('Mean divergence Background'); 
     189colorbar 
     190 
     191 
     192figure(3) 
     193subplot(3,1,3) 
     194imagesc(lon,lat,mean(Div2_ana,3),[-2 2]*1e-5); 
     195axis xy; 
     196axis(LIM); 
     197title('Mean divergence Analyse'); 
     198colorbar 
     199 
     200subplot(3,1,2); 
     201imagesc(lon,lat,mean(Div2_ana,3)-mean(Div2_fg,3),[-2 2]*1e-5); 
     202axis xy; 
     203axis(LIM); 
     204title('Mean difference'); 
     205 colorbar 

  • altifloat/src/Makefile

    r109 r147  
    1 GENOPT= -E 
     1GENOPT= -DINC_FIL -DINC_GRID -E 
    22YAOPT= +g -x floater_delta 
    33#YAOBIN=${YAODIR}/etc/YaoI 

  • altifloat/src/floater.h

    r144 r147  
    4040//YREAL dtfil=26042000; 
    4141//YREAL dtfil=2.604166666666667e+07; 
     42 
     43#ifdef INC_FIL 
     44#include "fil_def.h" 
     45#else 
    4246YREAL dtfil=200000; 
    4347#endif 
     48#endif 
    4449 
    4550//Options du run incr士ental 
    46 int nb_extiter=4; 
     51int nb_extiter=8; 
    4752short inc_save=0; 
    4853char dirsave[50]="../obs_float/"; 
    … …  
    5560extern int Yobs_insert_data (char *nmmod, int sortie, int iaxe, int jaxe, int kaxe,int pdt, YREAL val); 
    5661void run_inc(int argv, char *argc[]); 
     62void div_init(); 
    5763 
    5864void appli_start(int argc, char *argv[]){ 
    … …  
    8187      pmask[i][j]=0; 
    8288      } 
     89 
     90 
    8391 /* 
    8492char L1[50]="read_obs"; 
    … …  
    274282     Column 1 (%d) : time step (from 0 to jptfl-1) 
    275283     Column 2 (%d) : idfloat (from 0 to jpnfl-1) 
    276      Column 3 (%f) : grid point piret (lon) 
    277      Column 4 (%f) : grid point pjret (lat) 
     284     Column 3 (%f) : grid point piret (lat) 
     285     Column 4 (%f) : grid point pjret (lon) 
    278286  */ 
    279287  if (argv!=2) { 
    … …  
    476484                                              
    477485  fclose(fid); 
     486 
     487} 
     488 
     489void div_init() { 
     490  //Make pseudo-obs of divergence to constrain 
     491  int i,j,it; 
     492  for (it=0;it<jtlag;it++) 
     493    for (j=0;j<nlon;j++) 
     494      for (i=0;i<nlat;i++) { 
     495        Yobs_insert_data("cout_div",0,j,i,0,it,0); 
     496      } 
    478497 
    479498} 
    … …  
    535554  load_init(niti); 
    536555   
     556  //obs pour divergence 
     557  div_init(); 
     558  
    537559  //Réglage de l'ébauche 
    538560   load_eb(niti,nitn); 
    … …  
    563585      // Chargement de toutes les init et obs dans les bon pas de temps 
    564586      load_init(niti); 
     587 
     588   
     589      //obs pour divergence 
     590      div_init(); 
    565591 
    566592      // Réglage de l'ébauche 

  • altifloat/src/floater_delta.d

    r144 r147  
    11#define FILTER 
    22 
    3 defval jtlag 24//nombre de pas de temps Delta t/ delta t  
    4  
    5  
    6  
    7 defval jptfl 25 //Nombre total de pas de temps 
    8 defval jpnfl 3 //nombre total de flotteurs 
     3defval jtlag 48//nombre de pas de temps Delta t/ delta t  
     4 
     5 
     6 
     7defval jptfl 49 //Nombre total de pas de temps 
     8defval jpnfl 5 //nombre total de flotteurs 
    99 
    1010 
    1111//defval nlon 87 //nbre de points de grilles en longitude 
    1212//defval nlat 58 //nbre de points de grilles en latitude 
     13//defval nlon 489 
     14//defval nlat 438 
     15 
     16#ifdef INC_GRID 
     17#include "grid_def.d" 
     18#else 
    1319defval nlon 489 
    1420defval nlat 438 
     21#endif 
    1522 
    1623defval nfloat 5 //nbre de flotteurs max par assim 
    … …  
    1825#ifdef FILTER 
    1926defval K_FILTER //create a define K_FILTER in the source code 
     27 
     28#ifdef INC_FIL 
     29#include "fil_def.d" 
     30#else 
    2031defval OFIL 5  //order of the filter 
     32#endif 
    2133#else 
    2234defval OFIL 0 //Offtime of the main trajectory (0 since there is no filter) 
    … …  
    3648#endif 
    3749 
    38 //    name   M upt offt  dt    nbstep 
    39 traj T_float M 1    OFIL  1    jtlag //trajectoire pour r  
    40 traj T_euler M 0    OFIL  1    jtlag //trajectoire pour u et v 
     50//    name   M     upt   offt     dt    nbstep 
     51traj T_float M      1    OFIL     1     jtlag //trajectoire pour r  
     52traj T_euler M      0    OFIL     1     jtlag //trajectoire pour u et v 
    4153//for varibale u, put under nb step jtlag  
    4254 
    4355//deltas 
    44 traj T_euler_d  D 0    OFIL     1    1  
    45 traj T_float_d D 1    OFIL     1    jtlag //trajectoires type modele ou tangent. permet de controler quel type activer dans le i 
    46  
     56traj T_euler_d  D   0    OFIL     1    1  
     57traj T_float_d  D  1    OFIL     1    jtlag //trajectoires type modele ou tangent. permet de controler quel type activer dans le i 
     58//traj T_euler_d_t D  0    OFIL     1    jtlag 
    4759 
    4860#ifdef FILTER 
    … …  
    5971space S_euler_d M nlon nlat T_euler_d 
    6072space S_eulerlocate_d M nfloat nlon nlat T_float_d 
    61  
     73space S_euler_dt M nlon nlat T_float_d // for div calculation 
    6274 
    6375modul r_float space S_float input 4  output 2 tempo  
    … …  
    6880 
    6981 
    70 modul u       space S_euler noward   output 1  tempo //noward veut dire pas de calcul 
    71 modul v           space S_euler noward   output 1 tempo 
    72 modul locate  space S_eulerlocate input 2  output 1 tempo 
    73 modul mask    space S_euler noward output 1 
    74  
     82modul u       space S_euler       noward    output 1 tempo //noward veut dire pas de calcul 
     83modul v       space S_euler       noward    output 1 tempo 
     84modul locate  space S_eulerlocate input  2  output 1 tempo 
     85modul mask    space S_euler       noward    output 1 
     86modul Div     space S_euler       input  4  output 1 tempo 
     87modul Div_d   space S_euler_d     clonol Div 
    7588 
    7689//deltas 
    … …  
    102115 
    103116modul locate_d  space S_eulerlocate_d input 2  output 1 tempo 
    104 modul r_cout_d space S_float_d input 2  output 2 tempo cout 
     117modul r_cout_d  space S_float_d  input 2  output 2 tempo cout 
     118modul cout_div  space S_euler_dt input 1  output 1 tempo cout 
    105119 
    106120//Filter modules 
    … …  
    113127ctin r_float  1..2   from r_float 1..2 i t-1 // i par default est la premier dimension, ici i=1 floater 
    114128ctin r_float  3..4   from ur      1..2 i t    // comme c'est un calcul, il depend pas de la valeur initiale 
     129 
     130ctin Div      1     from  u       1    i+1  j   t  
     131ctin Div      2     from  u       1    i-1  j   t 
     132ctin Div      3     from  v       1    i    j-1 t 
     133ctin Div      4     from  v       1    i    j+1 t 
     134 
     135ctin Div_d     1     from  u_d     1    i+1 j    
     136ctin Div_d     2     from  u_d     1    i-1 j    
     137ctin Div_d     3     from  v_d     1    i   j-1  
     138ctin Div_d     4     from  v_d     1    i   j+1  
     139 
     140ctin cout_div 1      from  Div_d    1    i   j t 
    115141 
    116142ctin ur       1..8   from uinter 1..8 i t 
    … …  
    171197exec disp_ct_in //afficher les connexions 
    172198 
     199order modinspace S_euler 
     200      order YA1 YA2 
     201            Div 
     202      forder 
     203forder   
     204 
    173205#ifdef FILTER 
    174206order modinspace S_euler_df_t 
    … …  
    180212order modinspace S_euler_d 
    181213        order YA1 YA2 
    182             u_d v_d 
     214            u_d v_d Div_d 
    183215        forder 
    184216forder 
    … …  
    197229forder 
    198230 
     231order spaceintraj T_euler 
     232      S_euler 
     233forder 
    199234 
    200235order spaceintraj T_float 
    … …  
    216251forder 
    217252 
    218  
     253order modinspace S_euler_dt 
     254      order YA1 YA2 
     255            cout_div 
     256      forder 
     257forder 
    219258 
    220259#ifdef FILTER 
    … …  
    229268 
    230269order spaceintraj T_float_d 
    231         S_eulerlocate_d S_float_d 
    232 forder 
    233   
     270        S_eulerlocate_d S_float_d S_euler_dt 
     271forder 
     272  
     273 
    234274  
    235275//insert_fct arg load_aviso 
Note: See TracChangeset for help on using the changeset viewer.