Changeset 26836 for issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp

Timestamp:

01/29/22 07:39:27 (3 years ago)

Author:

vverjans

Message:

NEW: added class autoregressionlinearbasalforcings

File:

• issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp

@@ -45 +45 @@
                         if(VerboseSolution())_printf0_("        call BeckmannGoosse Floating melting rate module\n");
                         BeckmannGoosseFloatingiceMeltingRatex(femmodel);
+                        break;
+                case AutoregressionLinearFloatingMeltRateEnum:
+                        if(VerboseSolution())_printf0_("        call Autoregression Linear Floating melting rate module\n");
+                        AutoregressionLinearFloatingiceMeltingRatex(femmodel);
                         break;
                 default:
@@ -206 +210 @@
 }
 /*}}}*/
+void AutoregressionLinearFloatingiceMeltingRateInitx(FemModel* femmodel){/*{{{*/
+
+        /*Initialization step of AutoregressionLinearFloatingiceMeltingRatex*/
+   int M,N,Nphi,arorder,numbasins,my_rank;
+   IssmDouble starttime,tstep_ar,tinit_ar;
+   femmodel->parameters->FindParam(&numbasins,BasalforcingsLinearNumBasinsEnum);
+   femmodel->parameters->FindParam(&arorder,BasalforcingsAutoregressiveOrderEnum);
+   IssmDouble* beta0    = NULL;
+   IssmDouble* beta1    = NULL;
+   IssmDouble* phi      = NULL;
+   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+   femmodel->parameters->FindParam(&tstep_ar,BasalforcingsAutoregressionTimestepEnum);
+   femmodel->parameters->FindParam(&tinit_ar,BasalforcingsAutoregressionInitialTimeEnum);
+   femmodel->parameters->FindParam(&beta0,&M,BasalforcingsBeta0Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&beta1,&M,BasalforcingsBeta1Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&phi,&M,&Nphi,BasalforcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
+        
+        /*AR model spin-up with 0 noise to initialize BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum (688 = log(0.001)/log(0.99): decaying time of inluence of phi[0]=0.99 to 0.001 of beta_0*/
+   int nspin = 688;
+   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,BasalforcingsDeepwaterMeltingRateAutoregressionEnum);
+        }
+        /*Cleanup*/
+   xDelete<IssmDouble>(beta0);
+   xDelete<IssmDouble>(beta1);
+   xDelete<IssmDouble>(phi);
+}/*}}}*/
+void AutoregressionLinearFloatingiceMeltingRatex(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,BasalforcingsAutoregressionTimestepEnum);
+
+   /*Initialize module at first time step*/
+   if(time<=starttime+dt){AutoregressionLinearFloatingiceMeltingRateInitx(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 isdeepmeltingstochastic = false;
+   int M,N,Nphi,arorder,numbasins,my_rank;
+   femmodel->parameters->FindParam(&numbasins,BasalforcingsLinearNumBasinsEnum);
+        femmodel->parameters->FindParam(&arorder,BasalforcingsAutoregressiveOrderEnum);
+   IssmDouble tinit_ar;
+   IssmDouble* beta0          = NULL;
+   IssmDouble* beta1          = NULL;
+   IssmDouble* phi            = NULL;
+   IssmDouble* deepwaterel    = NULL;
+   IssmDouble* upperwaterel   = NULL;
+   IssmDouble* upperwatermelt = NULL;
+   IssmDouble* perturbation   = NULL;
+
+        /*Get autoregressive parameters*/
+   femmodel->parameters->FindParam(&tinit_ar,BasalforcingsAutoregressionInitialTimeEnum);
+   femmodel->parameters->FindParam(&beta0,&M,BasalforcingsBeta0Enum);               _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&beta1,&M,BasalforcingsBeta1Enum);               _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&phi,&M,&Nphi,BasalforcingsPhiEnum);             _assert_(M==numbasins); _assert_(Nphi==arorder);
+
+        /*Get basin-specific parameters*/
+   femmodel->parameters->FindParam(&deepwaterel,&M,BasalforcingsDeepwaterElevationEnum);            _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&upperwaterel,&M,BasalforcingsUpperwaterElevationEnum);          _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&upperwatermelt,&M,BasalforcingsUpperwaterMeltingRateEnum);      _assert_(M==numbasins);
+
+        /*Evaluate whether stochasticity on DeepwaterMeltingRate is requested*/
+        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]==BasalforcingsDeepwaterMeltingRateAutoregressionEnum) isdeepmeltingstochastic = 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 FloatingiceMeltingRate at vertices*/
+   for(Object* &object:femmodel->elements->objects){
+      Element* element = xDynamicCast<Element*>(object);
+      /*Compute autoregression*/
+      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,isdeepmeltingstochastic,BasalforcingsDeepwaterMeltingRateAutoregressionEnum);
+                element->BasinLinearFloatingiceMeltingRate(deepwaterel,upperwatermelt,upperwaterel,perturbation);
+        }
+}/*}}}*/