transient_core.cpp File Reference

#include <float.h>
#include "./cores.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"

Go to the source code of this file.

Functions
void	transient_core (FemModel *femmodel)

Function Documentation

transient_core()

void transient_core

(

FemModel *

femmodel

)

Definition at line 19 of file transient_core.cpp.

                                       {
 
   /*parameters: */
   IssmDouble finaltime,dt,yts,starttime;
   bool       isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isesa,isslr,iscoupler,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,iscontrol,isautodiff;
   bool       save_results,dakota_analysis;
   int        timestepping;
   int        output_frequency;
   int        sb_coupling_frequency;
   int        recording_frequency;
   int        domaintype,groundingline_migration,smb_model,amr_frequency,amr_restart;
   int        numoutputs;
   Analysis  *analysis          = NULL;
   char     **requested_outputs = NULL;
 
   /*intermediary: */
   int        step;
   IssmDouble time;
 
   /*first, figure out if there was a check point, if so, do a reset of the FemModel* femmodel structure. */
   femmodel->parameters->FindParam(&recording_frequency,SettingsRecordingFrequencyEnum);
   if(recording_frequency) femmodel->Restart();
 
   /*then recover parameters common to all solutions*/
   femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
   femmodel->parameters->FindParam(&step,StepEnum);
   femmodel->parameters->FindParam(&time,TimeEnum);
   femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
   femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
   femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
   femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum);
   femmodel->parameters->FindParam(&sb_coupling_frequency,SettingsSbCouplingFrequencyEnum);
   femmodel->parameters->FindParam(&timestepping,TimesteppingTypeEnum);
   femmodel->parameters->FindParam(&isstressbalance,TransientIsstressbalanceEnum);
   femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
   femmodel->parameters->FindParam(&issmb,TransientIssmbEnum);
   femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
   femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
   femmodel->parameters->FindParam(&isesa,TransientIsesaEnum);
   femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
   femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
   femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
   femmodel->parameters->FindParam(&ismovingfront,TransientIsmovingfrontEnum);
   femmodel->parameters->FindParam(&isoceancoupling,TransientIsoceancouplingEnum);
   femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
   femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
   femmodel->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
   femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
   femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum);
   femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
   if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
   femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
   if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum);
 
   #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
   if(amr_frequency){
     femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum);
     if(amr_restart) femmodel->ReMesh();
  }
   #endif
 
   #if defined(_HAVE_OCEAN_ )
   if(isoceancoupling) OceanExchangeDatax(femmodel,true);
   #endif
 
   IssmDouble  output_value;
   int         num_dependents;
   IssmPDouble *dependents;
   DataSet*    dependent_objects=NULL;
   IssmDouble  J=0.;
 
   if(iscontrol && isautodiff){
 
      femmodel->parameters->SetParam(starttime,TimeEnum);
      femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
      femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
 
      /*Go through our dependent variables, and compute the response:*/
      dependents=xNew<IssmPDouble>(num_dependents);
   }
 
   if(isslr) sealevelrise_core_geometry(femmodel);
 
   while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime.
 
      /*Time Increment*/
      switch(timestepping){
         case AdaptiveTimesteppingEnum:
            femmodel->TimeAdaptx(&dt);
            if(time+dt>finaltime) dt=finaltime-time;
            femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
            break;
         case FixedTimesteppingEnum:
            femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
            break;
         default:
            _error_("Time stepping \""<<EnumToStringx(timestepping)<<"\" not supported yet");
      }
      step+=1;
      time+=dt;
      femmodel->parameters->SetParam(time,TimeEnum);
      femmodel->parameters->SetParam(step,StepEnum);
 
      if(VerboseSolution()) _printf0_("iteration " << step << "/" << ceil((finaltime-time)/dt)+step << "  time [yr]: " << setprecision(4) << time/yts << " (time step: " << dt/yts << ")\n");
      if(step%output_frequency==0 || (time >= finaltime - (yts*DBL_EPSILON)) || step==1)
       save_results=true;
      else
       save_results=false;
      femmodel->parameters->SetParam(save_results,SaveResultsEnum);
 
   #if defined(_HAVE_OCEAN_ )
   if(isoceancoupling) OceanExchangeDatax(femmodel,false);
   #endif
 
      if(isthermal && domaintype==Domain3DEnum){
         if(issmb){
            bool isenthalpy;
            femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
            femmodel->parameters->FindParam(&smb_model,SmbEnum);
            if(isenthalpy){
               if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
                  femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
                  ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
               }
            }
            else{
               if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
                  femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
                  ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
               }
            }
         }
         thermal_core(femmodel);
      }
      /* Using Hydrology dc  coupled we need to compute smb in the hydrology inner time loop*/
      if(issmb) {
         if(VerboseSolution()) _printf0_("   computing smb\n");
         smb_core(femmodel);
      }
 
      if(ishydrology){
         if(VerboseSolution()) _printf0_("   computing hydrology\n");
         int hydrology_model;
         hydrology_core(femmodel);
         femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
         if(hydrology_model!=HydrologydcEnum && issmb)smb_core(femmodel);
      }
 
      if(isstressbalance && (step%sb_coupling_frequency==0 || step==1) ) {
         if(VerboseSolution()) _printf0_("   computing stress balance\n");
         stressbalance_core(femmodel);
      }
 
      if(isdamageevolution) {
         if(VerboseSolution()) _printf0_("   computing damage\n");
         damage_core(femmodel);
      }
 
      if(ismovingfront) {
         if(VerboseSolution()) _printf0_("   computing moving front\n");
         movingfront_core(femmodel);
      }
 
      /* from here on, prepare geometry for next time step*/
      //if(issmb) smb_core(femmodel);
 
      if(ismasstransport){
         if(VerboseSolution()) _printf0_("   computing mass transport\n");
         bmb_core(femmodel);
         masstransport_core(femmodel);
         femmodel->UpdateVertexPositionsx();
      }
 
      if(isgroundingline){
 
         /*Start profiler*/
         femmodel->profiler->Start(GROUNDINGLINECORE);
 
         if(VerboseSolution()) _printf0_("   computing new grounding line position\n");
         GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 
         femmodel->parameters->SetParam(MaskOceanLevelsetEnum,InputToExtrudeEnum);
         extrudefrombase_core(femmodel);
         femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
         extrudefrombase_core(femmodel);
         femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
         extrudefrombase_core(femmodel);
 
         /*Stop profiler*/
         femmodel->profiler->Stop(GROUNDINGLINECORE);
 
         if(save_results){
            int outputs[3] = {SurfaceEnum,BaseEnum,MaskOceanLevelsetEnum};
            femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
         }
      }
 
      if(isgia){
         if(VerboseSolution()) _printf0_("   computing glacial isostatic adjustment\n");
         #ifdef _HAVE_GIA_
         gia_core(femmodel);
         #else
         _error_("ISSM was not compiled with gia capabilities. Exiting");
         #endif
      }
 
      /*esa: */
      if(isesa) esa_core(femmodel);
 
      /*Sea level rise: */
      if(isslr) sealevelchange_core(femmodel);
 
      /*unload results*/
      if(VerboseSolution()) _printf0_("   computing requested outputs\n");
      femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs,save_results);
      if(isgroundingline && (groundingline_migration==AggressiveMigrationEnum || groundingline_migration==ContactEnum)){
         int outputs[1] = {MaskOceanLevelsetEnum};
         femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1,save_results);
      }
 
      if(save_results){
         if(VerboseSolution()) _printf0_("   saving temporary results\n");
         OutputResultsx(femmodel);
      }
 
      if(recording_frequency && (step%recording_frequency==0)){
         if(VerboseSolution()) _printf0_("   checkpointing model \n");
         femmodel->CheckPoint();
      }
 
      /*Adaptive mesh refinement*/
      if(amr_frequency){
 
         #if !defined(_HAVE_AD_)
         if(save_results) femmodel->WriteMeshInResults();
         if(step%amr_frequency==0 && time<finaltime){
            if(VerboseSolution()) _printf0_("   refining mesh\n");
            femmodel->ReMesh();//Do not refine the last step
         }
 
         #else
         _error_("AMR not suppored with AD");
         #endif
      }
 
      if (iscontrol && isautodiff) {
         /*Go through our dependent variables, and compute the response:*/
         for(int i=0;i<dependent_objects->Size();i++){
            DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
             dep->Responsex(&output_value,femmodel);
             dep->AddValue(output_value);
            }
         }
   }
 
   if(!iscontrol || !isautodiff) femmodel->RequestedDependentsx();
   if(iscontrol && isautodiff) femmodel->parameters->SetParam(dependent_objects,AutodiffDependentObjectsEnum);
 
   /*Free ressources:*/
   if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
}