Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26849) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26850) @@ -118,4 +118,9 @@ case CalvingTestEnum: break; + case CalvingParameterizationEnum: + iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum); + iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum); + iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); + break; default: _error_("Calving law "<AddObject(iomodel->CopyConstantObject("md.calving.speedfactor",CalvingTestSpeedfactorEnum)); + break; + case CalvingParameterizationEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.use_param",CalvingUseParamEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.scale_theta",CalvingScaleThetaEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.amp_alpha",CalvingAmpAlphaEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.midp",CalvingMidpointEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.nonlinearlaw",CalvingNonlinearLawEnum)); break; default: @@ -572,4 +585,40 @@ } } + else if(calvinglaw==CalvingParameterizationEnum){ + + /*Intermediaries*/ + IssmDouble thickness,bed,sealevel; + IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum); + + /*Loop over all elements of this partition*/ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + + int numnodes = element->GetNumberOfNodes(); + Gauss* gauss = element->NewGauss(); + Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); + Input* b_input = element->GetInput(BedEnum); _assert_(b_input); + Input* sl_input = element->GetInput(SealevelEnum); _assert_(sl_input); + + /*Potentially constrain nodes of this element*/ + for(int in=0;inGaussNode(element->GetElementType(),in); + Node* node=element->GetNode(in); + if(!node->IsActive()) continue; + + H_input->GetInputValue(&thickness,gauss); + b_input->GetInputValue(&bed,gauss); + sl_input->GetInputValue(&sealevel,gauss); + if(thicknessApplyConstraint(0,+1.); + } + else { + /* no ice, set no spc */ + node->DofInFSet(0); + } + } + delete gauss; + } + } else if(calvinglaw==CalvingHabEnum){ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26849) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26850) @@ -3441,4 +3441,7 @@ case CalvingCrevasseDepthEnum: this->CalvingCrevasseDepth(); + break; + case CalvingParameterizationEnum: + this->CalvingRateParameterization(); break; default: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26849) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26850) @@ -229,4 +229,5 @@ virtual void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0; virtual void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");}; + virtual void CalvingRateParameterization(void){_error_("not implemented yet");}; virtual void CalvingRateVonmises(void){_error_("not implemented yet");}; virtual void CalvingRateTest(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26849) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26850) @@ -837,4 +837,153 @@ delete gauss; } +} +/*}}}*/ +void Tria::CalvingRateParameterization(){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble calvingratex[NUMVERTICES]; + IssmDouble calvingratey[NUMVERTICES]; + IssmDouble calvingrate[NUMVERTICES]; + IssmDouble lambda1,lambda2,ex,ey,vx,vy,vel; + IssmDouble sigma_vm[NUMVERTICES]; + IssmDouble B,sigma_max,sigma_max_floating,sigma_max_grounded,n; + IssmDouble epse_2,groundedice,bed,sealevel; + IssmDouble mrate, rho_ice, rho_water, thickness, paramX, Hab; + int use_parameter=0; + int nonlinear_law=0; + IssmDouble theta, alpha, midp, gamma; + + /* Get node coordinates and dof list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Retrieve all inputs and parameters we will need*/ + Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); + Input* B_input = this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); + Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); + Input* bs_input = this->GetInput(BedEnum); _assert_(bs_input); + Input* H_input = this->GetInput(ThicknessEnum); _assert_(H_input); + Input* smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input); + Input* smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input); + Input* n_input = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* sl_input = this->GetInput(SealevelEnum); _assert_(sl_input); + Input* mrate_input = this->GetInput(CalvingMeltingrateEnum); _assert_(mrate_input); + + /* Ice and sea water density */ + this->FindParam(&rho_ice,MaterialsRhoIceEnum); + this->FindParam(&rho_water,MaterialsRhoSeawaterEnum); + + /* Use which parameter */ + this->FindParam(&use_parameter, CalvingUseParamEnum); + this->FindParam(&theta, CalvingScaleThetaEnum); + this->FindParam(&alpha, CalvingAmpAlphaEnum); + this->FindParam(&midp, CalvingMidpointEnum); + this->FindParam(&nonlinear_law, CalvingNonlinearLawEnum); + + /* Start looping on the number of vertices: */ + GaussTria* gauss=new GaussTria(); + for(int iv=0;ivGaussVertex(iv); + + /*Get velocity components and thickness*/ + B_input->GetInputValue(&B,gauss); + n_input->GetInputValue(&n,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + gr_input->GetInputValue(&groundedice,gauss); + bs_input->GetInputValue(&bed,gauss); + H_input->GetInputValue(&thickness,gauss); + smax_fl_input->GetInputValue(&sigma_max_floating,gauss); + smax_gr_input->GetInputValue(&sigma_max_grounded,gauss); + vel=sqrt(vx*vx+vy*vy)+1.e-14; + sl_input->GetInputValue(&sealevel,gauss); + mrate_input->GetInputValue(&mrate,gauss); + + /*Compute strain rate and viscosity: */ + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + + /*Get Eigen values*/ + Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]); + _assert_(!xIsNan(lambda1)); + _assert_(!xIsNan(lambda2)); + + /*Process Eigen values (only account for extension)*/ + lambda1 = max(lambda1,0.); + lambda2 = max(lambda2,0.); + + /*Calculate sigma_vm*/ + epse_2 = 1./2. *(lambda1*lambda1 + lambda2*lambda2); + sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n)); + + /*Tensile stress threshold*/ + if(groundedice<0) + sigma_max = sigma_max_floating; + else + sigma_max = sigma_max_grounded; + + switch (use_parameter) { + case 0: + /* bed elevation */ + paramX = bed; + break; + case 1: + /* Height above floatation */ + if (bed>sealevel) paramX = 0.0; + else paramX = thickness - (rho_water/rho_ice) * (sealevel-bed); + break; + case 2: + /* Thicknese */ + paramX = thickness; + break; + case 4: + /* bed elevation+linear curve fitting */ + paramX = bed; + break; + case -1: + /* use nothing, just the mrate*/ + break; + default: + _error_("The parameter is not supported yet!"); + } + + /* Compute the hyperbolic tangent function */ + if ((use_parameter>-0.5) & (use_parameter<3)) { + gamma = 0.5*theta*(1.0-tanh((paramX+midp)/alpha))+(1.0-theta); + if (gamma<0.0) gamma =0.0; + } + else if (use_parameter>=3) { + gamma = alpha*paramX + theta; + if (gamma > 1.0) gamma = 1.0; + if (gamma < 0.0) gamma = 0.0; + } + else gamma = 1; + + /*-------------------------------------------*/ + if (nonlinear_law) { + /*This von Mises type has too strong positive feedback with vel included + * calvingrate[iv] = (mrate+sigma_vm[iv]*vel/sigma_max)*gamma; + */ + Hab = thickness - (rho_water/rho_ice) * (sealevel-bed); + if (Hab < 0.) Hab = 0.; + if (bed > sealevel) Hab = 0.; + + calvingrate[iv] = (mrate+Hab/sigma_max)*gamma; + } + else { + calvingrate[iv] = mrate*gamma; + } + calvingratex[iv] = calvingrate[iv]*vx/vel; + calvingratey[iv] = calvingrate[iv]*vy/vel; + } + + /*Add input*/ + this->AddInput(CalvingratexEnum,&calvingratex[0],P1DGEnum); + this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); + this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); + this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum); + + /*Clean up and return*/ + delete gauss; } /*}}}*/ @@ -4256,4 +4405,10 @@ calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); break; + case CalvingParameterizationEnum: + lsf_slopex_input = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); + if(dim==2) lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); + calvingrate_input = this->GetInput(CalvingCalvingrateEnum); _assert_(calvingrate_input); + meltingrate_input = this->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); + break; default: _error_("Calving law "<GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + calvingrate_input->GetInputValue(&calvingrate,&gauss); + meltingrate_input->GetInputValue(&meltingrate,&gauss); + meltingrate = 0.; + + norm_dlsf=0.; + for(i=0;i1.e-10) + for(i=0;iElementOperationx(&Element::CalvingRateTest); break; + case CalvingParameterizationEnum: + femmodel->ElementOperationx(&Element::CalvingRateParameterization); + break; default: _error_("Caving law "<',0,'nan',1,'Inf',1); + md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'nan',1,'Inf',1); + md = checkfield(md,'fieldname','calving.min_thickness','>=',0,'NaN',1,'Inf',1,'numel',1); + md = checkfield(md,'fieldname','calving.use_param','values',[-1, 0, 1, 2, 3, 4]); + md = checkfield(md,'fieldname','calving.scale_theta','>=',0,'NaN',1,'Inf',1,'numel',1); + md = checkfield(md,'fieldname','calving.amp_alpha','<>',0,'NaN',1,'Inf',1,'numel',1); + md = checkfield(md,'fieldname','calving.midp','NaN',1,'Inf',1,'numel',1); + md = checkfield(md,'fieldname','calving.nonlinearlaw','values',[0, 1]); + end % }}} + function disp(self) % {{{ + disp(sprintf(' Calving test parameters:')); + fielddisplay(self,'stress_threshold_groundedice','sigma_max applied to grounded ice only [Pa]'); + fielddisplay(self,'stress_threshold_floatingice','sigma_max applied to floating ice only [Pa]'); + fielddisplay(self,'min_thickness','minimum thickness below which no ice is allowed [m]'); + fielddisplay(self,'use_param','0 - Use bed elevation, 1 - use heigh above floatation, 2 - use ice thickness'); + fielddisplay(self,'scale_theta','larger than 0, larger theta means more reduction for shallower ice'); + fielddisplay(self,'amp_alpha','alpha'); + fielddisplay(self,'midp','mid-point'); + fielddisplay(self,'nonlinearlaw','use a nonlinear law'); + + end % }}} + function marshall(self,prefix,md,fid) % {{{ + yts=md.constants.yts; + WriteData(fid,prefix,'name','md.calving.law','data',9,'format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_groundedice','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_floatingice','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','use_param','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','scale_theta','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','amp_alpha','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','midp','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','nonlinearlaw','format','Integer'); + end % }}} + end +end