Changeset 26744 for issm/trunk/src/c/analyses/StressbalanceVerticalAnalysis.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/analyses/StressbalanceVerticalAnalysis.cpp (modified) (6 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/analyses/StressbalanceVerticalAnalysis.cpp ¶

@@ -5 +5 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+
+//#define INWOOVZ 1 //For Inwoo Park 0: incompressible assumption(IA), 1: internal deformation(ID), 2: IA+ID
 
 /*Model processing*/
@@ -105 +107 @@
         iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
         iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
-        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
+        iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0);
         iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
         if(iomodel->domaintype!=Domain2DhorizontalEnum){
@@ -112 +114 @@
         }
         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
-        //iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
-
+        #ifdef INWOOVZ
+        iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
+        #endif
 
         /*Add basal forcings to compute melt rate*/
+        bool isstochastic;
         int basalforcing_model;
         iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
+        iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing");
         switch(basalforcing_model){
                 case FloatingMeltRateEnum:
                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
+                        if(isstochastic){
+            iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
+            iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.floatingice_melting_rate",BaselineBasalforcingsFloatingiceMeltingRateEnum);
+         }
                         break;
                 case LinearFloatingMeltRateEnum:
@@ -153 +162 @@
 void StressbalanceVerticalAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 
-        /*No specific parameters*/
+        /* specific parameters*/
 
 }/*}}}*/
@@ -163 +172 @@
                 femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum);
                 solutionsequence_linear(femmodel);
+}/*}}}*/
+void           StressbalanceVerticalAnalysis::PreCore(FemModel* femmodel){/*{{{*/
+        _error_("not implemented");
 }/*}}}*/
 ElementVector* StressbalanceVerticalAnalysis::CreateDVector(Element* element){/*{{{*/
@@ -576 +588 @@
         }
 
+        #ifdef INWOOVZ
+        IssmDouble*  thickness = xNew<IssmDouble>(numnodes);
+        IssmDouble*  base            = xNew<IssmDouble>(numnodes);
+        IssmDouble*  smb       = xNew<IssmDouble>(numnodes);
+
+        IssmDouble rheology_n=element->material->GetN();
+        IssmDouble isvelbc   = 10./(365.25*24*3600); /*10 m/yr*/
+        /*Set analytical vertical velocity field at slow flow region.*/
+        if(INWOOVZ == 1){
+                element->GetInputListOnNodes(&thickness[0],ThicknessEnum,0.);
+                element->GetInputListOnNodes(&base[0],BaseEnum,0.);
+                element->GetInputListOnNodes(&smb[0],SmbMassBalanceEnum,0.);
+                for(i=0;i<numnodes;i++){
+                        IssmDouble eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
+                        vz[i] = -(pow(eta,rheology_n+2.0)-1-(rheology_n+2.0)*(eta-1.))/(rheology_n+1.)*smb[i];
+                }       
+        }
+        else if(INWOOVZ== 2) {
+                element->GetInputListOnNodes(&thickness[0],ThicknessEnum,0.);
+                element->GetInputListOnNodes(&base[0],BaseEnum,0.);
+                element->GetInputListOnNodes(&smb[0],SmbMassBalanceEnum,0.);
+                for(i=0;i<numnodes;i++){
+                        if ( vx[i]*vx[i]+vy[i]*vy[i] < isvelbc*isvelbc ){
+                                IssmDouble eta = (base[i]+thickness[i]-xyz_list[i*3+2])/thickness[i];
+                                vz[i] = -(pow(eta,rheology_n+2.0)-1-(rheology_n+2.0)*(eta-1.))/(rheology_n+1.)*smb[i];
+                        }
+                }       
+        }
+        else if(INWOOVZ==0){
+                /*nothing to add*/
+        }
+        else{
+                _error_("not supported yet");
+        }
+
+        /*Cleanup*/
+        xDelete<IssmDouble>(smb);
+        xDelete<IssmDouble>(base);
+        xDelete<IssmDouble>(thickness);
+        #endif
+
         /*Now Compute vel*/
         for(i=0;i<numnodes;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);

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