Changeset 25464


Ignore:
Timestamp:
08/25/20 10:59:33 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moving inner loop of transient_core to transient_step

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/cores/transient_core.cpp

    r25340 r25464  
    1717#include "../solutionsequences/solutionsequences.h"
    1818
    19 void transient_core(FemModel* femmodel){
     19/*Prototypes*/
     20void transient_step(FemModel* femmodel);
     21
     22void transient_core(FemModel* femmodel){/*{{{*/
    2023
    2124        /*parameters: */
    2225        IssmDouble finaltime,dt,yts,starttime;
    23         bool       isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isesa,isslr,iscoupler,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,iscontrol,isautodiff;
     26        bool       isoceancoupling,iscontrol,isautodiff,isgroundingline,isslr;
    2427        bool       save_results,dakota_analysis;
    2528        int        timestepping;
    … …  
    2730        int        sb_coupling_frequency;
    2831        int        recording_frequency;
    29         int        domaintype,groundingline_migration,smb_model,amr_frequency,amr_restart;
     32        int        groundingline_migration,smb_model,amr_frequency,amr_restart;
    3033        int        numoutputs;
    31         Analysis  *analysis          = NULL;
    3234        char     **requested_outputs = NULL;
    3335
    … …  
    4143
    4244        /*then recover parameters common to all solutions*/
    43         femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    4445        femmodel->parameters->FindParam(&step,StepEnum);
    4546        femmodel->parameters->FindParam(&time,TimeEnum);
    … …  
    5152        femmodel->parameters->FindParam(&sb_coupling_frequency,SettingsSbCouplingFrequencyEnum);
    5253        femmodel->parameters->FindParam(&timestepping,TimesteppingTypeEnum);
    53         femmodel->parameters->FindParam(&isstressbalance,TransientIsstressbalanceEnum);
    54         femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
    55         femmodel->parameters->FindParam(&issmb,TransientIssmbEnum);
    56         femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
    57         femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
    58         femmodel->parameters->FindParam(&isesa,TransientIsesaEnum);
    5954        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
    60         femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
    6155        femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
    62         femmodel->parameters->FindParam(&ismovingfront,TransientIsmovingfrontEnum);
    6356        femmodel->parameters->FindParam(&isoceancoupling,TransientIsoceancouplingEnum);
    64         femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
    65         femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
    6657        femmodel->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
    67         femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
    6858        femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum);
    6959        femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
    … …  
    128118                femmodel->parameters->SetParam(save_results,SaveResultsEnum);
    129119
    130         #if defined(_HAVE_OCEAN_ )
    131         if(isoceancoupling) OceanExchangeDatax(femmodel,false);
    132         #endif
    133 
    134                 if(isthermal && domaintype==Domain3DEnum){
    135                         if(issmb){
    136                                 bool isenthalpy;
    137                                 femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
    138                                 femmodel->parameters->FindParam(&smb_model,SmbEnum);
    139                                 if(isenthalpy){
    140                                         if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
    141                                                 femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
    142                                                 ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
    143                                         }
    144                                 }
    145                                 else{
    146                                         if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
    147                                                 femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
    148                                                 ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
    149                                         }
    150                                 }
    151                         }
    152                         thermal_core(femmodel);
    153                 }
    154                 /* Using Hydrology dc  coupled we need to compute smb in the hydrology inner time loop*/
    155                 if(issmb) {
    156                         if(VerboseSolution()) _printf0_("   computing smb\n");
    157                         smb_core(femmodel);
    158                 }
    159 
    160                 if(ishydrology){
    161                         if(VerboseSolution()) _printf0_("   computing hydrology\n");
    162                         int hydrology_model;
    163                         hydrology_core(femmodel);
    164                         femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
    165                         if(hydrology_model!=HydrologydcEnum && issmb)smb_core(femmodel);
    166                 }
    167 
    168                 if(isstressbalance && (step%sb_coupling_frequency==0 || step==1) ) {
    169                         if(VerboseSolution()) _printf0_("   computing stress balance\n");
    170                         stressbalance_core(femmodel);
    171                 }
    172 
    173                 if(isdamageevolution) {
    174                         if(VerboseSolution()) _printf0_("   computing damage\n");
    175                         damage_core(femmodel);
    176                 }
    177 
    178                 if(ismovingfront)       {
    179                         if(VerboseSolution()) _printf0_("   computing moving front\n");
    180                         movingfront_core(femmodel);
    181                 }
    182 
    183                 /* from here on, prepare geometry for next time step*/
    184                 //if(issmb) smb_core(femmodel);
    185 
    186                 if(ismasstransport){
    187                         if(VerboseSolution()) _printf0_("   computing mass transport\n");
    188                         bmb_core(femmodel);
    189                         masstransport_core(femmodel);
    190                         femmodel->UpdateVertexPositionsx();
    191                 }
    192 
    193                 if(isgroundingline){
    194 
    195                         /*Start profiler*/
    196                         femmodel->profiler->Start(GROUNDINGLINECORE);
    197 
    198                         if(VerboseSolution()) _printf0_("   computing new grounding line position\n");
    199                         GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    200 
    201                         femmodel->parameters->SetParam(MaskOceanLevelsetEnum,InputToExtrudeEnum);
    202                         extrudefrombase_core(femmodel);
    203                         femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
    204                         extrudefrombase_core(femmodel);
    205                         femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
    206                         extrudefrombase_core(femmodel);
    207 
    208                         /*Stop profiler*/
    209                         femmodel->profiler->Stop(GROUNDINGLINECORE);
    210 
    211                         if(save_results){
    212                                 int outputs[3] = {SurfaceEnum,BaseEnum,MaskOceanLevelsetEnum};
    213                                 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
    214                         }
    215                 }
    216 
    217                 if(isgia){
    218                         if(VerboseSolution()) _printf0_("   computing glacial isostatic adjustment\n");
    219                         #ifdef _HAVE_GIA_
    220                         gia_core(femmodel);
    221                         #else
    222                         _error_("ISSM was not compiled with gia capabilities. Exiting");
    223                         #endif
    224                 }
    225 
    226                 /*esa: */
    227                 if(isesa) esa_core(femmodel);
    228 
    229                 /*Sea level rise: */
    230                 if(isslr) sealevelchange_core(femmodel);
     120                /*Run transient step!*/
     121                transient_step(femmodel);
    231122
    232123                /*unload results*/
    … …  
    263154                }
    264155
    265                 if (iscontrol && isautodiff) {
     156                if(iscontrol && isautodiff){
    266157                        /*Go through our dependent variables, and compute the response:*/
    267158                        for(int i=0;i<dependent_objects->Size();i++){
    … …  
    278169        /*Free ressources:*/
    279170        if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
    280 }
     171
     172}/*}}}*/
     173void transient_step(FemModel* femmodel){/*{{{*/
     174
     175        /*parameters: */
     176        bool isstressbalance,ismasstransport,issmb,isthermal,isgroundingline,isgia,isesa;
     177        bool isslr,ismovingfront,isdamageevolution,ishydrology,isoceancoupling;
     178        bool save_results;
     179        int  step,sb_coupling_frequency;
     180        int  domaintype,smb_model;
     181
     182        /*then recover parameters common to all solutions*/
     183        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
     184        femmodel->parameters->FindParam(&step,StepEnum);
     185        femmodel->parameters->FindParam(&sb_coupling_frequency,SettingsSbCouplingFrequencyEnum);
     186        femmodel->parameters->FindParam(&isstressbalance,TransientIsstressbalanceEnum);
     187        femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
     188        femmodel->parameters->FindParam(&issmb,TransientIssmbEnum);
     189        femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
     190        femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
     191        femmodel->parameters->FindParam(&isesa,TransientIsesaEnum);
     192        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
     193        femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
     194        femmodel->parameters->FindParam(&ismovingfront,TransientIsmovingfrontEnum);
     195        femmodel->parameters->FindParam(&isoceancoupling,TransientIsoceancouplingEnum);
     196        femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
     197        femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
     198
     199        #if defined(_HAVE_OCEAN_ )
     200        if(isoceancoupling) OceanExchangeDatax(femmodel,false);
     201        #endif
     202
     203        if(isthermal && domaintype==Domain3DEnum){
     204                if(issmb){
     205                        bool isenthalpy;
     206                        femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
     207                        femmodel->parameters->FindParam(&smb_model,SmbEnum);
     208                        if(isenthalpy){
     209                                if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
     210                                        femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
     211                                        ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
     212                                }
     213                        }
     214                        else{
     215                                if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
     216                                        femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
     217                                        ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
     218                                }
     219                        }
     220                }
     221                thermal_core(femmodel);
     222        }
     223        /* Using Hydrology dc  coupled we need to compute smb in the hydrology inner time loop*/
     224        if(issmb) {
     225                if(VerboseSolution()) _printf0_("   computing smb\n");
     226                smb_core(femmodel);
     227        }
     228
     229        if(ishydrology){
     230                if(VerboseSolution()) _printf0_("   computing hydrology\n");
     231                int hydrology_model;
     232                hydrology_core(femmodel);
     233                femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
     234                if(hydrology_model!=HydrologydcEnum && issmb)smb_core(femmodel);
     235        }
     236
     237        if(isstressbalance && (step%sb_coupling_frequency==0 || step==1) ) {
     238                if(VerboseSolution()) _printf0_("   computing stress balance\n");
     239                stressbalance_core(femmodel);
     240        }
     241
     242        if(isdamageevolution) {
     243                if(VerboseSolution()) _printf0_("   computing damage\n");
     244                damage_core(femmodel);
     245        }
     246
     247        if(ismovingfront)       {
     248                if(VerboseSolution()) _printf0_("   computing moving front\n");
     249                movingfront_core(femmodel);
     250        }
     251
     252        /* from here on, prepare geometry for next time step*/
     253
     254        if(ismasstransport){
     255                if(VerboseSolution()) _printf0_("   computing mass transport\n");
     256                bmb_core(femmodel);
     257                masstransport_core(femmodel);
     258                femmodel->UpdateVertexPositionsx();
     259        }
     260
     261        if(isgroundingline){
     262
     263                /*Start profiler*/
     264                femmodel->profiler->Start(GROUNDINGLINECORE);
     265
     266                if(VerboseSolution()) _printf0_("   computing new grounding line position\n");
     267                GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     268
     269                femmodel->parameters->SetParam(MaskOceanLevelsetEnum,InputToExtrudeEnum);
     270                extrudefrombase_core(femmodel);
     271                femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
     272                extrudefrombase_core(femmodel);
     273                femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
     274                extrudefrombase_core(femmodel);
     275
     276                /*Stop profiler*/
     277                femmodel->profiler->Stop(GROUNDINGLINECORE);
     278
     279                if(save_results){
     280                        int outputs[3] = {SurfaceEnum,BaseEnum,MaskOceanLevelsetEnum};
     281                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
     282                }
     283        }
     284
     285        if(isgia){
     286                if(VerboseSolution()) _printf0_("   computing glacial isostatic adjustment\n");
     287                #ifdef _HAVE_GIA_
     288                gia_core(femmodel);
     289                #else
     290                _error_("ISSM was not compiled with gia capabilities. Exiting");
     291                #endif
     292        }
     293
     294        /*esa: */
     295        if(isesa) esa_core(femmodel);
     296
     297        /*Sea level rise: */
     298        if(isslr) sealevelchange_core(femmodel);
     299
     300}/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.