Changes between Version 21 and Version 22 of 2020WP/KERNEL-06_techene_better_e3_management


Ignore:
Timestamp:
2020-10-26T18:32:43+01:00 (6 months ago)
Author:
techene
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • 2020WP/KERNEL-06_techene_better_e3_management

    v21 v22  
    4545'''KERNEL-06's version 1 [pre mid-merge 2020] : /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3'''  
    4646 
    47 ==== starting point 
     47When key_qco is not activated NEMO should be produce exactly same results as the trunk and passes SETTE version r13167 has been delivered for mid-merge party ! 
    4848 
    49 In NEMO r12377 scale factors (e3*) at u-v-w-uw-vw-f-points are interpolated from e3t at Kbb, Kmm, and Kaa. 
    50 The module in charge of scale factor management is src/OCE/DOM/domvvl.F90.  
    51   
    52 domvvl.F90 interpolation routine is called :  
    53 - at initialisation or restart at u-v-w-uw-vw-f-points 
    54 {{{ 
    55 nemogcm  
    56 --> nemo_init 
    57    --> dom_init 
    58       --> dom_vvl_init 
    59          --> dom_vvl_rst 
    60          --> dom_vvl_zgr 
    61             --> dom_vvl_interpol  
    62 }}} 
    63 - at each time step for "after" scale factor at u-v-points each time ssh[Naa] is computed 
    64 {{{ 
    65 nemogcm  
    66 --> stp 
    67    --> dom_vvl_sf_nxt 
    68       --> e3t[Kaa] 
    69       --> dom_vvl_interpol 
    70 }}} 
    71 - at each time step for "now" scale factor at u-v-points e3t is directly filtered 
    72 {{{ 
    73 nemogcm  
    74 --> stp 
    75    --> tra_atf 
    76       --> e3t[Kmm] 
    77    --> dyn_atf 
    78       --> e3t[Kmm] 
    79       --> dom_vvl_interpol 
    80 }}} 
    81 - at each time step after time swapping for "now" at f-point and "before" scale factor both at w-uw-vw-points 
    82 {{{ 
    83 nemogcm  
    84 --> stp 
    85    --> dom_vvl_sf_update 
    86       --> dom_vvl_interpol 
    87 }}} 
    88  
    89 ==== developments  
    90  
    91 In version 1 of source:/NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3 we implement changes progressively and validate step by step regarding GYRE_PISCES test case. We remove all the vvl routines usage thus we replace e3t/u/v/w/uw/vw at 3 time step + e3f (19) 3d tables storage and twice e3u/v "after" + e3u/v "now" + e3f/w/uw/vw "now" + e3w/uw/vw "before" (13) 3d tables interpolation at each step and e3u/v/f/w/uw/vw "now" + e3u/v/w/uw/vw "before" (11) 3d tables interpolation at initialisation or restart by r3t/u/v at 3 time steps and r3f i.e. (10) 2d table storage and on the fly light computation. For backward compatibility we introduce a cpp key key_qco in order to isolate new scale factor implementation from former vvl version. qco stand for for "quasi eulerian coordinate".  
    92 - new variables added in dom_oce and domain  
    93 {{{ 
    94 r3P with P = [t,u,v,f] are 2d in space ssh/hP0 
    95 }}} 
    96 - new module in DOM : dom_qco 
    97 {{{ 
    98 -- dom_qco_r3c ! compute ssh/ht0 at t-point and interpolate ssh/h.0 at u-v-f-points from ssh 
    99 -- dom_qco_zgr ! set ssh/ht0 at t-u-v-f-points for Kmm and at t-u-v-points for Kbb 
    100    --> dom_qco_r3c[Kmm] at t-u-v-f-points 
    101    --> dom_qco_r3c[Kbb] at t-u-v  -points 
    102 }}} 
    103 - new substitute in DOM : domzgr_substitute When key_qco is active e3. is no longer a variable but an expression. Each time e3. appears in a routine an include domzgr_substitute in the module enables to replace it by a (e3._0 ( 1 + r3. ) * mask.) like expression.  
    104 {{{ 
    105 #   define  e3t(i,j,k,t)   (e3t_0(i,j,k)*(1._wp+r3t(i,j,t)*tmask(i,j,k))) 
    106 #   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1._wp+r3u(i,j,t)*umask(i,j,k))) 
    107 #   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1._wp+r3v(i,j,t)*vmask(i,j,k))) 
    108 #   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1._wp+r3f(i,j)*fmask(i,j,k))) 
    109 #   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1._wp+r3t(i,j,t))) 
    110 #   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1._wp+r3u(i,j,t))) 
    111 #   define  e3vw(i,j,k,t)  (e3vw_0(i,j,k)*(1._wp+r3v(i,j,t))) 
    112 }}} 
    113 - finally we extended the e3. modification to h. r1_h. and gde.. 
    114 {{{ 
    115 #   define  ht(i,j)        (ht_0(i,j)+ssh(i,j,Kmm)) 
    116 #   define  hu(i,j,t)      (hu_0(i,j)*(1._wp+r3u(i,j,t))) 
    117 #   define  hv(i,j,t)      (hv_0(i,j)*(1._wp+r3v(i,j,t))) 
    118 #   define  r1_hu(i,j,t)   (r1_hu_0(i,j)/(1._wp+r3u(i,j,t))) 
    119 #   define  r1_hv(i,j,t)   (r1_hv_0(i,j)/(1._wp+r3v(i,j,t))) 
    120 #   define  gdept(i,j,k,t) (gdept_0(i,j,k)*(1._wp+r3t(i,j,t))) 
    121 #   define  gdepw(i,j,k,t) (gdepw_0(i,j,k)*(1._wp+r3t(i,j,t))) 
    122 #   define  gde3w(i,j,k)   (gdept_0(i,j,k)*(1._wp+r3t(i,j,Kmm))-ssh(i,j,Kmm)) 
    123 }}} 
    124  
    125 When this key_qco is not activated NEMO should be exactly the same as the trunk.  
    126  
    127  
    128 Important points :  
     49Modification of the code have been implemented step by step. This enabled to stress out some interesting features.  
    12950- e3. expression involves tables of distinct dimension then e3.(:,:,: ) call fails it may be necessary to introduce temporary variables (same for water height expression) 
    13051- "e3. =" is no longer possible 
    13152- e3t/u/v/f modifications did not introduce any difference in the results, e3w modification does because both approaches vvl and qco do not take into account of the bottom level in the same way 
    132 - in GYRE e3w_0 are not the half sum of e3u_0 so the way it is implemented in the reference version is not convinient 
     53- in GYRE e3w_0 are not the half sum of e3u_0 so the way it is implemented in the reference version is not convenient 
    13354- e3. substitution makes lines longer than 136 character this may be a problem for compilers (most have been checked but not all) 
    13455- ssh filtering has been displaced upper in order to provide filtered r3P in TOP asselin filtering  
    135 - when key_qco is not active we pass SETTE and this version r13167 has been delivered for mid-merge party ! Some silly allocating memory bugs found and a not that silly bug in implicit mode triggering for SPITZ12 configuration (Dt instead of 2Dt required). 
    13656 
    13757 
     
    13959'''KERNEL-06's version 2 [post mid-merge 2020] : dev_r13327_KERNEL-06_2_techene_e3''' 
    14060 
    141 In this branch we debug with respect to sette tests :  
    142 - non linear version of GYRE_PISCES shows restartability issue due to  
    143 BUG 1 : dom_qco_init.F90 dom_init is called before dyn_adv_init, then default value of ln_dynadv_vec is at false. r3P coefficients computation is not the same if we use flux of vector form. GYRE_PISCES runs in vector form so at the restart r3P are initialised with their flux form instead of their vector from. 
    144 RESOLUTION : cpp key added key_qcoTest_FluxForm performance tests have to be done if performances are not relevant the flux from qco switch will be removed. Note that a switch as been added in dynspg_ts to compute surface average consistently with qco flux form. En gros soit on fait la distinction a chaque interpolation lineaire de la ssh soit pas du tout, il faut juste être consistent. Si les performances flux form (il y a moins de calcul) ne sont pas significatives, on pourra abandonner le calcul sinon il faudra ajouter une option de restart. Routines from domqco.F90 and dynapg_ts.F90 havs to be modified. 
     61When key_qco is activated NEMO passes suitable SETTE tests i.e. custom GYRE_PISCES, ORCA2_ICE_PISCES, SPITZ12, AMM12. ISOMIP cannot run yet nor AGRIF.  
    14562 
    146 BUG 2.1 : in dynatf_qco.F90 barotrope uub and vvb are compute with the un-filtrered scale factor 
    147 BUG 2.2 : the way to compute uub and vvb with e3 from ssh/h0 is not exactly the same as the implementation of e3 filtered  
    148 RESOLUTION : use filtered r3P_f and be carefull the computation has to be done exactly as in domzgr_substitute i.e. including mask and parenthesis otherwise the restartability fails in dynatf_qco.F90. 
    149  
    150 BUG 3: when time splitting is not active ssh/h0 at f-points after the dynamic was not computed  
    151 RESOLUTION : an extra dom_qco_r3c after dyn_spg_ts in case ln_dynspg_ts = F in stpMLF.F90 
    152  
    153 - ORCA2_ICE_PISCES 
    154 WARNING in debug mode restartability fails it is the case even for the trunk @ 13327 need to be confirmed running tests on the current trunk when it will run...  
    155  
    156 - SPITZ12 OK 
    157  
    158 - AMM12 OK 
     63Pseudo merge with NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater consists mainly in modifying dynldf_lap_blp for adding symmetric tensor operator and dynvor for cleaning and diawri for F-point outputs and adding a test case with specific RK3 step.  
    15964 
    16065 
    161 ''...'' 
    16266 
    16367=== Documentation updates 
     
    16872}}} 
    16973 
    170 ''...'' 
     74Need for documentation on qco changes !  
     75Need for documentation on symetric tensor for viscosity operator 
     76Need for documentation on e3f computation in ENE and ENS in zps ? + dynspg_ts ! 
    17177 
    17278== Preview  
     
    17884 
    17985 
    180  
    181 Eventually, all the dom_vvl_interpol call are removed, each time e3 is called we use a substitute to replace e3 by e3_0 (1 + ssh / h_0). For backward compatibility a cpp key manages the use of the new version vs. the old version. We will duplicate modules such as step and domvvl into stepLF and domQE (QE stands for Quasi Eulerian) and create a subtitute module.  
    182  
    183 Integrated in mid merge trunk.  
    184  
    185 List the Fortran modules and subroutines to be created. 
    186 substitute.F90 
    187  
    188 Step 1 : Check the error for e3t, e3w between the current way to compute e3 at T-, W-point and the proposed way to compute e3 at T-, W-point. 
    189 - prints added with no change in the results 
    190 Step 2 : First we change only the core routine in domvvl which should be changed into domQE. 
    191 - add new variables, duplicate step into steplf and domvvl into domQE 
    192 - change interpolation routines into scaling routines in domQE 
    193 Step 3 : Then we change the Asselin filtering routine indeed because water forcing are applied locally. 
    194 - change Asselin routines (maybe not required since e3 scale with vertical with JC modif) 
    195 Step 4 : Finally we remove the interpol routine in the whole code 
    196 - remove interpolating routine in all the code (AGRIF, OFF,...)  
    197 - use a SUBSTITUTE when there are e3 CALL 
    198 - make some changes in step and domQE to have the whole thing consistent 
    199  
    20086''...'' 
    20187 
    202 == Tests 
    20388 
    204 {{{#!box width=50em info 
    205 [[Include(wiki:Developers/DevProcess#tests)]] 
    206 }}} 
    207  
    208 ''...'' 
    209 We want to track and maybe explain the differences observed at every steps.  
    210 Reference set up : For that we produce a reference data set with the trunk -r 12377 using the GYRE_PISCES configuration where top cpp_key has been removed. We run it on 120 time steps. The drag coefficient is zero. We XIOS output an averaged field every 5 days.  
    211  
    212 Step 1 : We print MAXVAL of error between both way to compute the vertical scale factors at each time step, note that we cancelled forcing (in the r12377 revision it should not change anything since water forcings such as run off and emp scale with the vertical). 
    213  
    214 error between proposed and former way to compute vertical scale factors at time kt = 1, 120, 85 
    215 e3t (1) 
    216 0.0000000000000000 
    217 3999.6591076268369 
    218 4.54747350886E-013 
    219 e3t (2) 
    220 5.68434188608E-014 
    221 5.11590769747E-013 
    222 4.54747350886E-013 
    223 e3w 
    224 '''4.64477238892E-007 
    225 6.13657050507E-006 
    226 5.27333801869E-006''' 
    227 gde3w 
    228 1.81898940354E-012 
    229 2.72848410531E-012 
    230 2.72848410531E-012 
    231  
    232 QUESTION : Why do we have such an error on the e3w scale factors ? It is not consistent with machine accuracy error. It seems to be related to the e3w_0 computation. How do we compute e3w_0 ?  
    233 OK SOLVED ! THIS IS A KIND OF ERROR IN THE CODE !!! DUE TO THE FACT THAT E3W_0(jk) != 0.5 * ( E3T_0(jk) + E3T_0(jk-1) )... 
    234  
    235  
    236 Step 2 : Change the code in domvvl turn into domqe. Duplicate step.F90 into steplf.F90 and call domqe routines inside. 
    237  
    238 We observe small errors but not errors at the truncature level as expected with the curent trunk version. This is due to the differences spotted above.  
    239 WE CAN NO LONGER USE THE TRUNK PRODUCTION AS A REFERENCE...  
    24089 
    24190