Changeset 25941

Timestamp:

01/18/21 03:37:32 (4 years ago)

Author:

inwoo

Message:

CHG: New feature controls vertical velocity options, such as incompressible assumption, internal deformation, and combination of both assumption.

Location:

issm/trunk-jpl/src

Files:

• c/analyses/StressbalanceVerticalAnalysis.cpp (modified) (7 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• m/classes/stressbalance.m (modified) (5 diffs)

issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp

@@ -112 +112 @@
         }
         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
-        //iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
+        iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
 
 
@@ -153 +153 @@
 void StressbalanceVerticalAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 
-        /*No specific parameters*/
+        /* specific parameters*/
+        parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.isvzsolution",StressbalanceIsVzSolutionEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.isvelbc",StressbalanceIsVelBCEnum));
 
 }/*}}}*/
@@ -516 +518 @@
         IssmDouble   rho_ice,g;
 
+        int          isvzsolution;
+        IssmDouble   isvelbc, rheology_n, eta;
+
         /*Get the approximation and do nothing if the element in FS or None*/
         element->GetInputValue(&approximation,ApproximationEnum);
@@ -521 +526 @@
                 return;
         }
+
+        /*Get parameters*/
+        element->FindParam(&isvzsolution,StressbalanceIsVzSolutionEnum);
+        element->FindParam(&isvelbc,StressbalanceIsVelBCEnum);
 
         /*Get dof list and vertices coordinates: */
@@ -535 +544 @@
         IssmDouble*  pressure  = xNew<IssmDouble>(numnodes);
         IssmDouble*  surface   = xNew<IssmDouble>(numnodes);
+        IssmDouble*  thickness = xNew<IssmDouble>(numnodes);
+        IssmDouble*  base            = xNew<IssmDouble>(numnodes);
+        IssmDouble*  smb       = xNew<IssmDouble>(numnodes);
 
         /*Use the dof list to index into the solution vector vz: */
@@ -575 +587 @@
                 }
         }
+
+        /*Set analytical vertical velocity field at slow flow region.*/
+        if (isvzsolution == 1){
+                rheology_n=element->material->GetN();
+                element->GetInputListOnNodes(&thickness[0],ThicknessEnum,0.);
+                element->GetInputListOnNodes(&base[0],BaseEnum,0.);
+                element->GetInputListOnNodes(&smb[0],SmbMassBalanceEnum,0.);
+                for(i=0;i<numnodes;i++){
+                        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 (isvzsolution == 2) {
+                rheology_n=element->material->GetN();
+                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 ){
+                                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];
+                        }
+                }       
+        }
+
 
         /*Now Compute vel*/
@@ -600 +637 @@
 
         /*Free ressources:*/
+        xDelete<IssmDouble>(smb);
+        xDelete<IssmDouble>(base);
+        xDelete<IssmDouble>(thickness);
         xDelete<IssmDouble>(surface);
         xDelete<IssmDouble>(pressure);

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -441 +441 @@
         StressbalanceRiftPenaltyThresholdEnum,
         StressbalanceShelfDampeningEnum,
+        StressbalanceIsVzSolutionEnum,
+        StressbalanceIsVelBCEnum,
         ThermalIsdrainicecolumnEnum,
         ThermalIsdynamicbasalspcEnum,

issm/trunk-jpl/src/m/classes/stressbalance.m

@@ -23 +23 @@
                 loadingforce           = NaN;
                 requested_outputs      = {};
+                isvzsolution           = NaN;
+                isvelbc                = NaN;   
         end
         methods
@@ -75 +77 @@
                          self.rift_penalty_lock=10;
 
+                         %Set vertical velocity option
+                         self.isvzsolution = 0;
+                         self.isvelbc      = 10.;
+
                          %output default:
                          self.requested_outputs={'default'};
@@ -96 +102 @@
                         md = checkfield(md,'fieldname','stressbalance.referential','size',[md.mesh.numberofvertices 6]);
                         md = checkfield(md,'fieldname','stressbalance.loadingforce','size',[md.mesh.numberofvertices 3]);
+                        md = checkfield(md,'fieldname','stressbalance.isvzsolution','numel',[1],'values',[0 1 2]);
+                        md = checkfield(md,'fieldname','stressbalance.isvelbc','numel',[1],'>=',0.0);
                         md = checkfield(md,'fieldname','stressbalance.requested_outputs','stringrow',1);
                         if ~any(isnan(md.stressbalance.vertex_pairing)),
@@ -158 +166 @@
                         fielddisplay(self,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset');
                         fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized');
+
+                        disp(sprintf('\n      %s','Estimation of vertical velocity field options:'));
+                        fielddisplay(self,'isvzsolution','0: incompressible assumption(IA), 1: internal deformation(ID), 2: IA+ID');
+                        fielddisplay(self,'isvelbc','If turned on isvzsolution as 2, region (vel < isvelbc) calculates for ID.');
 
                         disp(sprintf('\n      %s','Other:'));
@@ -187 +199 @@
                         WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer');
                         WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1);
+                        WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','isvzsolution','format','Integer','numel',1,'values',[0 1 2]);
+                        WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','isvelbc','format','Double','numel',1,'scale',1./yts);
 
                         if size(self.loadingforce,2)==3,