Index: /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp (revision 21462) +++ /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp (revision 21463) @@ -120,6 +120,6 @@ 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,4 +147,5 @@ parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum)); }/*}}}*/ @@ -178,7 +179,4 @@ /*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,7 +243,4 @@ /*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,9 +301,9 @@ 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;ivalues[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,4 +349,10 @@ 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(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;iAddInput(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(values); xDelete(thickness); @@ -395,4 +403,5 @@ xDelete(doflist); xDelete(eff_pressure); + xDelete(head_old); }/*}}}*/ void HydrologySommersAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ @@ -468,7 +477,5 @@ /*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,5 +544,5 @@ /*Divide by connectivity*/ newgap = newgap/totalweights; - IssmDouble mingap = 1e-5; + IssmDouble mingap = 1e-6; if(newgapRequestedOutputsx(&femmodel->results,&outputs[0],3); + int outputs[4] = {HydrologyHeadEnum,HydrologyGapHeightEnum,EffectivePressureEnum,HydrologyBasalFluxEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4); /*unload results*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21462) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21463) @@ -174,4 +174,6 @@ HydrologyReynoldsEnum, HydrologyNeumannfluxEnum, + HydrologyRelaxationEnum, + HydrologyBasalFluxEnum, InversionControlParametersEnum, InversionControlScalingFactorsEnum, Index: /issm/trunk-jpl/src/m/classes/hydrologysommers.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologysommers.m (revision 21462) +++ /issm/trunk-jpl/src/m/classes/hydrologysommers.m (revision 21463) @@ -15,4 +15,5 @@ spchead = NaN; neumannflux = NaN; + relaxation = 0; end methods @@ -30,4 +31,6 @@ 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,5 +49,6 @@ 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,4 +63,5 @@ 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,4 +79,5 @@ 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