Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 17363) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 17364) @@ -15,5 +15,5 @@ /*}}}*/ void LevelsetAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ - int stabilization,finiteelement; + int finiteelement; /*Finite element type*/ @@ -33,4 +33,5 @@ iomodel->FetchDataToInput(elements,VxEnum); iomodel->FetchDataToInput(elements,VyEnum); + iomodel->FetchDataToInput(elements,MasstransportCalvingrateEnum); } @@ -93,8 +94,11 @@ /*Intermediaries */ const int dim = 2; // solve for LSF in horizontal plane only + int i, row, col; IssmDouble kappa; IssmDouble Jdet, dt, D_scalar; IssmDouble h,hx,hy,hz; - IssmDouble vel,vx,vy,bx,by; + IssmDouble vel,v[dim]; + IssmDouble calvingrate, c[dim]; + IssmDouble dlsf[dim], norm_dlsf, normal[dim]; IssmDouble* xyz_list = NULL; @@ -107,5 +111,5 @@ IssmDouble* B = xNew(dim*numnodes); IssmDouble* Bprime = xNew(dim*numnodes); - IssmDouble D[dim][dim], K[dim][dim]; + IssmDouble D[dim][dim]; /*Retrieve all inputs and parameters*/ @@ -114,4 +118,7 @@ Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* lsf_slopex_input = element->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); + Input* lsf_slopey_input = element->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); + Input* calvingrate_input = element->GetInput(MasstransportCalvingrateEnum); _assert_(calvingrate_input); /* Start looping on the number of gaussian points: */ @@ -136,12 +143,23 @@ GetB(B,element,xyz_list,gauss); GetBprime(Bprime,element,xyz_list,gauss); - vx_input->GetInputValue(&vx,gauss); // in 3D case, add mesh velocity - vy_input->GetInputValue(&vy,gauss); - bx=0.; // horizontal mass change velocities (melt/refreeze/calving) FIXME: insert values from model here - by=0.; - D[0][0]=D_scalar*(vx+bx); - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=D_scalar*(vy+by); + vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity + vy_input->GetInputValue(&v[1],gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],gauss); + lsf_slopey_input->GetInputValue(&dlsf[1],gauss); + calvingrate_input->GetInputValue(&calvingrate,gauss); + + norm_dlsf=0.; + for(i=0;i1.e-10) + for(i=0;iElementSizes(&hx,&hy,&hz); - vel=sqrt(vx*vx + vy*vy) + 1e-14; - h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) ); //FIXME: is this correct? + h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); //FIXME: is this correct? kappa=h*vel/2.; //FIXME: insert suitable value for kappa - //GetBprime(Bprime,element,xyz_list,gauss); // recalculation of Bprime needed? - D[0][0]=D_scalar*kappa; - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=D_scalar*kappa; + for(row=0;rowElementSizes(&hx,&hy,&hz); - vel=sqrt(vx*vx + vy*vy )+1.e-14; - h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) ); - K[0][0]=h/(2.*vel)*vx*vx; K[0][1]=h/(2.*vel)*vx*vy; - K[1][0]=h/(2.*vel)*vy*vx; K[1][1]=h/(2.*vel)*vy*vy; - for(int i=0;ivalues[0],1); @@ -203,7 +220,8 @@ /*Intermediaries */ - int i, ig; + const int dim = 2; + int i, ig, k; IssmDouble Jdet,dt; - IssmDouble phi; + IssmDouble lsf; IssmDouble* xyz_list = NULL; @@ -230,6 +248,8 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - levelset_input->GetInputValue(&phi,gauss); - for(i=0;ivalues[i]+=Jdet*gauss->weight*phi*basis[i]; + + /* old function value */ + levelset_input->GetInputValue(&lsf,gauss); + for(i=0;ivalues[i]+=Jdet*gauss->weight*lsf*basis[i]; } @@ -320,3 +340,136 @@ }/*}}}*/ - +void LevelsetAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/ + + /* Intermediaries */ + int i,k; + IssmDouble dmaxp=0.,dmaxm=0,val=0.; + + /*Initialize vector with number of vertices*/ + int numvertices=femmodel->vertices->NumberOfVertices(); + Element* element; + + Vector* vec_dist_zerolevelset = NULL; + GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum); + + /* set NaN on elements intersected by zero levelset */ + for(i=0;ielements->Size();i++){ + element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + if(element->IsZeroLevelset(MaskIceLevelsetEnum)) + for(k=0;kGetNumberOfVertices();k++) + vec_dist_zerolevelset->SetValue(element->vertices[k]->Sid(),NAN,INS_VAL); + } + + /* set distance on elements intersected by zero levelset */ + for(i=0;ielements->Size();i++){ + element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + if(element->IsZeroLevelset(MaskIceLevelsetEnum)){ + SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element); + + /* Get maximum distance to interface along vertices */ + for(k=0;kGetNumberOfVertices();k++){ + vec_dist_zerolevelset->GetValue(&val,element->vertices[k]->Sid()); + if((val>0.) && (val>dmaxp)) + dmaxp=val; + else if((val<0.) && (valGetValue(&val,i); + if(val==1.) //FIXME: improve check + vec_dist_zerolevelset->SetValue(i,3.*dmaxp,INS_VAL); + else if(val==-1.) + vec_dist_zerolevelset->SetValue(i,3.*dmaxm,INS_VAL); + } + + /*Assemble vector and serialize */ + vec_dist_zerolevelset->Assemble(); + IssmDouble* dist_zerolevelset=vec_dist_zerolevelset->ToMPISerial(); + InputUpdateFromVectorx(femmodel,dist_zerolevelset,MaskIceLevelsetEnum,VertexSIdEnum); + + /*Clean up and return*/ + delete vec_dist_zerolevelset; + delete dist_zerolevelset; +}/*}}}*/ +void LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector* vec_signed_dist, Element* element){/*{{{*/ + + if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) + return; + + /* Intermediaries */ + const int dim=3; + int i,d; + int numvertices=element->GetNumberOfVertices(); + IssmDouble s0[dim], s1[dim], v[dim]; + IssmDouble dist,lsf_old; + + IssmDouble* lsf = xNew(numvertices); + IssmDouble* sign_lsf = xNew(numvertices); + IssmDouble* signed_dist = xNew(numvertices); + IssmDouble* xyz_list = NULL; + IssmDouble* xyz_list_zero = NULL; + + /* retrieve inputs and parameters */ + element->GetVerticesCoordinates(&xyz_list); + element->GetInputListOnVertices(lsf,MaskIceLevelsetEnum); + + /* get sign of levelset function */ + for(i=0;i=0.?1.:-1.); + + element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum); + for(d=0;dGetValue(&lsf_old, element->vertices[i]->Sid()); + if(xIsNan(lsf_old) || fabs(signed_dist[i])SetValue(element->vertices[i]->Sid(),signed_dist[i],INS_VAL); + } + + xDelete(lsf); + xDelete(sign_lsf); + xDelete(signed_dist); + xDelete(xyz_list); + xDelete(xyz_list_zero); +}/*}}}*/ +IssmDouble LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/ + // returns distance d of point q to straight going through points s0, s1 + // d=|a x b|/|b| + // with a=q-s0, b=s1-s0 + + /* Intermediaries */ + const int dim=2; + int i; + IssmDouble a[dim], b[dim]; + IssmDouble norm_b; + + for(i=0;i0.); + + return fabs(a[0]*b[1]-a[1]*b[0])/norm_b; +}/*}}}*/ + Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 17363) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 17364) @@ -11,5 +11,5 @@ class LevelsetAnalysis: public Analysis{ - public: +public: /*Model processing*/ int DofsPerNode(int** doflist,int meshtype,int approximation); @@ -27,9 +27,11 @@ ElementVector* CreatePVector(Element* element); void GetSolutionFromInputs(Vector* solution,Element* element); - void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void UpdateConstraints(FemModel* femmodel); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); - + void SetDistanceOnIntersectedElements(FemModel* femmodel); + void SetDistanceToZeroLevelsetElement(Vector* vec_signed_dist, Element* element); + IssmDouble GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1); }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17363) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17364) @@ -463,4 +463,5 @@ name==LevelsetfunctionSlopeYEnum || name==LevelsetfunctionPicardEnum || + name==MasstransportCalvingrateEnum || name==GradientEnum || name==OldGradientEnum || Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 17363) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 17364) @@ -93,4 +93,11 @@ } + if(islevelset){ + /* set distance on elements intersected by zero level set */ + LevelsetAnalysis* lsfanalysis = new LevelsetAnalysis(); + lsfanalysis->SetDistanceOnIntersectedElements(femmodel); + delete lsfanalysis; + } + while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime. @@ -125,9 +132,8 @@ if(islevelset){ if(VerboseSolution()) _printf0_(" computing movement of ice boundaries\n"); - - /* get slope of lsf for computation of normal on ice domain*/ + /* smoothen slope of lsf for computation of normal on ice domain*/ levelsetfunctionslope_core(femmodel); - /* extrapolate along normal */ + /* extrapolate velocities onto domain with no ice */ Analysis* extanalysis = new ExtrapolationAnalysis(); const int nvars=2; @@ -139,4 +145,5 @@ delete extanalysis; + /* solve level set equation */ analysis = new LevelsetAnalysis(); analysis->Core(femmodel); @@ -149,5 +156,4 @@ int outputs[1] = {IceMaskNodeActivationEnum}; femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1); - } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17363) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17364) @@ -217,4 +217,5 @@ MasstransportPenaltyFactorEnum, MasstransportSpcthicknessEnum, + MasstransportCalvingrateEnum, MasstransportStabilizationEnum, MasstransportVertexPairingEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17363) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17364) @@ -225,4 +225,5 @@ case MasstransportPenaltyFactorEnum : return "MasstransportPenaltyFactor"; case MasstransportSpcthicknessEnum : return "MasstransportSpcthickness"; + case MasstransportCalvingrateEnum : return "MasstransportCalvingrate"; case MasstransportStabilizationEnum : return "MasstransportStabilization"; case MasstransportVertexPairingEnum : return "MasstransportVertexPairing"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17363) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17364) @@ -228,4 +228,5 @@ else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum; + else if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum; else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum; else if (strcmp(name,"MasstransportVertexPairing")==0) return MasstransportVertexPairingEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; else if (strcmp(name,"Surface")==0) return SurfaceEnum; - else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum; + if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum; + else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum; else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum; else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"Results")==0) return ResultsEnum; else if (strcmp(name,"GenericParam")==0) return GenericParamEnum; - else if (strcmp(name,"AdolcParam")==0) return AdolcParamEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"BoolInput")==0) return BoolInputEnum; + if (strcmp(name,"AdolcParam")==0) return AdolcParamEnum; + else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; else if (strcmp(name,"Contour")==0) return ContourEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"VzMesh")==0) return VzMeshEnum; else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; - else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; + if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; + else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"Sset")==0) return SsetEnum; else if (strcmp(name,"Verbose")==0) return VerboseEnum; - else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum; + if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; + else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum; else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum; else if (strcmp(name,"XY")==0) return XYEnum; Index: /issm/trunk-jpl/src/m/classes/masstransport.m =================================================================== --- /issm/trunk-jpl/src/m/classes/masstransport.m (revision 17363) +++ /issm/trunk-jpl/src/m/classes/masstransport.m (revision 17364) @@ -7,4 +7,5 @@ properties (SetAccess=public) spcthickness = NaN; + calvingrate = NaN; isfreesurface = 0; min_thickness = 0; @@ -88,4 +89,8 @@ md = checkfield(md,'fieldname','masstransport.spcthickness','forcing',1); + if(ismember(LevelsetAnalysisEnum(), analyses) & md.transient.islevelset) + md = checkfield(md,'fieldname','masstransport.calvingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0); + end + md = checkfield(md,'fieldname','masstransport.isfreesurface','values',[0 1]); md = checkfield(md,'fieldname','masstransport.hydrostatic_adjustment','values',{'Absolute' 'Incremental'}); @@ -98,4 +103,5 @@ disp(sprintf(' Masstransport solution parameters:')); fielddisplay(obj,'spcthickness','thickness constraints (NaN means no constraint) [m]'); + fielddisplay(obj,'calvingrate','calving rate at given location [m/a]'); fielddisplay(obj,'isfreesurface','do we use free surfaces (FS only) are mass conservation'); fielddisplay(obj,'min_thickness','minimum ice thickness allowed [m]'); @@ -110,5 +116,9 @@ end % }}} function marshall(obj,md,fid) % {{{ + + yts=365.*24.*3600.; + WriteData(fid,'object',obj,'fieldname','spcthickness','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'fieldname','calvingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts); WriteData(fid,'object',obj,'fieldname','isfreesurface','format','Boolean'); WriteData(fid,'object',obj,'fieldname','min_thickness','format','Double'); Index: /issm/trunk-jpl/src/m/classes/masstransport.py =================================================================== --- /issm/trunk-jpl/src/m/classes/masstransport.py (revision 17363) +++ /issm/trunk-jpl/src/m/classes/masstransport.py (revision 17364) @@ -15,4 +15,5 @@ def __init__(self): # {{{ self.spcthickness = float('NaN') + self.calvingrate = float('NaN') self.isfreesurface = 0 self.min_thickness = 0 @@ -30,4 +31,5 @@ string=' Masstransport solution parameters:' string="%s\n%s"%(string,fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]')) + string="%s\n%s"%(string,fielddisplay(self,'calvingrate','calving rate at given location [m/a]')) string="%s\n%s"%(string,fielddisplay(self,'isfreesurface','do we use free surfaces (FS only) are mass conservation')) string="%s\n%s"%(string,fielddisplay(self,'min_thickness','minimum ice thickness allowed [m]')) @@ -68,4 +70,6 @@ md = checkfield(md,'fieldname','masstransport.spcthickness','forcing',1) + if LevelsetAnalysisEnum() in analyses and md.transient.islevelset: + md = checkfield(md,'fieldname','masstransport.calvingrate','NaN',1,'size',[md.mesh.numberofvertices],'>=',0) md = checkfield(md,'fieldname','masstransport.isfreesurface','values',[0,1]) md = checkfield(md,'fieldname','masstransport.hydrostatic_adjustment','values',['Absolute','Incremental']) @@ -77,5 +81,9 @@ # }}} def marshall(self,md,fid): # {{{ + + yts=365.*24.*3600. + WriteData(fid,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) + WriteData(fid,'object',self,'fieldname','calvingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts) WriteData(fid,'object',self,'fieldname','isfreesurface','format','Boolean') WriteData(fid,'object',self,'fieldname','min_thickness','format','Double') Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17363) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 17364) @@ -217,4 +217,5 @@ def MasstransportPenaltyFactorEnum(): return StringToEnum("MasstransportPenaltyFactor")[0] def MasstransportSpcthicknessEnum(): return StringToEnum("MasstransportSpcthickness")[0] +def MasstransportCalvingrateEnum(): return StringToEnum("MasstransportCalvingrate")[0] def MasstransportStabilizationEnum(): return StringToEnum("MasstransportStabilization")[0] def MasstransportVertexPairingEnum(): return StringToEnum("MasstransportVertexPairing")[0] Index: /issm/trunk-jpl/src/m/enum/MasstransportCalvingrateEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MasstransportCalvingrateEnum.m (revision 17364) +++ /issm/trunk-jpl/src/m/enum/MasstransportCalvingrateEnum.m (revision 17364) @@ -0,0 +1,11 @@ +function macro=MasstransportCalvingrateEnum() +%MASSTRANSPORTCALVINGRATEENUM - Enum of MasstransportCalvingrate +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=MasstransportCalvingrateEnum() + +macro=StringToEnum('MasstransportCalvingrate'); Index: /issm/trunk-jpl/test/NightlyRun/test336.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test336.py (revision 17363) +++ /issm/trunk-jpl/test/NightlyRun/test336.py (revision 17364) @@ -35,4 +35,5 @@ R=0.2*(xmax-xmin) md.mask.ice_levelset=D-R +md.masstransport.calvingrate=0.*numpy.ones((md.mesh.numberofvertices,1)) md=solve(md,TransientSolutionEnum())