Changeset 21463

Timestamp:

01/03/17 14:21:29 (8 years ago)

Author:

aleahsommers

Message:

NEW: Added under-relaxation capability and added basal flux as an output

Location:

issm/trunk-jpl/src

Files:

• c/analyses/HydrologySommersAnalysis.cpp (modified) (11 diffs)
• c/cores/hydrology_core.cpp (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• m/classes/hydrologysommers.m (modified) (5 diffs)

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

@@ -120 +120 @@
         iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
-        //iomodel->FetchDataToInput(elements,"md.hydrology.eff_pressure",EffectivePressureEnum);
         iomodel->FindConstant(&frictionlaw,"md.friction.law");
+
         /*Friction law variables*/
         switch(frictionlaw){
@@ -147 +147 @@
         parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
         parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
+   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum));
 }/*}}}*/
 
@@ -178 +179 @@
         /*Get conductivity from inputs*/
         IssmDouble conductivity = GetConductivity(element);
-//if(element->Id()==1){
-//      printf("Conductivity at CreateKMatrix: %g \n",conductivity);
-//}
 
         /* Start  looping on the number of gaussian points: */
@@ -245 +243 @@
         /*Get conductivity from inputs*/
         IssmDouble conductivity = GetConductivity(element);
-//if(element->Id()==1){
-//      printf("Conductivity in CreatePVector: %g \n",conductivity);
-//}
 
         /*Build friction element, needed later: */
@@ -306 +301 @@
                 meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
                 _assert_(meltrate>0.);
-                
+        
                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
                  (
                   meltrate*(1/rho_water-1/rho_ice)
-                  +A*pow(fabs(pressure_ice - pressure_water_old),n-1)*(pressure_ice - pressure_water_old)*gap
+                  +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap
                   -beta*sqrt(vx*vx+vy*vy)
                   +ieb
@@ -354 +349 @@
         element->GetInputListOnNodes(&bed[0],BaseEnum);
 
+        /*Get head from previous time-step and under-relaxation coefficient to use in under-relaxation for nonlinear convergence*/
+   IssmDouble* head_old  = xNew<IssmDouble>(numnodes); 
+        element->GetInputListOnNodes(&head_old[0],HydrologyHeadEnum);
+   IssmDouble relaxation; 
+        element->FindParam(&relaxation,HydrologyRelaxationEnum);
+
         /*Use the dof list to index into the solution vector: */
         for(int i=0;i<numnodes;i++){
@@ -367 +368 @@
                         values[i] = bed[i];
                 }
+
+                /*Under-relaxation*/
+           values[i] = head_old[i] - relaxation*(head_old[i]-values[i]);
 
                 /*Calculate effective pressure*/
@@ -385 +389 @@
         IssmDouble conductivity = GetConductivity(element);
 
-        IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/(2.*NU);
+        IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/NU;
         element->AddInput(HydrologyReynoldsEnum,&reynolds,P0Enum);
 
-        /*Free ressources:*/
+        /*Calculate basal flux for output*/
+        IssmDouble q = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1]);
+        element->AddInput(HydrologyBasalFluxEnum,&q,P1Enum);
+
+        /*Free resources:*/
         xDelete<IssmDouble>(values);
         xDelete<IssmDouble>(thickness);
@@ -395 +403 @@
         xDelete<int>(doflist);
         xDelete<IssmDouble>(eff_pressure);
+   xDelete<IssmDouble>(head_old);
 }/*}}}*/
 void           HydrologySommersAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
@@ -468 +477 @@
         /*Get conductivity from inputs*/
         IssmDouble conductivity = GetConductivity(element);
-//if(element->Id()==1){
-//      printf("Conductivity at gap update: %g \n",conductivity);
-//}
+
         /*Build friction element, needed later: */
         Friction* friction=new Friction(element,2);
@@ -537 +544 @@
         /*Divide by connectivity*/
         newgap = newgap/totalweights;
-        IssmDouble mingap = 1e-5;
+        IssmDouble mingap = 1e-6;
         if(newgap<mingap) newgap=mingap;
 

issm/trunk-jpl/src/c/cores/hydrology_core.cpp

@@ -89 +89 @@
                 if(save_results){
                         if(VerboseSolution()) _printf0_("   saving results \n");
-                        int outputs[3] = {HydrologyHeadEnum,HydrologyGapHeightEnum,EffectivePressureEnum};
-                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
+                        int outputs[4] = {HydrologyHeadEnum,HydrologyGapHeightEnum,EffectivePressureEnum,HydrologyBasalFluxEnum};
+                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4);
                         
                         /*unload results*/

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

@@ -174 +174 @@
         HydrologyReynoldsEnum,
         HydrologyNeumannfluxEnum,
+   HydrologyRelaxationEnum,
+        HydrologyBasalFluxEnum,
         InversionControlParametersEnum,
         InversionControlScalingFactorsEnum,

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

@@ -15 +15 @@
                 spchead         = NaN;
                 neumannflux     = NaN;
+                relaxation      = 0;
         end
         methods
@@ -30 +31 @@
                 end % }}}
                 function self = setdefaultparameters(self) % {{{
+              % Set under-relaxation parameter to be 1 (no under-relaxation of nonlinear iteration)     
+                        self.relaxation=1;
                 end % }}}
                 function md = checkconsistency(self,md,solution,analyses) % {{{
@@ -46 +49 @@
                         md = checkfield(md,'fieldname','hydrology.reynolds','>',0,'size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1);
                         md = checkfield(md,'fieldname','hydrology.neumannflux','timeseries',1,'NaN',1,'Inf',1);
-                        md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices 1]);
+                        md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices 1]);        
+         md = checkfield(md,'fieldname','hydrology.relaxation','>=',0); 
                 end % }}}
                 function disp(self) % {{{
@@ -59 +63 @@
                         fielddisplay(self,'neumannflux','water flux applied along the model boundary (m^2/s)');
                         fielddisplay(self,'spchead','water head constraints (NaN means no constraint) (m)');
+                        fielddisplay(self,'relaxation','under-relaxation coefficient for nonlinear iteration');
                 end % }}}
                 function marshall(self,prefix,md,fid) % {{{
@@ -74 +79 @@
                         WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','neumannflux','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
                         WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','spchead','format','DoubleMat','mattype',1);
+         WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','relaxation','format','Double');
                 end % }}}
         end