Changeset 26744 for issm/trunk/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp

Timestamp:

12/22/21 10:39:44 (3 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 26742

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/c (modified) (1 prop)
• src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (modified) (2 diffs)

issm/trunk/src/c

@@ -20 +20 @@
 issm
 kriging
+issm_slc
 issm_slr
 issm_ocean

@@ -20 +20 @@
 issm
 kriging
+issm_slc
 issm_slr
 issm_ocean

issm/trunk/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp

@@ -6 +6 @@
 #include "../../shared/shared.h"
 #include "../../toolkits/toolkits.h"
+#include "../../shared/Random/random.h"
 
-void FrontalForcingsx(FemModel* femmodel){
+void FrontalForcingsx(FemModel* femmodel){/*{{{*/
 
         /*Recover melt_parameterization*/
         int melt_parameterization;
         femmodel->parameters->FindParam(&melt_parameterization,FrontalForcingsParamEnum);
-        
+
         /*Calculate melting rate*/
         switch(melt_parameterization){
                 case FrontalForcingsDefaultEnum:
                         break;
+                case FrontalForcingsRignotAutoregressionEnum:
+                        Thermalforcingautoregressionx(femmodel);
+                        /*Do not break here, call IcefrontAreax(),RignotMeltParameterizationx()*/
                 case FrontalForcingsRignotEnum:
                         femmodel->IcefrontAreax();
@@ -24 +28 @@
                         _error_("Frontal forcings "<<EnumToStringx(melt_parameterization)<<" not supported yet");
         }
-}
+}/*}}}*/
+void ThermalforcingautoregressionInitx(FemModel* femmodel){/*{{{*/
+
+   /*Initialization step of Thermalforcingautoregressionx*/
+   int M,N,Nphi,arorder,numbasins,my_rank;
+   IssmDouble starttime,tstep_ar,tinit_ar;
+   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
+   femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
+   IssmDouble* beta0    = NULL;
+   IssmDouble* beta1    = NULL;
+   IssmDouble* phi      = NULL;
+   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+   femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum);
+   femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
+   femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
+
+   /*AR model spin-up with 0 noise to initialize ThermalforcingValuesAutoregressionEnum*/
+        int nspin{2*arorder+5};
+   for(Object* &object:femmodel->elements->objects){
+      Element* element      = xDynamicCast<Element*>(object); //generate element object
+      element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,FrontalForcingsRignotAutoregressionEnum);
+   }
+
+   /*Cleanup*/
+   xDelete<IssmDouble>(beta0);
+   xDelete<IssmDouble>(beta1);
+   xDelete<IssmDouble>(phi);
+}/*}}}*/
+void Thermalforcingautoregressionx(FemModel* femmodel){/*{{{*/
+
+   /*Get time parameters*/
+   IssmDouble time,dt,starttime,tstep_ar;
+   femmodel->parameters->FindParam(&time,TimeEnum);
+   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+   femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum);
+
+   /*Initialize module at first time step*/
+   if(time<=starttime+dt){ThermalforcingautoregressionInitx(femmodel);}
+   /*Determine if this is a time step for the AR model*/
+   bool isstepforar = false;
+
+   #ifndef _HAVE_AD_
+   if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt) isstepforar = true;
+   #else
+   _error_("not implemented yet");
+   #endif
+
+   /*Load parameters*/
+        bool isstochastic;
+   bool istfstochastic = false;
+        int M,N,Nphi,arorder,numbasins,my_rank;
+   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
+   femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
+   IssmDouble tinit_ar;
+   IssmDouble* beta0      = NULL;
+   IssmDouble* beta1      = NULL;
+   IssmDouble* phi        = NULL;
+
+        femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
+   femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
+
+        femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
+        if(isstochastic){
+                int  numstochasticfields;
+                int* stochasticfields;
+                femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
+                femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
+                for(int i=0;i<numstochasticfields;i++){
+                        if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum) istfstochastic = true;
+                }
+                xDelete<int>(stochasticfields);
+        }
+   /*Time elapsed with respect to AR model initial time*/
+   IssmDouble telapsed_ar = time-tinit_ar;
+
+   /*Loop over each element to compute Thermal Forcing at vertices*/
+   for(Object* &object:femmodel->elements->objects){
+      Element* element = xDynamicCast<Element*>(object);
+      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum);
+   }
+
+   /*Cleanup*/
+   xDelete<IssmDouble>(beta0);
+   xDelete<IssmDouble>(beta1);
+   xDelete<IssmDouble>(phi);
+}/*}}}*/

@@ -20 +20 @@
 issm
 kriging
+issm_slc
 issm_slr
 issm_ocean