Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 22954) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 22955) @@ -509,6 +509,4 @@ if SEALEVELRISE issm_sources += ./cores/sealevelrise_core.cpp\ - ./cores/sealevelrise_core_eustatic.cpp\ - ./cores/sealevelrise_core_noneustatic.cpp\ ./analyses/SealevelriseAnalysis.cpp endif Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 22954) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 22955) @@ -124,4 +124,5 @@ bool isoceancoupling; bool issmb; + int basalforcingsmodel; /*Fetch data needed: */ @@ -132,4 +133,5 @@ iomodel->FindConstant(&isoceancoupling,"md.transient.isoceancoupling"); iomodel->FindConstant(&issmb,"md.transient.issmb"); + iomodel->FindConstant(&basalforcingsmodel,"md.basalforcings.model"); /*Finite element type*/ @@ -156,6 +158,22 @@ iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum); iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum); + if (basalforcingsmodel==SpatialLinearFloatingMeltRateEnum){ + iomodel->FetchDataToInput(elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum); + iomodel->FetchDataToInput(elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum); + iomodel->FetchDataToInput(elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum); + } iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum); iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum); + + /*Initialize cumdeltalthickness input: unfortunately, we don't have femmodel, so we need to iterate on + *elements, otherwise we would have used InputUpdateFromConstantx: */ + counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->InputUpdateFromConstant(0.0,SealevelriseCumDeltathicknessEnum); + counter++; + } + } /*Get what we need for ocean-induced basal melting*/ @@ -706,6 +724,8 @@ basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); IssmDouble* newthickness = xNew(numnodes); - IssmDouble* deltathickness = xNew(numnodes); + IssmDouble* cumdeltathickness = xNew(numnodes); + IssmDouble* deltathickness = xNew(numnodes); IssmDouble* newbase = xNew(numnodes); + IssmDouble* bed = xNew(numnodes); IssmDouble* newsurface = xNew(numnodes); IssmDouble* oldthickness = xNew(numnodes); @@ -725,5 +745,5 @@ } - /*Get previous base, thickness, surfac and current sealevel:*/ + /*Get previous base, thickness, surfac and current sealevel and bed:*/ basalelement->GetInputListOnNodes(&oldbase[0],BaseEnum); basalelement->GetInputListOnNodes(&oldsurface[0],SurfaceEnum); @@ -731,7 +751,12 @@ basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum); basalelement->GetInputListOnNodes(&sealevel[0],SealevelEnum); - - /*What is the delta thickness forcing the sea-level rise core: */ - for(i=0;iGetInputListOnNodes(&bed[0],BedEnum); + basalelement->GetInputListOnNodes(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum); + + /*What is the delta thickness forcing the sea-level rise core: cumulated over time, hence the +=:*/ + for(i=0;i0.){ //this is an ice sheet: just add thickness to base. - newsurface[i] = oldbase[i]+newthickness[i]; //surface = oldbase + newthickness - newbase[i] = oldbase[i]; //same base: do nothing + /*Update! actually, the bed has moved too!:*/ + newsurface[i] = bed[i]+newthickness[i]; //surface = bed + newthickness + newbase[i] = bed[i]; //new base at new bed } else{ //this is an ice shelf: hydrostatic equilibrium*/ @@ -760,10 +786,12 @@ /*Add input to the element: */ element->AddBasalInput(ThicknessEnum,newthickness,P1Enum); - element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum); element->AddBasalInput(SurfaceEnum,newsurface,P1Enum); element->AddBasalInput(BaseEnum,newbase,P1Enum); + element->AddBasalInput(SealevelriseCumDeltathicknessEnum,cumdeltathickness,P1Enum); + element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum); /*Free ressources:*/ xDelete(newthickness); + xDelete(cumdeltathickness); xDelete(newbase); xDelete(newsurface); @@ -774,4 +802,6 @@ xDelete(phi); xDelete(sealevel); + xDelete(bed); + xDelete(doflist); if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 22954) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 22955) @@ -19,4 +19,6 @@ }/*}}}*/ void SealevelriseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + int geodetic=0; /*Update elements: */ @@ -31,9 +33,40 @@ /*Create inputs: */ + iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum); iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); - iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); - iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum); + //those only if we have requested geodetic computations: + iomodel->FetchData(&geodetic,"md.slr.geodetic"); + if (geodetic){ + char* masktype=NULL; + iomodel->FetchData(&masktype,"md.mask.type"); + if (strcmp(masktype,"maskpsl")==0){ + iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); + iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum); + } + xDelete(masktype); + } iomodel->FetchDataToInput(elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum); + iomodel->FetchDataToInput(elements,"md.slr.spcthickness",SealevelriseSpcthicknessEnum); iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0); + iomodel->FetchDataToInput(elements,"md.geometry.bed",BedEnum); + iomodel->FetchDataToInput(elements,"md.slr.Ngia",SealevelNGiaRateEnum); + iomodel->FetchDataToInput(elements,"md.slr.Ugia",SealevelUGiaRateEnum); + + /*Initialize cumdeltalthickness and sealevel rise rate input: unfortunately, we don't have femmodel, so we + * need to iterate on elements, otherwise we would have used InputUpdateFromConstantx: */ + counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->InputUpdateFromConstant(0.0,SealevelriseCumDeltathicknessEnum); + element->InputUpdateFromConstant(0.0,SealevelNEsaRateEnum); + element->InputUpdateFromConstant(0.0,SealevelUEsaRateEnum); + element->InputUpdateFromConstant(0.0,SealevelRSLRateEnum); + element->InputUpdateFromConstant(0.0,SealevelEustaticMaskEnum); + element->InputUpdateFromConstant(0.0,SealevelEustaticOceanMaskEnum); + counter++; + } + } + iomodel->FetchDataToInput(elements,"md.slr.steric_rate",SealevelriseStericRateEnum); @@ -69,5 +102,7 @@ parameters->AddObject(iomodel->CopyConstantObject("md.slr.abstol",SealevelriseAbstolEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.maxiter",SealevelriseMaxiterEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.slr.loop_increment",SealevelriseLoopIncrementEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.slr.horiz",SealevelriseHorizEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum)); @@ -79,4 +114,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum)); iomodel->FetchData(&elastic,"md.slr.elastic"); @@ -249,33 +285,5 @@ }/*}}}*/ void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - - IssmDouble *deltaS = NULL; - IssmDouble *S = NULL; - int* sidlist = NULL; - int numvertices; - - numvertices= element->GetNumberOfVertices(); - sidlist=xNew(numvertices); - - element->GetVerticesSidList(sidlist); - - deltaS = xNew(numvertices); - for(int i=0;i(numvertices); - element->GetInputListOnVertices(S,SealevelEnum,0); - - /*Add deltaS to S:*/ - for (int i=0;iAddInput(SealevelEnum,S,P1Enum); - - /*Free ressources:*/ - xDelete(sidlist); - xDelete(deltaS); - xDelete(S); + _error_("not implemeneted yet!"); }/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22954) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22955) @@ -100,5 +100,4 @@ dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]); } - /*Build strain rate tensor*/ eps[0][0] = dvx[0]; eps[0][1] = .5*(dvx[1]+dvy[0]); eps[0][2] = .5*(dvx[2]+dvz[0]); @@ -1050,4 +1049,27 @@ gauss->GaussVertex(iv); input->GetInputValue(&pvalue[iv],gauss); + } + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +void Element::GetInputListOnVerticesAtTime(IssmDouble* pvalue, int enumtype, IssmDouble time){/*{{{*/ + + /*Recover input*/ + Input* input=this->GetInput(enumtype); + if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); + + /*Fetch number vertices for this element*/ + int numvertices = this->GetNumberOfVertices(); + + /*Checks in debugging mode*/ + _assert_(pvalue); + + /* Start looping on the number of vertices: */ + Gauss*gauss=this->NewGauss(); + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&pvalue[iv],gauss,time); } @@ -1333,4 +1355,39 @@ } + /*Clean up*/ + xDelete(doflist); + xDelete(values); + +} +/*}}}*/ +void Element::GetVectorFromInputs(Vector* vector,int input_enum,int type, IssmDouble time){/*{{{*/ + + /*Fetch number vertices for this element and allocate arrays*/ + int numvertices = this->GetNumberOfVertices(); + int numnodes = this->GetNumberOfNodes(); + int* doflist = NULL; + IssmDouble* values = NULL; + + switch(type){ + case VertexPIdEnum: + doflist = xNew(numvertices); + values = xNew(numvertices); + /*Fill in values*/ + this->GetVertexPidList(doflist); + this->GetInputListOnVerticesAtTime(values,input_enum,time); + vector->SetValues(numvertices,doflist,values,INS_VAL); + break; + case VertexSIdEnum: + doflist = xNew(numvertices); + values = xNew(numvertices); + /*Fill in values*/ + this->GetVerticesSidList(doflist); + this->GetInputListOnVerticesAtTime(values,input_enum,time); + vector->SetValues(numvertices,doflist,values,INS_VAL); + break; + default: + _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); + } + /*Clean up*/ xDelete(doflist); @@ -2039,16 +2096,26 @@ name==MaterialsRheologyEsEnum || name==MaterialsRheologyEsbarEnum || - name==SealevelEnum || - name==SealevelUmotionEnum || - name==SealevelNmotionEnum || - name==SealevelEmotionEnum || - name==SealevelAbsoluteEnum || - name==SealevelEustaticEnum || - name==SealevelriseDeltathicknessEnum || - name==EsaUmotionEnum || - name==EsaNmotionEnum || - name==EsaEmotionEnum || - name==EsaXmotionEnum || - name==EsaYmotionEnum || + name==SealevelEnum || + name==SealevelUEsaEnum || + name==SealevelUEsaRateEnum || + name==SealevelNEsaEnum || + name==SealevelNEsaRateEnum || + name==SealevelUNorthEsaEnum || + name==SealevelUEastEsaEnum || + name==SealevelRSLEustaticEnum || + name==SealevelRSLEustaticRateEnum || + name==SealevelRSLEnum || + name==SealevelRSLRateEnum || + name==SealevelEustaticMaskEnum || + name==SealevelEustaticOceanMaskEnum || + name==SealevelUGiaEnum || + name==SealevelUGiaRateEnum || + name==SealevelNGiaEnum || + name==SealevelNGiaRateEnum || + name==SealevelriseDeltathicknessEnum || + name==SealevelriseCumDeltathicknessEnum || + name==EsaUmotionEnum || + name==EsaNmotionEnum || + name==EsaEmotionEnum || name==EsaStrainratexxEnum || name==EsaStrainratexyEnum || @@ -2125,4 +2192,32 @@ this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); xDelete(base); + xDelete(values); + +}/*}}}*/ +void Element::SpatialLinearFloatingiceMeltingRate(){/*{{{*/ + + int numvertices = this->GetNumberOfVertices(); + IssmDouble* deepwatermelt = xNew(numvertices); + IssmDouble* deepwaterel = xNew(numvertices); + IssmDouble* upperwaterel = xNew(numvertices); + IssmDouble* base = xNew(numvertices); + IssmDouble* values = xNew(numvertices); + + this->GetInputListOnVertices(base,BaseEnum); + this->GetInputListOnVertices(deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum); + this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum); + this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum); + + for(int i=0;iupperwaterel[i]) values[i]=0; + else if (base[i]AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); + xDelete(base); + xDelete(deepwaterel); + xDelete(deepwatermelt); + xDelete(upperwaterel); xDelete(values); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22954) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22955) @@ -90,4 +90,5 @@ void GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype); + void GetInputListOnVerticesAtTime(IssmDouble* pvalue,int enumtype,IssmDouble time); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); void GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug); @@ -104,4 +105,5 @@ void GetSolutionFromInputsOneDof(Vector* solution,int solutionenum); void GetVectorFromInputs(Vector* vector, int name_enum, int type); + void GetVectorFromInputs(Vector* vector, int name_enum, int type,IssmDouble time); void GetVertexPidList(int* pidlist); void GetVerticesConnectivityList(int* connectivitylist); @@ -134,4 +136,5 @@ bool IsWaterInElement(); void LinearFloatingiceMeltingRate(); + void SpatialLinearFloatingiceMeltingRate(); void MantlePlumeGeothermalFlux(); void MarshallElement(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction,int numanalyses); @@ -336,5 +339,5 @@ virtual void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; virtual void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0; - virtual void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0; + virtual void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22954) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22955) @@ -201,5 +201,5 @@ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 22954) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 22955) @@ -182,5 +182,5 @@ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 22954) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 22955) @@ -189,5 +189,5 @@ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22954) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22955) @@ -1354,5 +1354,5 @@ } - else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum){ + else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){ /*Check that not all nodes are grounded or floating*/ if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded @@ -1905,4 +1905,5 @@ if(!IsIceInElement())return 0.; + if(!(this->inputs->Max(MaskIceLevelsetEnum)<0)) return 0; int domaintype; @@ -4676,5 +4677,4 @@ /*}}}*/ void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){/*{{{*/ - /*early return if we are not on an ice cap OR ocean:*/ if(!(this->inputs->Max(MaskIceLevelsetEnum)<0) && !IsWaterInElement()){ @@ -4778,4 +4778,5 @@ IssmDouble late,longe,re; IssmDouble lati,longi,ri; + bool notfullygrounded=false; /*elastic green function:*/ @@ -4785,4 +4786,7 @@ /*ice properties: */ IssmDouble rho_ice,rho_water,rho_earth; + + /*constants:*/ + IssmDouble constant=0; /*Initialize eustatic component: do not skip this step :):*/ @@ -4795,8 +4799,24 @@ /*early return if we are not on an ice cap:*/ - if(!(this->inputs->Max(MaskIceLevelsetEnum)<0)){ + if(!(this->inputs->Max(MaskIceLevelsetEnum)<=0)){ + constant=0; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum)); *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! return; } + + /*early return if we are fully floating: */ + if (this->inputs->Max(MaskGroundediceLevelsetEnum)<=0){ + constant=0; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum)); + *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! + return; + } + + /*If we are here, we are on ice that is fully grounded or half-way to floating: */ + if ((this->inputs->Min(MaskGroundediceLevelsetEnum))<0){ + notfullygrounded=true; //used later on. + } + + /*Inform mask: */ + constant=1; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum)); /*recover material parameters: */ @@ -4859,11 +4879,44 @@ /*}}}*/ - /*Compute area of element:*/ + /*Compute area of element. Scale it by grounded fraction if not fully grounded: */ area=GetAreaSpherical(); - + if(notfullygrounded){ + IssmDouble phi=0; + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias + area*=phi; + } + /*Compute ice thickness change: */ Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); - deltathickness_input->GetInputAverage(&I); + + /*If we are fully grounded, take the average over the element: */ + if(!notfullygrounded)deltathickness_input->GetInputAverage(&I); + else{ + IssmDouble total_weight=0; + bool mainlyfloating = true; + int point1; + IssmDouble fraction1,fraction2; + + /*Recover portion of element that is grounded*/ + this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); + + /* Start looping on the number of gaussian points and average over these gaussian points: */ + total_weight=0; + I=0; + for(int ig=gauss->begin();igend();ig++){ + IssmDouble Ig=0; + gauss->GaussPoint(ig); + deltathickness_input->GetInputValue(&Ig,gauss); + I+=Ig*gauss->weight; + total_weight+=gauss->weight; + } + I=I/total_weight; + delete gauss; + } /*Compute eustatic compoent:*/ @@ -4927,4 +4980,5 @@ IssmDouble minlong=400; IssmDouble maxlong=-20; + IssmDouble constant=0; /*precomputed elastic green functions:*/ @@ -4947,5 +5001,9 @@ /*early return if we are not on the ocean:*/ - if (!IsWaterInElement())return; + if (!IsWaterInElement()){ + constant=0; this->inputs->AddInput(new TriaInput(SealevelEustaticOceanMaskEnum,&constant,P0Enum)); + return; + } + constant=1; this->inputs->AddInput(new TriaInput(SealevelEustaticOceanMaskEnum,&constant,P0Enum)); /*recover computational flags: */ @@ -5058,5 +5116,5 @@ } /*}}}*/ -void Tria::SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/ /*diverse:*/ @@ -5082,4 +5140,9 @@ IssmDouble* N_elastic= NULL; IssmDouble* E_elastic= NULL; + DoubleVecParam* U_parameter = NULL; + DoubleVecParam* H_parameter = NULL; + IssmDouble* U_values=NULL; + IssmDouble* N_values=NULL; + IssmDouble* E_values=NULL; /*optimization:*/ @@ -5092,4 +5155,7 @@ /*early return if we are not on the ocean or on an ice cap:*/ if(!(this->inputs->Max(MaskIceLevelsetEnum)<0) && !IsWaterInElement()) return; + + /*early return if we are fully floating: */ + if (this->inputs->Max(MaskGroundediceLevelsetEnum)<=0)return; /*recover computational flags: */ @@ -5156,8 +5222,10 @@ /*recover elastic Green's functions for displacement:*/ - DoubleVecParam* U_parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter); - DoubleVecParam* H_parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter); + U_parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter); U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M); - H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M); + if(horiz){ + H_parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter); + H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M); + } /*From Sg, recover water sea level rise:*/ @@ -5171,11 +5239,15 @@ /*initialize: */ U_elastic=xNewZeroInit(gsize); - N_elastic=xNewZeroInit(gsize); - E_elastic=xNewZeroInit(gsize); + if(horiz){ + N_elastic=xNewZeroInit(gsize); + E_elastic=xNewZeroInit(gsize); + } int* indices=xNew(gsize); - IssmDouble* U_values=xNewZeroInit(gsize); - IssmDouble* N_values=xNewZeroInit(gsize); - IssmDouble* E_values=xNewZeroInit(gsize); + U_values=xNewZeroInit(gsize); + if(horiz){ + N_values=xNewZeroInit(gsize); + E_values=xNewZeroInit(gsize); + } IssmDouble alpha; IssmDouble delPhi,delLambda; @@ -5202,33 +5274,47 @@ } dx = x_element-x; dy = y_element-y; dz = z_element-z; - N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + if(horiz){ + N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + } /*Elastic component (from Eq 17 in Adhikari et al, GMD 2015): */ int index=reCast(alpha/PI*(M-1)); U_elastic[i] += U_elastic_precomputed[index]; - N_elastic[i] += H_elastic_precomputed[index]*N_azim; - E_elastic[i] += H_elastic_precomputed[index]*E_azim; + if(horiz){ + N_elastic[i] += H_elastic_precomputed[index]*N_azim; + E_elastic[i] += H_elastic_precomputed[index]*E_azim; + } /*Add all components to the pUp solution vectors:*/ if(this->inputs->Max(MaskIceLevelsetEnum)<0){ U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i]; - N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i]; - E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i]; + if(horiz){ + N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i]; + E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i]; + } } else if(IsWaterInElement()) { U_values[i]+=3*rho_water/rho_earth*area/eartharea*S*U_elastic[i]; - N_values[i]+=3*rho_water/rho_earth*area/eartharea*S*N_elastic[i]; - E_values[i]+=3*rho_water/rho_earth*area/eartharea*S*E_elastic[i]; + if(horiz){ + N_values[i]+=3*rho_water/rho_earth*area/eartharea*S*N_elastic[i]; + E_values[i]+=3*rho_water/rho_earth*area/eartharea*S*E_elastic[i]; + } } } pUp->SetValues(gsize,indices,U_values,ADD_VAL); - pNorth->SetValues(gsize,indices,N_values,ADD_VAL); - pEast->SetValues(gsize,indices,E_values,ADD_VAL); + if(horiz){ + pNorth->SetValues(gsize,indices,N_values,ADD_VAL); + pEast->SetValues(gsize,indices,E_values,ADD_VAL); + } /*free ressources:*/ xDelete(indices); - xDelete(U_values); xDelete(N_values); xDelete(E_values); - xDelete(U_elastic); xDelete(N_elastic); xDelete(E_elastic); + xDelete(U_values); + xDelete(U_elastic); + if(horiz){ + xDelete(N_values); xDelete(E_values); + xDelete(N_elastic); xDelete(E_elastic); + } return; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22954) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22955) @@ -162,5 +162,5 @@ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea); - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea); + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22954) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22955) @@ -1454,5 +1454,5 @@ /*Figure out maximum across the cluster: */ ISSM_MPI_Reduce(&maxvy,&node_maxvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); maxvy=node_maxvy; @@ -1478,5 +1478,5 @@ /*Figure out maximum across the cluster: */ ISSM_MPI_Reduce(&maxvz,&node_maxvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); maxvz=node_maxvz; @@ -1502,5 +1502,5 @@ /*Figure out minimum across the cluster: */ ISSM_MPI_Reduce(&minvel,&node_minvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); minvel=node_minvel; @@ -1526,5 +1526,5 @@ /*Figure out minimum across the cluster: */ ISSM_MPI_Reduce(&minvx,&node_minvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); minvx=node_minvx; @@ -1550,5 +1550,5 @@ /*Figure out minimum across the cluster: */ ISSM_MPI_Reduce(&minvy,&node_minvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); minvy=node_minvy; @@ -1574,5 +1574,5 @@ /*Figure out minimum across the cluster: */ ISSM_MPI_Reduce(&minvz,&node_minvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); minvz=node_minvz; @@ -1619,5 +1619,5 @@ omega_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); - /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ + /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; J+=weight*1/2*pow(dp[0]*dp[0]+dp[1]*dp[1],2)*Jdet*gauss->weight; @@ -1955,5 +1955,5 @@ IssmDouble* values = xNewZeroInit(nlines*ncols); IssmDouble* allvalues = xNew(nlines*ncols); - + /*Fill-in matrix*/ for(int j=0;jSize();j++){ @@ -1964,5 +1964,5 @@ ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); xDelete(values); - + if(save_results)results->AddResult(new GenericExternalResult(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time)); xDelete(allvalues); @@ -2012,5 +2012,5 @@ int domaintype; this->parameters->FindParam(&domaintype,DomainTypeEnum); - + /*1: go throug all elements of this partition and figure out how many * segments we have (corresopnding to levelset = 0)*/ @@ -2132,11 +2132,11 @@ case VelEnum: this->ElementResponsex(responses,VelEnum); break; case FrictionCoefficientEnum: NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break; - default: + default: if(response_descriptor_enum>=Outputdefinition1Enum && response_descriptor_enum <=Outputdefinition100Enum){ IssmDouble double_result = OutputDefinitionsResponsex(this,response_descriptor_enum); *responses=double_result; } - else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); - break; + else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); + break; } @@ -2269,5 +2269,5 @@ thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); - /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ + /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; } @@ -2460,7 +2460,7 @@ analysis->UpdateConstraints(this); delete analysis; - + /*Second, constraints might be time dependent: */ - SpcNodesx(nodes,constraints,parameters,analysis_type); + SpcNodesx(nodes,constraints,parameters,analysis_type); /*Now, update degrees of freedoms: */ @@ -2473,5 +2473,5 @@ IssmDouble *surface = NULL; IssmDouble *bed = NULL; - + if(VerboseSolution()) _printf0_(" updating vertices positions\n"); @@ -2512,5 +2512,5 @@ /*AMR*/ -#if !defined(_HAVE_ADOLC_) +#if !defined(_HAVE_ADOLC_) void FemModel::ReMesh(void){/*{{{*/ @@ -2522,5 +2522,5 @@ int newnumberofvertices = -1; int newnumberofelements = -1; - bool* my_elements = NULL; + bool* my_elements = NULL; int* my_vertices = NULL; int elementswidth = this->GetElementsWidth();//just tria elements in this version @@ -2528,5 +2528,5 @@ bool isgroundingline; - /*Branch to specific amr depending on requested method*/ + /*Branch to specific amr depending on requested method*/ this->parameters->FindParam(&amrtype,AmrTypeEnum); switch(amrtype){ @@ -2541,5 +2541,5 @@ default: _error_("not implemented yet"); } - + /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/ this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices); @@ -2572,5 +2572,5 @@ /*As the domain is 2D, it is not necessary to create nodes for this analysis*/ - if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; + if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes); @@ -2602,5 +2602,5 @@ //SetCurrentConfiguration(analysis_type); - this->analysis_counter=i; + this->analysis_counter=i; /*Now, plug analysis_counter and analysis_type inside the parameters: */ this->parameters->SetParam(this->analysis_counter,AnalysisCounterEnum); @@ -2619,5 +2619,5 @@ ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters); - if(i==0){ + if(i==0){ VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices } @@ -2663,5 +2663,5 @@ /*Insert bedrock from mismip+ setup*/ /*This was used to Misomip project/simulations*/ - + if(VerboseSolution())_printf0_(" call Mismip bedrock adjust module\n"); @@ -2680,5 +2680,5 @@ by = 500./(1.+exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+exp((2./4000.)*(y-80000./2.+24000.))); r[i] = max(bx+by,-720.); - } + } /*insert new bedrock*/ element->AddInput(BedEnum,&r[0],P1Enum); @@ -2693,5 +2693,5 @@ if(VerboseSolution())_printf0_(" call adjust base and thickness module\n"); - + int numvertices = this->GetElementsWidth(); IssmDouble rho_water,rho_ice,density,base_float; @@ -2705,5 +2705,5 @@ for(int el=0;elelements->Size();el++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(el)); - + element->GetInputListOnVertices(&s[0],SurfaceEnum); element->GetInputListOnVertices(&r[0],BedEnum); @@ -2717,14 +2717,14 @@ base_float = rho_ice*s[i]/(rho_ice-rho_water); if(r[i]>base_float){ - b[i] = r[i]; - } + b[i] = r[i]; + } else { - b[i] = base_float; - } + b[i] = base_float; + } if(fabs(sl[i])>0) _error_("Sea level value "<AddInput(BaseEnum,&b[0],P1Enum); } - + /*Delete*/ xDelete(phi); @@ -2762,5 +2762,5 @@ Vector* input_interpolations = NULL; IssmDouble* input_interpolations_serial = NULL; - int* pos = NULL; + int* pos = NULL; IssmDouble value = 0; @@ -2782,5 +2782,5 @@ P0input_interp = xNew(numP0inputs); P1input_enums = xNew(numP1inputs); - P1input_interp = xNew(numP1inputs); + P1input_interp = xNew(numP1inputs); } numP0inputs = 0; @@ -2814,6 +2814,6 @@ } } - - /*Get P0 and P1 inputs over the elements*/ + + /*Get P0 and P1 inputs over the elements*/ pos = xNew(elementswidth); vP0inputs= new Vector(numberofelements*numP0inputs); @@ -2821,14 +2821,14 @@ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - + /*Get P0 inputs*/ for(int j=0;j(element->GetInput(P0input_enums[j])); + TriaInput* input=xDynamicCast(element->GetInput(P0input_enums[j])); input->GetInputAverage(&value); pos[0]=element->Sid()*numP0inputs+j; - /*Insert input in the vector*/ + /*Insert input in the vector*/ vP0inputs->SetValues(1,pos,&value,INS_VAL); } - + /*Get P1 inputs*/ for(int j=0;jvertices[1]->Sid()*numP1inputs+j; pos[2]=element->vertices[2]->Sid()*numP1inputs+j; - /*Insert input in the vector*/ - vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL); + /*Insert input in the vector*/ + vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL); } } @@ -2847,14 +2847,14 @@ P0inputs=vP0inputs->ToMPISerial(); P1inputs=vP1inputs->ToMPISerial(); - + /*Assign pointers*/ - *pnumP0inputs = numP0inputs; - *pP0inputs = P0inputs; + *pnumP0inputs = numP0inputs; + *pP0inputs = P0inputs; *pP0input_enums = P0input_enums; *pP0input_interp = P0input_interp; - *pnumP1inputs = numP1inputs; - *pP1inputs = P1inputs; + *pnumP1inputs = numP1inputs; + *pP1inputs = P1inputs; *pP1input_enums = P1input_enums; - *pP1input_interp = P1input_interp; + *pP1input_interp = P1input_interp; /*Cleanup*/ @@ -2867,5 +2867,5 @@ /*}}}*/ void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/ - + int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh int newnumberofelements = newfemmodel_elements->Size(); //just on the new partition @@ -2883,5 +2883,5 @@ int* P1input_enums = NULL; int* P1input_interp = NULL; - IssmDouble* values = NULL; + IssmDouble* values = NULL; IssmDouble* vector = NULL; IssmDouble* x = NULL;//global, entire old mesh @@ -2895,5 +2895,5 @@ IssmDouble* newyc = NULL;//just on the new partition int* newelementslist = NULL;//just on the new partition - int* sidtoindex = NULL;//global vertices sid to partition index + int* sidtoindex = NULL;//global vertices sid to partition index /*Get old P0 and P1 inputs (entire mesh)*/ @@ -2928,15 +2928,15 @@ /*Insert P0 and P1 inputs into the new elements (just on the new partition)*/ - values=xNew(elementswidth); + values=xNew(elementswidth); for(int i=0;iSize();i++){//just on the new partition Element* element=xDynamicCast(newfemmodel_elements->GetObjectByOffset(i)); /*newP0inputs is just on the new partition*/ for(int j=0;jAddInput(new DoubleInput(P0input_enums[j],newP0inputs[i*numP0inputs+j])); break; - case IntInputEnum: + case IntInputEnum: element->AddInput(new IntInput(P0input_enums[j],reCast(newP0inputs[i*numP0inputs+j]))); break; @@ -2956,5 +2956,5 @@ } } - + /*Cleanup*/ xDelete(P0inputs); @@ -2995,5 +2995,5 @@ if(!this->elements || !this->vertices || !this->results || !this->parameters) return; - + parameters->FindParam(&step,StepEnum); parameters->FindParam(&time,TimeEnum); @@ -3013,5 +3013,5 @@ this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshYEnum, y,numberofvertices,1,step,time)); - + /*Cleanup*/ xDelete(x); @@ -3074,5 +3074,5 @@ /*Assemble*/ vmasklevelset->Assemble(); - + /*Serialize and set output*/ (*pmasklevelset)=vmasklevelset->ToMPISerial(); @@ -3088,5 +3088,5 @@ /*newelementslist is in Matlab indexing*/ - + /*Creating connectivity table*/ int* connectivity=NULL; @@ -3099,10 +3099,10 @@ connectivity[vertexid-1]+=1;//Matlab to C indexing } - } + } /*Create vertex and insert in vertices*/ for(int i=0;iid=i+1; newvertex->sid=i; @@ -3115,6 +3115,6 @@ newvertex->connectivity=connectivity[i]; newvertex->clone=false;//itapopo check this - vertices->AddObject(newvertex); - } + vertices->AddObject(newvertex); + } } @@ -3143,5 +3143,5 @@ } else newtria->element_type_list=NULL; - + /*Element hook*/ int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials) @@ -3149,5 +3149,5 @@ /*retrieve vertices ids*/ int* vertex_ids=xNew(elementswidth); - for(int j=0;j(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing + for(int j=0;j(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing /*Setting the hooks*/ newtria->numanalyses =this->nummodels; @@ -3161,6 +3161,6 @@ /*Clean up*/ xDelete(vertex_ids); - elements->AddObject(newtria); - } + elements->AddObject(newtria); + } } @@ -3172,12 +3172,12 @@ for(int i=0;iAddObject(new Matice(i+1,i,MaticeEnum)); - } - } - + materials->AddObject(new Matice(i+1,i,MaticeEnum)); + } + } + /*Add new constant material property to materials, at the end: */ Matpar *newmatpar=static_cast(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy()); newmatpar->SetMid(newnumberofelements+1); - materials->AddObject(newmatpar);//put it at the end of the materials + materials->AddObject(newmatpar);//put it at the end of the materials } /*}}}*/ @@ -3186,8 +3186,8 @@ int lid=0; for(int j=0;jid=nodecounter+j+1; @@ -3195,12 +3195,12 @@ newnode->lid=lid++; newnode->analysis_enum=analysis_enum; - + /*Initialize coord_system: Identity matrix by default*/ for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0; for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0; - + /*indexing:*/ newnode->indexingupdate=true; - + Analysis* analysis=EnumToAnalysis(analysis_enum); int *doftypes=NULL; @@ -3216,5 +3216,5 @@ /*Stressbalance Horiz*/ if(analysis_enum==StressbalanceAnalysisEnum){ - // itapopo this code is rarely used. + // itapopo this code is rarely used. /*Coordinate system provided, convert to coord_system matrix*/ //XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]); @@ -3231,5 +3231,5 @@ if(!femmodel_vertices) _error_("GetMesh: vertices are NULL."); if(!femmodel_elements) _error_("GetMesh: elements are NULL."); - + int numberofvertices = femmodel_vertices->NumberOfVertices(); int numberofelements = femmodel_elements->NumberOfElements(); @@ -3237,5 +3237,5 @@ IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; + IssmDouble* z = NULL; int* elementslist = NULL; int* elem_vertices = NULL; @@ -3246,5 +3246,5 @@ /*Get vertices coordinates*/ VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ; - + /*Get element vertices*/ elem_vertices = xNew(elementswidth); @@ -3261,5 +3261,5 @@ vid3->SetValue(element->sid,elem_vertices[2],INS_VAL); } - + /*Assemble*/ vid1->Assemble(); @@ -3271,5 +3271,5 @@ id2 = vid2->ToMPISerial(); id3 = vid3->ToMPISerial(); - + /*Construct elements list*/ elementslist=xNew(numberofelements*elementswidth); @@ -3280,5 +3280,5 @@ elementslist[elementswidth*i+2] = reCast(id3[i])+1; //InterpMesh wants Matlab indexing } - + /*Assign pointers*/ *px = x; @@ -3301,5 +3301,5 @@ if(!femmodel_vertices) _error_("GetMesh: vertices are NULL."); if(!femmodel_elements) _error_("GetMesh: elements are NULL."); - + int numberofvertices = femmodel_vertices->Size(); //number of vertices of this partition int numbertotalofvertices = femmodel_vertices->NumberOfVertices(); //number total of vertices (entire mesh) @@ -3308,9 +3308,9 @@ IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; + IssmDouble* z = NULL; int* elementslist = NULL; int* sidtoindex = NULL; int* elem_vertices = NULL; - + /*Get vertices coordinates of this partition*/ sidtoindex = xNewZeroInit(numbertotalofvertices);//entire mesh, all vertices @@ -3318,8 +3318,8 @@ y = xNew(numberofvertices);//just this partitio; z = xNew(numberofvertices);//just this partitio; - + /*Go through in this partition (vertices)*/ for(int i=0;iGetObjectByOffset(i); + Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i); /*Attention: no spherical coordinates*/ x[i]=vertex->GetX(); @@ -3334,5 +3334,5 @@ elementslist = xNew(numberofelements*elementswidth); if(numberofelements*elementswidth<0) _error_("numberofelements negative."); - + for(int i=0;i(femmodel_elements->GetObjectByOffset(i)); @@ -3341,6 +3341,6 @@ elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing - } - + } + /*Assign pointers*/ *px = x; @@ -3348,5 +3348,5 @@ *pz = z; *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type. - *psidtoindex = sidtoindex; //it is ncessary to insert inputs + *psidtoindex = sidtoindex; //it is ncessary to insert inputs /*Cleanup*/ @@ -3359,7 +3359,7 @@ /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/ if(analysis_enum!=StressbalanceAnalysisEnum) return; - + int numberofnodes_analysistype= this->nodes->NumberOfNodes(analysis_enum); - int dofpernode = 2; //vx and vy + int dofpernode = 2; //vx and vy int numberofcols = dofpernode*2; //to keep dofs and flags in the vspc vector int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh @@ -3385,5 +3385,5 @@ newy=xNew(newnumberofvertices);//just the new partition for(int i=0;iGetObjectByOffset(i); + Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i); /*Attention: no spherical coordinates*/ newx[i]=vertex->GetX(); @@ -3393,5 +3393,5 @@ /*Get spcvx and spcvy of old mesh*/ for(int i=0;iconstraints->Size();i++){ - + Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i); if(!constraint->InAnalysis(analysis_enum)) _error_("AMR create constraints for "<GetDof(); int node = spcstatic->GetNodeId(); - IssmDouble spcvalue = spcstatic->GetValue(); + IssmDouble spcvalue = spcstatic->GetValue(); int nodeindex = node-1; - + /*vx and vx flag insertion*/ if(dof==0) {//vx vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL); //vx vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag - } + } /*vy and vy flag insertion*/ if(dof==1){//vy @@ -3423,5 +3423,5 @@ spc,numberofvertices,numberofcols, newx,newy,newnumberofvertices,NULL); - + /*Now, insert the interpolated constraints in the data set (constraints)*/ count=0; @@ -3440,5 +3440,5 @@ /*spcvy*/ if(!xIsNan(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){ - newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); + newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; @@ -3497,11 +3497,11 @@ bool *my_elements = NULL; int *my_vertices = NULL; - - _assert_(newnumberofvertices>0); - _assert_(newnumberofelements>0); + + _assert_(newnumberofvertices>0); + _assert_(newnumberofelements>0); epart=xNew(newnumberofelements); npart=xNew(newnumberofvertices); index=xNew(elementswidth*newnumberofelements); - + for (int i=0;i(newnumberofvertices); @@ -3533,12 +3533,12 @@ for(int i=0;i(epart); - xDelete(npart); + xDelete(npart); xDelete(index); } /*}}}*/ void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/ - + int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements int numberofvertices = this->vertices->NumberOfVertices(); IssmDouble weight = 0.; - IssmDouble* tauxx = NULL; - IssmDouble* tauyy = NULL; - IssmDouble* tauxy = NULL; + IssmDouble* tauxx = NULL; + IssmDouble* tauyy = NULL; + IssmDouble* tauxy = NULL; IssmDouble* totalweight = NULL; IssmDouble* deviatoricstressxx = xNew(elementswidth); @@ -3573,8 +3573,8 @@ Vector* vectauxy = new Vector(numberofvertices); Vector* vectotalweight = new Vector(numberofvertices); - + /*Update the Deviatoric Stress tensor over the elements*/ this->DeviatoricStressx(); - + /*Calculate the Smoothed Deviatoric Stress tensor*/ for(int i=0;ielements->Size();i++){ @@ -3621,5 +3621,5 @@ /*Divide for the total weight*/ for(int i=0;i0); + _assert_(totalweight[i]>0); tauxx[i] = tauxx[i]/totalweight[i]; tauyy[i] = tauyy[i]/totalweight[i]; @@ -3646,5 +3646,5 @@ void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/ - /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. + /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/ @@ -3666,5 +3666,5 @@ /*Get smoothed deviatoric stress tensor*/ this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy); - + /*Integrate the error over elements*/ for(int i=0;ielements->Size();i++){ @@ -3674,5 +3674,5 @@ element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum); element->GetVerticesSidList(elem_vertices); - + /*Integrate*/ element->GetVerticesCoordinates(&xyz_list); @@ -3689,7 +3689,7 @@ ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n]; } - error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 - } - /*Set the error in the global vector*/ + error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 + } + /*Set the error in the global vector*/ sid=element->Sid(); error = sqrt(error);//sqrt(e^2) @@ -3705,5 +3705,5 @@ /*Serialize and set output*/ (*pelementerror)=velementerror->ToMPISerial(); - + /*Cleanup*/ xDelete(smoothedtauxx); @@ -3749,5 +3749,5 @@ Tria* triaelement = xDynamicCast(element); weight = triaelement->GetArea();//the tria area is a choice for the weight - + /*dH/dx*/ vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL); @@ -3817,5 +3817,5 @@ /*Get smoothed deviatoric stress tensor*/ this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy); - + /*Integrate the error over elements*/ for(int i=0;ielements->Size();i++){ @@ -3903,5 +3903,5 @@ IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; + IssmDouble* z = NULL; IssmDouble* xyz_list = NULL; IssmDouble x1,y1,x2,y2,x3,y3; @@ -3912,5 +3912,5 @@ //element->GetVerticesSidList(elem_vertices); int sid = element->Sid(); - element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinates(&xyz_list); x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; @@ -3945,9 +3945,9 @@ _error_("level set type not implemented yet!"); } - + /*Outputs*/ IssmDouble* zerolevelset_points = NULL; int npoints = 0; - + /*Intermediaries*/ int elementswidth = this->GetElementsWidth(); @@ -3962,5 +3962,5 @@ int count,sid; IssmDouble xc,yc,x1,y1,x2,y2,x3,y3; - + /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/ for(int i=0;ielements->Size();i++){ @@ -3969,13 +3969,13 @@ element->GetVerticesSidList(elem_vertices); sid= element->Sid(); - element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinates(&xyz_list); x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; xc = NAN; - yc = NAN; + yc = NAN; Tria* tria = xDynamicCast(element); if(tria->IsIceInElement()){/*verify if there is ice in the element*/ - if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || + if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || abs(levelset[0]*levelset[1])Size(); - + /*Figure out max of ns: */ ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); @@ -4162,5 +4162,5 @@ } } - + /*One last time: */ pUp->Assemble(); @@ -4181,9 +4181,9 @@ int ns,nsmax; - + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ ns = elements->Size(); - - /*First, figure out the surface area of Earth: */ + + /*First, figure out the surface area of Earth: */ for(int i=0;i(elements->GetObjectByOffset(i)); @@ -4209,5 +4209,5 @@ } } - + /*One last time: */ pUp->Assemble(); @@ -4226,5 +4226,5 @@ #endif #ifdef _HAVE_SEALEVELRISE_ -void FemModel::SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/ +void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ /*serialized vectors:*/ @@ -4262,16 +4262,16 @@ if(i(elements->GetObjectByOffset(i)); - element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); + element->SealevelriseEustatic(pRSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); eustatic_cpu+=eustatic_cpu_e; } - if(i%100==0)pSgi->Assemble(); + if(i%loop==0)pRSLgi->Assemble(); } if(VerboseConvergence())_printf0_("\n"); - + /*One last time: */ - pSgi->Assemble(); + pRSLgi->Assemble(); /*Sum all eustatic components from all cpus:*/ @@ -4285,16 +4285,16 @@ } /*}}}*/ -void FemModel::SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution){/*{{{*/ +void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ /*serialized vectors:*/ - IssmDouble* Sg_old=NULL; - + IssmDouble* RSLg_old=NULL; + IssmDouble eartharea=0; IssmDouble eartharea_cpu=0; int ns,nsmax; - + /*Serialize vectors from previous iteration:*/ - Sg_old=pSg_old->ToMPISerial(); + RSLg_old=pRSLg_old->ToMPISerial(); /*Go through elements, and add contribution from each element to the deflection vector wg:*/ @@ -4306,5 +4306,5 @@ eartharea_cpu += element->GetAreaSpherical(); } - + ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); @@ -4317,20 +4317,20 @@ for(int i=0;i(elements->GetObjectByOffset(i)); - element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea); - } - if(i%100==0)pSgo->Assemble(); - } - if(verboseconvolution)if(VerboseConvergence())_printf_("\n"); - + element->SealevelriseNonEustatic(pRSLgo,RSLg_old,latitude,longitude,radius,eartharea); + } + if(i%loop==0)pRSLgo->Assemble(); + } + if(verboseconvolution)if(VerboseConvergence())_printf0_("\n"); + /*Free ressources:*/ - xDelete(Sg_old); -} -/*}}}*/ -void FemModel::SealevelriseRotationalFeedback(Vector* pSgo_rot, Vector* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ + xDelete(RSLg_old); +} +/*}}}*/ +void FemModel::SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ /*serialized vectors:*/ - IssmDouble* Sg_old=NULL; + IssmDouble* RSLg_old=NULL; IssmDouble eartharea=0; IssmDouble eartharea_cpu=0; @@ -4341,5 +4341,5 @@ /*Serialize vectors from previous iteration:*/ - Sg_old=pSg_old->ToMPISerial(); + RSLg_old=pRSLg_old->ToMPISerial(); /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ @@ -4355,5 +4355,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea); + element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,eartharea); moi_list_cpu[0] += moi_list[0]; moi_list_cpu[1] += moi_list[1]; @@ -4399,36 +4399,37 @@ (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi))); - pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times + pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times } /*Assemble mesh velocity*/ - pSgo_rot->Assemble(); + pRSLgo_rot->Assemble(); /*Assign output pointers:*/ - *pIxz=moi_list[0]; - *pIyz=moi_list[1]; - *pIzz=moi_list[2]; + if(pIxz)*pIxz=moi_list[0]; + if(pIyz)*pIyz=moi_list[1]; + if(pIzz)*pIzz=moi_list[2]; /*Free ressources:*/ - xDelete(Sg_old); - -} -/*}}}*/ -void FemModel::SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ + xDelete(RSLg_old); + + +} +/*}}}*/ +void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ /*serialized vectors:*/ - IssmDouble* Sg=NULL; - + IssmDouble* RSLg=NULL; + IssmDouble eartharea=0; IssmDouble eartharea_cpu=0; int ns,nsmax; - + /*Serialize vectors from previous iteration:*/ - Sg=pSg->ToMPISerial(); + RSLg=pRSLg->ToMPISerial(); /*Go through elements, and add contribution from each element to the deflection vector wg:*/ ns = elements->Size(); - + /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ for(int i=0;i(elements->GetObjectByOffset(i)); - element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea); - } - if(i%100==0){ + element->SealevelriseGeodetic(pUp,pNorth,pEast,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz); + } + if(i%loop==0){ pUp->Assemble(); - pNorth->Assemble(); - pEast->Assemble(); - } - } - + if (horiz){ + pNorth->Assemble(); + pEast->Assemble(); + } + } + } + /*One last time: */ pUp->Assemble(); - pNorth->Assemble(); - pEast->Assemble(); + if(horiz){ + pNorth->Assemble(); + pEast->Assemble(); + } + if(VerboseConvergence())_printf0_("\n"); /*Free ressources:*/ - xDelete(Sg); - xDelete(latitude); - xDelete(longitude); - xDelete(radius); - xDelete(xx); - xDelete(yy); - xDelete(zz); -} -/*}}}*/ -IssmDouble FemModel::SealevelriseOceanAverage(Vector* Sg) { /*{{{*/ - - IssmDouble* Sg_serial=NULL; + xDelete(RSLg); +} +/*}}}*/ +IssmDouble FemModel::SealevelriseOceanAverage(Vector* RSLg) { /*{{{*/ + + IssmDouble* RSLg_serial=NULL; IssmDouble oceanvalue,oceanvalue_cpu; IssmDouble oceanarea,oceanarea_cpu; /*Serialize vectors from previous iteration:*/ - Sg_serial=Sg->ToMPISerial(); + RSLg_serial=RSLg->ToMPISerial(); /*Initialize:*/ @@ -4488,15 +4489,15 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); oceanarea_cpu += element->OceanArea(); - oceanvalue_cpu += element->OceanAverage(Sg_serial); + oceanvalue_cpu += element->OceanAverage(RSLg_serial); } ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - + ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); /*Free ressources:*/ - xDelete(Sg_serial); - + xDelete(RSLg_serial); + return oceanvalue/oceanarea; } @@ -4514,5 +4515,5 @@ int* eplzigzag_counter = NULL; int eplflip_lock; - + HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis(); @@ -4521,8 +4522,8 @@ mask=new Vector(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); recurence=new Vector(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); - this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); - this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); + this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); + this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum); - + for (int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); @@ -4542,6 +4543,6 @@ /*Assemble and serialize*/ mask->Assemble(); - serial_mask=mask->ToMPISerial(); - + serial_mask=mask->ToMPISerial(); + xDelete(eplzigzag_counter); xDelete(serial_rec); @@ -4585,83 +4586,8 @@ int sum_counter; ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); counter=sum_counter; *pEplcount = counter; if(VerboseSolution()) _printf0_(" Number of active nodes in EPL layer: "<< counter <<"\n"); - - /*Update dof indexings*/ - this->UpdateConstraintsx(); - -} -/*}}}*/ -void FemModel::HydrologyIDSupdateDomainx(IssmDouble* pIDScount){ /*{{{*/ - - bool isthermal; - Vector* mask = NULL; - Vector* active = NULL; - IssmDouble* serial_mask = NULL; - IssmDouble* serial_active = NULL; - - HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis(); - parameters->FindParam(&isthermal,TransientIsthermalEnum); - - /*When solving a thermal model we update the thawed nodes*/ - if(isthermal){ - /*Step 1: update mask, the mask correspond to thawed nodes (that have a meltingrate)*/ - mask=new Vector(this->nodes->NumberOfNodes(HydrologyDCInefficientAnalysisEnum)); - - for (int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - inefanalysis->HydrologyIDSGetMask(mask,element); - } - /*Assemble and serialize*/ - mask->Assemble(); - serial_mask=mask->ToMPISerial(); - delete mask; - } - /*for other cases we just grab the mask from the initialisation value*/ - else{ - GetVectorFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum); - } - /*Update Mask and elementize*/ - InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskThawedNodeEnum,NodeSIdEnum); - xDelete(serial_mask); - inefanalysis->ElementizeIdsMask(this); - - /*get node mask coherent with element mask*/ - active=new Vector(nodes->NumberOfNodes(HydrologyDCInefficientAnalysisEnum)); - for (int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - inefanalysis->HydrologyIdsGetActive(active,element); - } - - /*Assemble and serialize*/ - active->Assemble(); - serial_active=active->ToMPISerial(); - delete active; - - /*Update node activation accordingly*/ - int counter =0; - for (int i=0;iSize();i++){ - Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); - if(node->InAnalysis(HydrologyDCInefficientAnalysisEnum)){ - if(serial_active[node->Sid()]==1.){ - node->Activate(); - if(!node->IsClone()) counter++; - } - else{ - node->Deactivate(); - } - } - } - - xDelete(serial_active); - delete inefanalysis; - int sum_counter; - ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - counter=sum_counter; - *pIDScount = counter; - if(VerboseSolution()) _printf0_(" Number of active nodes in IDS layer: "<< counter <<"\n"); /*Update dof indexings*/ @@ -4706,5 +4632,5 @@ int sum_counter; ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); counter=sum_counter; *pL2count = counter; @@ -4791,5 +4717,5 @@ } /*}}}*/ -#ifdef _HAVE_JAVASCRIPT_ +#ifdef _HAVE_JAVASCRIPT_ FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/ /*configuration: */ @@ -4806,10 +4732,10 @@ /*From command line arguments, retrieve different filenames needed to create the FemModel: */ solution_type=StringToEnumx(solution); - + /*Create femmodel from input files: */ profiler->Start(MPROCESSOR); this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL); profiler->Stop(MPROCESSOR); - + /*Save communicator in the parameters dataset: */ this->parameters->AddObject(new GenericParam(incomm,FemModelCommEnum)); @@ -4826,5 +4752,5 @@ size_t* poutputbuffersize; - + /*Before we delete the profiler, report statistics for this run: */ profiler->Stop(TOTAL); //final tagging @@ -4839,5 +4765,5 @@ ); _printf0_("\n"); - + /*Before we close the output file, recover the buffer and size:*/ outputbufferparam = xDynamicCast*>(this->parameters->FindParamObject(OutputBufferPointerEnum)); @@ -4879,5 +4805,5 @@ /*Open output file once for all and add output file descriptor to parameters*/ - output_fid=open_memstream(&outputbuffer,&outputsize); + output_fid=open_memstream(&outputbuffer,&outputsize); if(output_fid==NULL)_error_("could not initialize output stream"); this->parameters->SetParam(output_fid,OutputFilePointerEnum); @@ -4887,4 +4813,5 @@ }/*}}}*/ #endif + #if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_) @@ -4920,5 +4847,5 @@ /*Initialize hmaxvertices with NAN*/ hmaxvertices_serial=xNew(numberofvertices); - for(int i=0;iamrbamg->groundingline_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskGroundediceLevelsetEnum); @@ -4927,5 +4854,5 @@ if(this->amrbamg->deviatoricerror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum); } - + if(my_rank==0){ this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); @@ -4947,8 +4874,8 @@ newelementslist=xNew(newnumberofelements*this->GetElementsWidth()); } - ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); /*Assign output pointers*/ @@ -4983,5 +4910,5 @@ /*Create bamg data structures for bamg*/ this->amrbamg = new AmrBamg(); - + /*Get amr parameters*/ this->parameters->FindParam(&hmin,AmrHminEnum); @@ -5007,5 +4934,5 @@ /*Re-create original mesh and put it in bamg structure (only cpu 0)*/ - if(my_rank==0){ + if(my_rank==0){ this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements); } @@ -5021,5 +4948,5 @@ if(!hmaxvertices) _error_("hmaxvertices is NULL!\n"); - + /*Intermediaries*/ int numberofvertices = this->vertices->NumberOfVertices(); @@ -5028,5 +4955,5 @@ switch(levelset_type){ - case MaskGroundediceLevelsetEnum: + case MaskGroundediceLevelsetEnum: threshold = this->amrbamg->groundingline_distance; resolution = this->amrbamg->groundingline_resolution; @@ -5042,5 +4969,5 @@ this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type); if(!verticedistance) _error_("verticedistance is NULL!\n"); - + /*Fill hmaxVertices*/ for(int i=0;ivertices->NumberOfVertices(); IssmDouble* maxlength = xNew(numberofelements); - IssmDouble* error_vertices = xNewZeroInit(numberofvertices); + IssmDouble* error_vertices = xNewZeroInit(numberofvertices); IssmDouble* error_elements = NULL; IssmDouble* x = NULL; @@ -5077,5 +5004,5 @@ /*Fill variables*/ switch(errorestimator_type){ - case ThicknessErrorEstimatorEnum: + case ThicknessErrorEstimatorEnum: threshold = this->amrbamg->thicknesserror_threshold; groupthreshold = this->amrbamg->thicknesserror_groupthreshold; @@ -5102,10 +5029,10 @@ case ThicknessErrorEstimatorEnum: this->amrbamg->thicknesserror_maximum = maxerror;break; case DeviatoricStressErrorEstimatorEnum: this->amrbamg->deviatoricerror_maximum = maxerror;break; - } + } } /*Get mesh*/ this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index); - + /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/ for(int i=0;ivertices->NumberOfVertices(); @@ -5177,6 +5104,6 @@ /*Get vertices coordinates*/ VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; - - /*Get points which level set is zero (center of elements with zero level set)*/ + + /*Get points which level set is zero (center of elements with zero level set)*/ this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type); @@ -5187,7 +5114,7 @@ for(int j=0;jamr->groundingline_distance>0) this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskGroundediceLevelsetEnum); if(this->amr->icefront_distance>0) this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum); - if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror); - if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror); - + if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror); + if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror); + if(my_rank==0){ this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror, - &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); + &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); newz=xNewZeroInit(newnumberofvertices); if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz."); @@ -5242,10 +5169,10 @@ newelementslist=xNew(newnumberofelements*this->GetElementsWidth()); } - ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); - - /*Assign the pointers*/ + ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Assign the pointers*/ (*pnewelementslist) = newelementslist; //Matlab indexing (*pnewx) = newx; @@ -5263,5 +5190,5 @@ /*}}}*/ void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/ - + /*Define variables*/ int my_rank = IssmComm::GetRank(); @@ -5272,5 +5199,5 @@ IssmDouble* z = NULL; int* elements = NULL; - + /*Initialize field as NULL for now*/ this->amr = NULL; @@ -5280,10 +5207,10 @@ this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); - /*Create initial mesh (coarse mesh) in neopz data structure*/ + /*Create initial mesh (coarse mesh) in neopz data structure*/ /*Just CPU #0 should keep AMR object*/ /*Initialize refinement pattern*/ this->SetRefPatterns(); this->amr = new AdaptiveMeshRefinement(); - this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) + this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) /*Get amr parameters*/ this->parameters->FindParam(&this->amr->level_max,AmrLevelMaxEnum); @@ -5298,5 +5225,5 @@ this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum); this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum); - if(my_rank==0){ + if(my_rank==0){ this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements); } @@ -5319,5 +5246,5 @@ /*Output*/ IssmDouble* elementdistance; - + /*Intermediaries*/ int numberofelements = this->elements->NumberOfElements(); @@ -5328,6 +5255,6 @@ IssmDouble xc,yc,x1,y1,x2,y2,x3,y3; int numberofpoints; - - /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/ + + /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/ this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type); @@ -5335,5 +5262,5 @@ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); int sid = element->Sid(); - element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinates(&xyz_list); x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; @@ -5345,9 +5272,9 @@ for(int j=0;jSetValue(sid,mindistance,INS_VAL); xDelete(xyz_list); - } + } /*Assemble*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22954) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22955) @@ -135,12 +135,12 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pX, Vector* pY, IssmDouble* xx, IssmDouble* yy); - void EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); + void EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pX, Vector* pY, IssmDouble* xx, IssmDouble* yy); + void EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); - void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution); - void SealevelriseRotationalFeedback(Vector* pSgo_rot, Vector* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); - void SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); + void SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop); + void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop); + void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); + void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); IssmDouble SealevelriseOceanAverage(Vector* Sg); #endif Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 22954) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 22955) @@ -118,4 +118,7 @@ void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");}; void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");}; + void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");}; + void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; + void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; /*}}}*/ /*Numerics: {{{*/ Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 22954) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 22955) @@ -54,6 +54,11 @@ void damage_core(FemModel* femmodel); void sealevelrise_core(FemModel* femmodel); +void geodetic_core(FemModel* femmodel); +void steric_core(FemModel* femmodel); Vector* sealevelrise_core_eustatic(FemModel* femmodel); -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic); +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* RSLg_eustatic); +void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* RSLg); +void sealevelrise_core_viscous(Vector** pU_gia,Vector** pN_gia,FemModel* femmodel,Vector* RSLg); +void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg); IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel); @@ -68,4 +73,6 @@ void TransferForcing(FemModel* femmodel,int forcingenum); void TransferSealevel(FemModel* femmodel,int forcingenum); +void EarthMassTransport(FemModel* femmodel); +void slrconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs); //solution configuration Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 22954) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 22955) @@ -10,23 +10,467 @@ #include "../solutionsequences/solutionsequences.h" +/*cores:*/ void sealevelrise_core(FemModel* femmodel){ /*{{{*/ - Vector *Sg = NULL; - Vector *Sg_absolute = NULL; - Vector *Sg_eustatic = NULL; - Vector *slr_initial = NULL; + + /*Parameters, variables:*/ + bool save_results; + bool isslr=0; + int solution_type; + + /*Retrieve parameters:*/ + femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ + if(solution_type==SealevelriseSolutionEnum)isslr=1; + + /*Should we be here?:*/ + if(!isslr)return; + + /*Verbose: */ + if(VerboseSolution()) _printf0_(" computing sea level rise\n"); + + /*set SLR configuration: */ + femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); + + /*Run geodetic:*/ + geodetic_core(femmodel); + + /*Run steric core for sure:*/ + steric_core(femmodel); + + /*Save results: */ + if(save_results){ + int numoutputs = 0; + char **requested_outputs = NULL; + + if(VerboseSolution()) _printf0_(" saving results\n"); + femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum); + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} + } + + /*requested dependents: */ + if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx(); +} +/*}}}*/ +void geodetic_core(FemModel* femmodel){ /*{{{*/ + + + /*variables:*/ + Vector *RSLg = NULL; + Vector *RSLg_rate = NULL; + Vector *RSLg_eustatic = NULL; + Vector *U_esa = NULL; + Vector *U_esa_rate = NULL; + Vector *N_esa = NULL; + Vector *N_esa_rate = NULL; + Vector *U_north_esa = NULL; + Vector *U_east_esa = NULL; + Vector *N_gia= NULL; + Vector *U_gia= NULL; + Vector *N_gia_rate= NULL; + Vector *U_gia_rate= NULL; + + /*parameters:*/ + bool iscoupler; + int solution_type; + int configuration_type; + int modelid,earthid; + bool istransientmasstransport; + int frequency,count; + int horiz; + int geodetic=0; + IssmDouble dt; + + /*Should we even be here?:*/ + femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return; + + /*Verbose: */ + if(VerboseSolution()) _printf0_(" computing geodetic sea level rise\n"); + + /*retrieve more parameters:*/ + femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); + femmodel->parameters->FindParam(&frequency,SealevelriseGeodeticRunFrequencyEnum); + femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + + if(iscoupler){ + femmodel->parameters->FindParam(&modelid,ModelIdEnum); + femmodel->parameters->FindParam(&earthid,EarthIdEnum); + femmodel->parameters->FindParam(&istransientmasstransport,TransientIsmasstransportEnum); + } + else{ + /* we are here, we are not running in a coupler, so we will indeed compute SLR, + * so make sure we are identified as being the Earth.:*/ + modelid=1; earthid=1; + /* in addition, if we are running solution_type SealevelriseSolutionEnum, make sure we + * run, irresepective of the time settings:*/ + count=frequency; + } + + /*If we are running in coupled mode, the Earth model needs to run its own mass transport (if + * not already done by the mass trasnport module. For ice caps, they rely on the transient mass + * transport module exclusively:*/ + if(iscoupler) if(modelid==earthid) if(!istransientmasstransport) EarthMassTransport(femmodel); + + /*increment counter, or call solution core if count==frequency:*/ + if (countparameters->SetParam(count,SealevelriseRunCountEnum); + return; + } + + /*call sea-level rise sub cores:*/ + if(iscoupler){ + /*transfer cumulated deltathickness forcing from ice caps to earth model: */ + TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum); + + /*we have accumulated thicknesses, dump them in deltathcikness: */ + if(modelid==earthid)InputDuplicatex(femmodel,SealevelriseCumDeltathicknessEnum,SealevelriseDeltathicknessEnum); + } + + /*run cores:*/ + if(modelid==earthid){ + + /*call eustatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */ + RSLg_eustatic=sealevelrise_core_eustatic(femmodel); + + /*call non-eustatic core (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */ + RSLg=sealevelrise_core_noneustatic(femmodel,RSLg_eustatic); + + /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */ + sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg); + + /*compute viscosus (GIA) geodetic signatures:*/ + sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg); + + /*compute sea-level rise (low-order spherical harmonics coefficients) diagnostics:*/ + sealevelrise_diagnostics(femmodel,RSLg); + + /*recover N_esa = U_esa + RSLg:*/ + N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1); + + /*transform these values into rates (as we only run this once each frequency turn:*/ + N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency)); + U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency)); + N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency)); + U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency)); + RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency)); + + /*get some results into elements:{{{*/ + InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,U_gia_rate,SealevelUGiaRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_gia_rate,SealevelNGiaRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,U_gia,SealevelUGiaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_gia,SealevelNGiaEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,RSLg_eustatic,SealevelRSLEustaticEnum,VertexSIdEnum); + + if (horiz){ + InputUpdateFromVectorx(femmodel,U_north_esa,SealevelUNorthEsaEnum,VertexSIdEnum); // north motion + InputUpdateFromVectorx(femmodel,U_east_esa,SealevelUEastEsaEnum,VertexSIdEnum); // east motion + } /*}}}*/ + } + + + if(iscoupler){ + /*transfer sea level back to ice caps:*/ + TransferSealevel(femmodel,SealevelNEsaRateEnum); + TransferSealevel(femmodel,SealevelNGiaRateEnum); + TransferSealevel(femmodel,SealevelUEsaRateEnum); + TransferSealevel(femmodel,SealevelUGiaRateEnum); + + //reset cumdeltathickness to 0: + InputUpdateFromConstantx(femmodel,0.0,SealevelriseCumDeltathicknessEnum); + } + + /*reset counter to 1:*/ + femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter. + + /*free ressources:{{{*/ + delete RSLg; + delete RSLg_rate; + delete RSLg_eustatic; + delete U_esa; + delete U_esa_rate; + delete N_esa; + delete N_esa_rate; + + if(horiz){ + delete U_north_esa; + delete U_east_esa; + } + delete N_gia; + delete U_gia; + delete N_gia_rate; + delete U_gia_rate; + /*}}}*/ + +} +/*}}}*/ +void steric_core(FemModel* femmodel){ /*{{{*/ + + /*variables:*/ + Vector *bedrock = NULL; + Vector *SL = NULL; Vector *steric_rate_g = NULL; - Vector *U_radial = NULL; - Vector *U_north = NULL; - Vector *U_east = NULL; - bool save_results,isslr,iscoupler; - int configuration_type; - int solution_type; - int numoutputs = 0; - char **requested_outputs = NULL; - - /*additional parameters: */ + Vector *U_esa_rate= NULL; + Vector *N_esa_rate= NULL; + Vector *U_gia_rate= NULL; + Vector *N_gia_rate= NULL; + + /*parameters: */ + bool isslr=0; + int solution_type; + IssmDouble dt; + int geodetic=0; + + /*Retrieve parameters:*/ + femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); + + /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ + if(solution_type==SealevelriseSolutionEnum)isslr=1; + + /*Should we be here?:*/ + if(!isslr)return; + + /*Verbose: */ + if(VerboseSolution()) _printf0_(" computing steric sea level rise\n"); + + /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/ + GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum); + GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum); + GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexSIdEnum); + if(geodetic){ + GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&U_gia_rate,femmodel,SealevelUGiaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&N_gia_rate,femmodel,SealevelNGiaRateEnum,VertexSIdEnum); + } + + /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate * dt*/ + if(geodetic){ + SL->AXPY(N_gia_rate,dt); + SL->AXPY(N_esa_rate,dt); + } + SL->AXPY(steric_rate_g,dt); + + /*compute new bedrock position: */ + if(geodetic){ + bedrock->AXPY(U_esa_rate,dt); + bedrock->AXPY(U_gia_rate,dt); + } + + /*update element inputs:*/ + InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum); + + /*Free ressources:*/ + delete bedrock; + delete SL; + delete steric_rate_g; + if(geodetic){ + delete U_esa_rate; + delete U_gia_rate; + delete N_esa_rate; + delete N_gia_rate; + } +} +/*}}}*/ +Vector* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/ + + /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ + + Vector *RSLgi = NULL; + IssmDouble RSLgi_oceanaverage = 0; + + /*parameters: */ + int configuration_type; int gsize; bool spherical=true; + IssmDouble *latitude = NULL; + IssmDouble *longitude = NULL; + IssmDouble *radius = NULL; + int loop; + + /*outputs:*/ + IssmDouble eustatic; + + /*recover parameters:*/ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); + + /*first, recover lat,long and radius vectors from vertices: */ + VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*Initialize:*/ + RSLgi = new Vector(gsize); + + /*call the eustatic main module: */ + femmodel->SealevelriseEustatic(RSLgi,&eustatic, latitude, longitude, radius,loop); //this computes + + /*we need to average RSLgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */ + RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi); + + /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the + * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/ + RSLgi->Shift(-eustatic-RSLgi_oceanaverage); + + /*save eustatic value for results: */ + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic)); + + + /*clean up and return:*/ + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + return RSLgi; +}/*}}}*/ +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* RSLg_eustatic){ /*{{{*/ + + /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. + non eustatic core of the SLR solution */ + + + Vector *RSLg = NULL; + Vector *RSLg_old = NULL; + + Vector *RSLgo = NULL; //ocean convolution of the perturbation to gravity potential. + Vector *RSLgo_rot= NULL; // rotational feedback + IssmDouble RSLgo_oceanaverage = 0; //average of RSLgo over the ocean. + + /*parameters: */ + int count; + bool save_results; + int gsize; + int configuration_type; + bool spherical=true; + bool converged=true; + bool rotation=true; + bool verboseconvolution=true; + int max_nonlinear_iterations; + IssmDouble eps_rel; + IssmDouble eps_abs; + IssmDouble *latitude = NULL; + IssmDouble *longitude = NULL; + IssmDouble *radius = NULL; + IssmDouble eustatic; + int loop; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&max_nonlinear_iterations,SealevelriseMaxiterEnum); + femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum); + femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); + + /*computational flag: */ + femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum); + + /*first, recover lat,long and radius vectors from vertices: */ + VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*Initialize:*/ + RSLg = new Vector(gsize); + RSLg->Assemble(); + RSLg_eustatic->Copy(RSLg); //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic. + + RSLg_old = new Vector(gsize); + RSLg_old->Assemble(); + + count=1; + converged=false; + + /*Start loop: */ + for(;;){ + + //save pointer to old sea level rise + delete RSLg_old; RSLg_old=RSLg; + + /*Initialize solution vector: */ + RSLg = new Vector(gsize); RSLg->Assemble(); + RSLgo = new Vector(gsize); RSLgo->Assemble(); + + /*call the non eustatic module: */ + femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop); + + /*assemble solution vector: */ + RSLgo->Assemble(); + + if(rotation){ + /*call rotational feedback module: */ + RSLgo_rot = new Vector(gsize); RSLgo_rot->Assemble(); + femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,NULL,NULL,NULL,latitude,longitude,radius); + RSLgo_rot->Assemble(); + + RSLgo->AXPY(RSLgo_rot,1); + } + + /*we need to average RSLgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */ + RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo); + + /*RSLg is the sum of the eustatic term, and the ocean terms: */ + RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1); + RSLg->Shift(-RSLgo_oceanaverage); + + /*convergence criterion:*/ + slrconvergence(&converged,RSLg,RSLg_old,eps_rel,eps_abs); + + /*free ressources: */ + delete RSLgo; + + /*Increase count: */ + count++; + if(converged==true){ + break; + } + if(count>=max_nonlinear_iterations){ + _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); + converged=true; + break; + } + + /*some minor verbosing adjustment:*/ + if(count>1)verboseconvolution=false; + + } + if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n"); + + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + delete RSLg_old; + + return RSLg; +} /*}}}*/ +void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg){ /*{{{*/ + + Vector *U_esa = NULL; + Vector *U_north_esa = NULL; + Vector *U_east_esa = NULL; + + /*parameters: */ + int configuration_type; + int gsize; + bool spherical=true; + IssmDouble *latitude = NULL; IssmDouble *longitude = NULL; @@ -35,108 +479,81 @@ IssmDouble *yy = NULL; IssmDouble *zz = NULL; + int loop; + int horiz; + + /*retrieve some parameters:*/ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum); + femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum); + + /*find size of vectors:*/ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*intialize vectors:*/ + U_esa = new Vector(gsize); + if (horiz){ + U_north_esa = new Vector(gsize); + U_east_esa = new Vector(gsize); + } + + /*retrieve geometric information: */ + VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); + + /*call the elastic main modlule:*/ + femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,latitude,longitude,radius,xx,yy,zz,loop,horiz); + + /*Assign output pointers:*/ + *pU_esa=U_esa; + if(horiz){ + *pU_east_esa=U_east_esa; + *pU_north_esa=U_north_esa; + } + + /*Free ressources: */ + delete longitude; + delete latitude; + delete radius; + delete xx; + delete yy; + delete zz; +} +/*}}}*/ +void sealevelrise_core_viscous(Vector** pU_gia, Vector** pN_gia,FemModel* femmodel,Vector* RSLg){ /*{{{*/ + + /*variables:*/ + Vector *U_gia = NULL; + Vector *N_gia = NULL; + + /*parameters:*/ + int frequency; IssmDouble dt; - /*Recover some parameters: */ - femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); - femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); - - /*first, recover lat,long and radius vectors from vertices: */ - VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); - - /*recover x,y,z vectors from vertices: */ - VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); - - /*Figure out size of g-set deflection vector and allocate solution vector: */ - gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); - - /*several cases here, depending on value of iscoupler and isslr: - solution_type == SealevelriseSolutionEnum) we are running sea level rise core (no coupler) - ( !iscoupler & !isslr) we are not interested in being here :) - ( !iscoupler & isslr) we are running in uncoupled mode - ( iscoupler & isslr) we are running in coupled mode, and better be earth - ( iscoupler & !isslr) we are running in coupled mode, and better be an ice cap - */ - - if(solution_type==SealevelriseSolutionEnum){ - isslr=1; - iscoupler=0; - } - - /*early return: */ - if( !iscoupler & !isslr) return; //we are not interested in being here :) - - /*In what follows we assume we are all running slr, either in coupled, or uncoupled mode:*/ - if(VerboseSolution()) _printf0_(" computing sea level rise\n"); - - /*set configuration: */ - if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); - - /*transfer deltathickness forcing from ice caps to earth model: */ - if(iscoupler) TransferForcing(femmodel,SealevelriseDeltathicknessEnum); - - /*call sea-level rise sub cores:*/ - if(isslr){ - Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS. - - /* The following is for reproducing Farrell&Clark1976 Fig. 1, and aimed for the workshop! - * md.slr.sealevel is considered as the "initial sea-level" where 1 m slr is distributed uniformly over the ocean - * similar strategy may work for computin SAL due to "internal mass distribution" associated with ocean dynamics - * if it breaks the code, it will be reverted. And, we will strategize how we want to accomodate */ - GetVectorFromInputsx(&slr_initial,femmodel,SealevelEnum,VertexSIdEnum); - Sg_eustatic->AXPY(slr_initial,1); - - Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //ocean loading tems (2nd and 5th terms on the RHS of Farrel and Clark) - - /*compute other geodetic signatures, such as absolute sea level chagne, components of 3-D crustal motion: */ - /*Initialize:*/ - U_radial = new Vector(gsize); - U_north = new Vector(gsize); - U_east = new Vector(gsize); - Sg_absolute = new Vector(gsize); - - /*call the geodetic main modlule:*/ - femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz); - - /*Now deal with steric ocean expansion by just shifting Sg by a spatial rate pattern : */ - femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexSIdEnum); - Sg->AXPY(steric_rate_g,dt); - - /*compute: absolute sea level change = relative sea level change + vertical motion*/ - Sg->Copy(Sg_absolute); Sg_absolute->AXPY(U_radial,1); - - /*get results into elements:*/ - InputUpdateFromVectorx(femmodel,U_radial,SealevelUmotionEnum,VertexSIdEnum); // radial displacement - InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum); // north motion - InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum); // east motion - InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelAbsoluteEnum,VertexSIdEnum); //absolute sea level - InputUpdateFromVectorx(femmodel,Sg,SealevelEnum,VertexSIdEnum); //relative sea level - - if(save_results){ - if(VerboseSolution()) _printf0_(" saving results\n"); - femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum); - femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); - } - - if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx(); - - /*Free ressources:*/ - delete Sg; - delete Sg_eustatic; - delete U_radial; - delete U_north; - delete U_east; - delete Sg_absolute; - delete steric_rate_g; - if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} - } - - /*transfer sea-level back to ice caps: */ - if(iscoupler)TransferSealevel(femmodel,SealevelEnum); -} + /*retrieve some parameters:*/ + femmodel->parameters->FindParam(&frequency,SealevelriseGeodeticRunFrequencyEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + + /*recover GIA rates:*/ + GetVectorFromInputsx(&U_gia,femmodel,SealevelUGiaRateEnum,VertexSIdEnum); + GetVectorFromInputsx(&N_gia,femmodel,SealevelNGiaRateEnum,VertexSIdEnum); + + /*we just loaded rates, that's not what's being asked, scale by time:*/ + U_gia->Scale(frequency*dt); + N_gia->Scale(frequency*dt); + + /*Assign output pointers:*/ + *pU_gia=U_gia; + *pN_gia=N_gia; + +} /*}}}*/ +void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg){ /*{{{*/ + + /*compute spherical harmonics deg 1 and deg 2 coefficeints:*/ + +} +/*}}}*/ + +/*support routines:*/ void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/ @@ -211,5 +628,5 @@ *First, build the global delta thickness vector in the earth model: */ nv=femmodel->vertices->NumberOfVertices(); - forcingglobal= new Vector(nv); + GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum); /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/ @@ -232,5 +649,5 @@ /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular * ice cap: */ - forcingglobal->SetValues(M,index,forcingfromcap,INS_VAL); + forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL); xDelete(index); } @@ -372,2 +789,90 @@ } /*}}}*/ +void EarthMassTransport(FemModel* femmodel){ /*{{{*/ + + IssmDouble time,dt; + Vector *oldthickness = NULL; + Vector *newthickness = NULL; + Vector *deltathickness = NULL; + Vector *cumdeltathickness = NULL; + int nv; + + if(VerboseSolution()) _printf0_(" computing earth mass transport\n"); + + /*This mass transport module for the Earth is because we might have thickness variations as spcs + * specified in the md.slr class, outside of what we will get from the icecaps. That's why we get t + * the thickness variations from SealevelriseSpcthicknessEnum.*/ + + /*No mass transport module was called, so we are just going to retrieve the geometry thickness + * at this time step, at prior time step, and plug the difference as deltathickness: */ + femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + nv=femmodel->vertices->NumberOfVertices(); + + GetVectorFromInputsx(&newthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum); + GetVectorFromInputsx(&oldthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum,time-dt); + + /*compute deltathickness: */ + deltathickness = new Vector(nv); + newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1); + + /*plug into elements:*/ + InputUpdateFromVectorx(femmodel,deltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum); + + /*add to cumulated delta thickness: */ + GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); + cumdeltathickness->AXPY(deltathickness,1); + InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum); + + /*free ressources:*/ + delete oldthickness; + delete newthickness; + delete deltathickness; + delete cumdeltathickness; + +} /*}}}*/ +void slrconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ + + bool converged=true; + IssmDouble ndS,nS; + Vector *dRSLg = NULL; + + //compute norm(du) and norm(u) if requested + dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0); + ndS=dRSLg->Norm(NORM_TWO); + + if (xIsNan(ndS)) _error_("convergence criterion is NaN!"); + + if(!xIsNan(eps_rel)){ + nS=RSLg_old->Norm(NORM_TWO); + if (xIsNan(nS)) _error_("convergence criterion is NaN!"); + } + + + //clean up + delete dRSLg; + + //print + if(!xIsNan(eps_rel)){ + if((ndS/nS) " << eps_rel*100 << " %\n"); + converged=false; + } + } + if(!xIsNan(eps_abs)){ + if(ndS " << eps_abs << " \n"); + converged=false; + } + } + + /*assign output*/ + *pconverged=converged; + +} /*}}}*/ Index: sm/trunk-jpl/src/c/cores/sealevelrise_core_eustatic.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core_eustatic.cpp (revision 22954) +++ (revision ) @@ -1,59 +1,0 @@ -/*!\file: sealevelrise_core_eustatic.cpp - * \brief: eustatic core of the SLR solution (terms that are constant with respect to sea-level) - */ - -#include "./cores.h" -#include "../toolkits/toolkits.h" -#include "../classes/classes.h" -#include "../shared/shared.h" -#include "../modules/modules.h" -#include "../solutionsequences/solutionsequences.h" - -Vector* sealevelrise_core_eustatic(FemModel* femmodel){ - - Vector *Sgi = NULL; - IssmDouble Sgi_oceanaverage = 0; - - /*parameters: */ - int configuration_type; - int gsize; - bool spherical=true; - IssmDouble *latitude = NULL; - IssmDouble *longitude = NULL; - IssmDouble *radius = NULL; - - /*outputs:*/ - IssmDouble eustatic; - - /*recover parameters:*/ - femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - - /*first, recover lat,long and radius vectors from vertices: */ - VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); - - /*Figure out size of g-set deflection vector and allocate solution vector: */ - gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); - - /*Initialize:*/ - Sgi = new Vector(gsize); - - /*call the eustatic main module: */ - femmodel->SealevelriseEustatic(Sgi,&eustatic, latitude, longitude, radius); //this computes - - /*we need to average Sgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */ - Sgi_oceanaverage=femmodel->SealevelriseOceanAverage(Sgi); - - /*Sg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the - * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/ - Sgi->Shift(-eustatic-Sgi_oceanaverage); - - /*save eustatic value for results: */ - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelEustaticEnum,-eustatic)); - - /*clean up and return:*/ - xDelete(latitude); - xDelete(longitude); - xDelete(radius); - return Sgi; -} - Index: sm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp (revision 22954) +++ (revision ) @@ -1,177 +1,0 @@ -/*!\file: sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. - * \brief: non eustatic core of the SLR solution - */ - -#include "./cores.h" -#include "../toolkits/toolkits.h" -#include "../classes/classes.h" -#include "../shared/shared.h" -#include "../modules/modules.h" -#include "../solutionsequences/solutionsequences.h" - -void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs); - -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* Sg_eustatic){ /*{{{*/ - - Vector *Sg = NULL; - Vector *Sg_old = NULL; - - Vector *Sgo = NULL; //ocean convolution of the perturbation to gravity potential. - Vector *Sgo_rot= NULL; // rotational feedback - IssmDouble Sgo_oceanaverage = 0; //average of Sgo over the ocean. - - /*parameters: */ - int count; - bool save_results; - int gsize; - int configuration_type; - bool spherical=true; - bool converged=true; - bool rotation=true; - bool verboseconvolution=true; - int max_nonlinear_iterations; - IssmDouble eps_rel; - IssmDouble eps_abs; - IssmDouble *latitude = NULL; - IssmDouble *longitude = NULL; - IssmDouble *radius = NULL; - IssmDouble eustatic; - IssmDouble Ixz; - IssmDouble Iyz; - IssmDouble Izz; - - /*Recover some parameters: */ - femmodel->parameters->FindParam(&max_nonlinear_iterations,SealevelriseMaxiterEnum); - femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum); - femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum); - femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - - /*computational flag: */ - femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum); - - /*first, recover lat,long and radius vectors from vertices: */ - VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); - - /*Figure out size of g-set deflection vector and allocate solution vector: */ - gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); - - /*Initialize:*/ - Sg = new Vector(gsize); - Sg->Assemble(); - Sg_eustatic->Copy(Sg); //first initialize Sg with the eustatic component computed in sealevelrise_core_eustatic. - - Sg_old = new Vector(gsize); - Sg_old->Assemble(); - - count=1; - converged=false; - - /*Start loop: */ - for(;;){ - - //save pointer to old sea level rise - delete Sg_old; Sg_old=Sg; - - /*Initialize solution vector: */ - Sg = new Vector(gsize); Sg->Assemble(); - Sgo = new Vector(gsize); Sgo->Assemble(); - - /*call the non eustatic module: */ - femmodel->SealevelriseNonEustatic(Sgo, Sg_old, latitude, longitude, radius,verboseconvolution); - - /*assemble solution vector: */ - Sgo->Assemble(); - - if(rotation){ - /*call rotational feedback module: */ - Sgo_rot = new Vector(gsize); Sgo_rot->Assemble(); - femmodel->SealevelriseRotationalFeedback(Sgo_rot,Sg_old,&Ixz,&Iyz,&Izz,latitude,longitude,radius); - Sgo_rot->Assemble(); - - /*save changes in inertia tensor as results: */ - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz)); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz)); - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,Izz)); - - Sgo->AXPY(Sgo_rot,1); - delete Sgo_rot; - } - - /*we need to average Sgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */ - Sgo_oceanaverage=femmodel->SealevelriseOceanAverage(Sgo); - - /*Sg is the sum of the eustatic term, and the ocean terms: */ - Sg_eustatic->Copy(Sg); Sg->AXPY(Sgo,1); delete Sgo; - Sg->Shift(-Sgo_oceanaverage); - - /*convergence criterion:*/ - slrconvergence(&converged,Sg,Sg_old,eps_rel,eps_abs); - - /*Increase count: */ - count++; - if(converged==true){ - break; - } - if(count>=max_nonlinear_iterations){ - _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); - converged=true; - break; - } - - /*some minor verbosing adjustment:*/ - if(count>1)verboseconvolution=false; - } - if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n"); - - xDelete(latitude); - xDelete(longitude); - xDelete(radius); - delete Sg_old; - - return Sg; -} /*}}}*/ - -void slrconvergence(bool* pconverged, Vector* Sg,Vector* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ - - bool converged=true; - IssmDouble ndS,nS; - Vector *dSg = NULL; - - //compute norm(du) and norm(u) if requested - dSg=Sg_old->Duplicate(); Sg_old->Copy(dSg); dSg->AYPX(Sg,-1.0); - ndS=dSg->Norm(NORM_TWO); - - if(!xIsNan(eps_rel)){ - nS=Sg_old->Norm(NORM_TWO); - } - - if (xIsNan(ndS) || xIsNan(nS)) _error_("convergence criterion is NaN!"); - - //clean up - delete dSg; - - //print - if(!xIsNan(eps_rel)){ - if((ndS/nS) " << eps_rel*100 << " %\n"); - converged=false; - } - } - if(!xIsNan(eps_abs)){ - if(ndS " << eps_abs << " \n"); - converged=false; - } - } - - /*assign output*/ - *pconverged=converged; - -} /*}}}*/ - Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 22954) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 22955) @@ -358,5 +358,5 @@ /*Sea level rise: */ - if(isslr || iscoupler) sealevelrise_core(femmodel); + if(isslr) sealevelrise_core(femmodel); /*unload results*/ Index: /issm/trunk-jpl/src/m/classes/maskpsl.m =================================================================== --- /issm/trunk-jpl/src/m/classes/maskpsl.m (revision 22954) +++ /issm/trunk-jpl/src/m/classes/maskpsl.m (revision 22955) @@ -10,4 +10,5 @@ ocean_levelset = NaN; land_levelset = NaN; + glacier_levelset = NaN; end methods (Static) @@ -68,6 +69,8 @@ fielddisplay(self,'ocean_levelset','is the vertex on the ocean ? yes if = 1, no if = 0'); fielddisplay(self,'land_levelset','is the vertex on the land ? yes if = 1, no if = 0'); + fielddisplay(self,'glacier_levelset','is the vertex on a glacier ? yes if = 1, no if = 0'); end % }}} function marshall(self,prefix,md,fid) % {{{ + WriteData(fid,prefix,'name','md.mask.type','data',class(md.mask),'format','String'); WriteData(fid,prefix,'object',self,'class','mask','fieldname','groundedice_levelset','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'class','mask','fieldname','ice_levelset','format','DoubleMat','mattype',1); Index: /issm/trunk-jpl/src/m/classes/mesh3dsurface.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dsurface.m (revision 22954) +++ /issm/trunk-jpl/src/m/classes/mesh3dsurface.m (revision 22955) @@ -103,5 +103,5 @@ end % }}} function disp(obj) % {{{ - disp(sprintf(' 2D tria Mesh (horizontal):')); + disp(sprintf(' 3D tria Mesh (surface):')); disp(sprintf('\n Elements and vertices:')); Index: /issm/trunk-jpl/src/m/classes/model.js =================================================================== --- /issm/trunk-jpl/src/m/classes/model.js (revision 22954) +++ /issm/trunk-jpl/src/m/classes/model.js (revision 22955) @@ -83,13 +83,12 @@ this.levelset = new levelset(); this.calving = new calving(); + this.gia = new gia(); this.love = new fourierlove(); - this.gia = new giaivins(); - this.esa = new esa(); this.autodiff = new autodiff(); this.inversion = new inversion(); this.qmu = new qmu(); - this.amr = new amr(); - + this.amr = new amr(); + this.radaroverlay = new radaroverlay(); this.results = {}; this.outputdefinition = new outputdefinition(); @@ -764,5 +763,5 @@ this.inversion = 0; this.qmu = 0; - this.amr = 0; + this.amr = 0; this.results = 0; Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 22954) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 22955) @@ -38,15 +38,13 @@ steadystate = 0; transient = 0; - levelset = 0; + levelse = 0; calving = 0; - love = 0; - gia = 0; + gia = 0; esa = 0; - + love = 0; autodiff = 0; inversion = 0; qmu = 0; - amr = 0; - + amr = 0; results = 0; outputdefinition = 0; @@ -128,9 +126,4 @@ %2016 October 11 if isa(md.esa,'double'); md.esa=esa(); end - %2017 Dec 21st (needs to be here) - if isempty(md.settings) - disp('Warning: md.settings had to be reset, make sure to adjust md.settings.output_frequency and other fields'); - md.settings = issmsettings(); - end %2017 February 10th if md.settings.solver_residue_threshold==0, @@ -147,4 +140,9 @@ disp('Warning: calvingdev is now calvingvonmises'); md.calving=calvingvonmises(md.calving); + end + %2017 Dec 21st (needs to be here) + if isempty(md.settings) + disp('Warning: md.settings had to be reset, make sure to adjust md.settings.output_frequency and other fields'); + md.settings = issmsettings(); end @@ -328,4 +326,8 @@ end + %lat long + if numel(md.mesh.lat) ==md.mesh.numberofvertices, md.mesh.lat=project2d(md,md.mesh.lat,1); end + if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end + %outputdefinitions for i=1:length(md.outputdefinition.definitions) @@ -361,5 +363,5 @@ end % }}} - function md2 = extract(md,area) % {{{ + function md2 = extract(md,area,varargin) % {{{ %extract - extract a model according to an Argus contour or flag list % @@ -383,6 +385,9 @@ md1=md; + %recover optoins: + options=pairoptions(varargin{:}); + %some checks - if ((nargin~=2) | (nargout~=1)), + if ((nargin<2) | (nargout~=1)), help extract error('extract error message: bad usage'); @@ -396,10 +401,12 @@ %kick out all elements with 3 dirichlets - spc_elem=find(~flag_elem); - spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); - flag=ones(md1.mesh.numberofvertices,1); - flag(spc_node)=0; - pos=find(sum(flag(md1.mesh.elements),2)==0); - flag_elem(pos)=0; + if getfieldvalue(options,'spccheck',1) + spc_elem=find(~flag_elem); + spc_node=sort(unique(md1.mesh.elements(spc_elem,:))); + flag=ones(md1.mesh.numberofvertices,1); + flag(spc_node)=0; + pos=find(sum(flag(md1.mesh.elements),2)==0); + flag_elem(pos)=0; + end %extracted elements and nodes lists @@ -827,4 +834,5 @@ md.calving=extrude(md.calving,md); md.hydrology = extrude(md.hydrology,md); + md.slr = extrude(md.slr,md); %connectivity @@ -1167,5 +1175,6 @@ md.levelset = levelset(); md.calving = calving(); - md.gia = giaivins(); + md.gia = giaivins(); + md.esa = esa(); md.love = fourierlove(); md.esa = esa(); @@ -1173,5 +1182,5 @@ md.inversion = inversion(); md.qmu = qmu(); - md.amr = amr(); + md.amr = amr(); md.radaroverlay = radaroverlay(); md.results = struct(); @@ -1341,7 +1350,7 @@ disp(sprintf('%19s: %-22s -- %s','levelset' ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)')); disp(sprintf('%19s: %-22s -- %s','calving' ,['[1x1 ' class(self.calving) ']'],'parameters for calving')); - disp(sprintf('%19s: %-22s -- %s','gia' ,['[1x1 ' class(self.gia) ']'],'parameters for gia solution')); - disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution')); + disp(sprintf('%19s: %-22s -- %s','gia' ,['[1x1 ' class(self.gia) ']'],'parameters for gia solution')); disp(sprintf('%19s: %-22s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution')); + disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution')); disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters')); disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods')); @@ -1443,9 +1452,13 @@ end % }}} - function savemodeljs(md,modelname,websiteroot) % {{{ - + function savemodeljs(md,modelname,websiteroot,varargin) % {{{ + %the goal of this routine is to save the model as a javascript array that can be included in any html %file: + options=pairoptions(varargin{:}); + optimization=getfieldvalue(options,'optimize',0); + + %disp: disp(['saving model ''' modelname ''' in file ' websiteroot '/js/' modelname '.js']); @@ -1466,4 +1479,12 @@ end + %some optimization: + if optimization==1, + %optimize for plotting only: + if ~ismember(field,{'geometry','mesh','mask'}), + continue; + end + end + %Check that current field is an object if ~isobject(md.(field)) Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 22954) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 22955) @@ -111,8 +111,7 @@ self.levelset = levelset() self.calving = calving() + self.gia = giaivins() self.love = fourierlove() - self.gia = giaivins() self.esa = esa() - self.autodiff = autodiff() self.inversion = inversion() @@ -128,5 +127,5 @@ def properties(self): # {{{ # ordered list of properties since vars(self) is random - return ['mesh', + return ['mesh', 'mask', 'geometry', Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.m =================================================================== --- /issm/trunk-jpl/src/m/classes/sealevelmodel.m (revision 22954) +++ /issm/trunk-jpl/src/m/classes/sealevelmodel.m (revision 22955) @@ -18,4 +18,6 @@ settings = 0; private = 0; + range = 0; + mergedcaps = 0; %}}} end @@ -43,10 +45,10 @@ for i=1:length(slm.icecaps), if slm.icecaps{i}.transient.iscoupler==0, - error(sprintf('sealevelmodel checkconsistenty error: icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name)); + warning(sprintf('sealevelmodel checkconsistenty error: icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name)); end end if slm.earth.transient.iscoupler==0, - error('sealevelmodel checkconsistenty error: earth model should have the transient coupler option turned on!'); + warning('sealevelmodel checkconsistenty error: earth model should have the transient coupler option turned on!'); end @@ -54,8 +56,22 @@ for i=1:length(slm.icecaps), if slm.icecaps{i}.mesh.numberofvertices ~= length(slm.earth.slr.transitions{i}), - error('sealevelmodel checkconsistenty issue with size of transition vectors!'); + error(['sealevelmodel checkconsistenty issue with size of transition vector for ice cap: ' num2str(i) ' name: ' slm.icecaps{i}.miscellaneous.name]); + end + end + + %check that run_frequency is the same everywhere: + for i=1:length(slm.icecaps), + if slm.icecaps{i}.slr.geodetic_run_frequency~=slm.earth.slr.geodetic_run_frequency, + error(sprintf('sealevelmodel checkconsistenty error: icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name)); end end + %make sure steric_rate is the same everywhere: + for i=1:length(slm.icecaps), + md= slm.icecaps{i}; + if ~isempty(find(md.slr.steric_rate - slm.earth.slr.steric_rate(slm.earth.slr.transitions{i}))), + error(sprintf('steric rate on ice cap %s is not the same as for the earth\n',md.miscellaneous.name)); + end + end end @@ -70,4 +86,5 @@ slm.private = private(); slm.cluster = generic(); + slm.range = {}; end %}}} @@ -78,4 +95,91 @@ disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of cpus...)')); disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields')); + disp(sprintf('%19s: %-22s -- %s','range' ,['[1x1 ' class(self.range) ']'],'ranges')); + end % }}} + function self=mergeresults(self) % {{{ + champs=fields(self.icecaps{1}.results.TransientSolution); + for i=1:length(self.mergedcaps)/2, + md=self.mergedcaps{2*(i-1)+1}; trans=self.mergedcaps{2*(i-1)+2}; + icecaps=self.icecaps(self.range{2*(i-1)+2}); + for j=1:length(self.icecaps{1}.results.TransientSolution), + for k=1:length(champs), + if strcmpi(class(icecaps{1}.results.TransientSolution(j).(champs{k})),'double'), + %vertex or element? + if length(icecaps{1}.results.TransientSolution(j).(champs{k}))==icecaps{1}.mesh.numberofvertices, + md.results.TransientSolution(j).(champs{k})=zeros(md.mesh.numberofvertices,1); + for l=1:length(trans), + resultcap=icecaps{l}.results.TransientSolution(j).(champs{k}); + md.results.TransientSolution(j).(champs{k})(trans{l})=resultcap; + end + else + if strcmpi(champs{k},'IceVolume') | strcmpi(champs{k},'IceVolumeAboveFloatation') , + md.results.TransientSolution(j).(champs{k})=0; + for l=1:length(trans), + resultcap=icecaps{l}.results.TransientSolution(j).(champs{k}); + md.results.TransientSolution(j).(champs{k})= md.results.TransientSolution(j).(champs{k})+resultcap; + end + elseif strcmpi(champs{k},'time'), + md.results.TransientSolution(j).(champs{k})= icecaps{1}.results.TransientSolution(j).(champs{k}); + else + continue; + end + end + else + continue; + end + end + end + self.mergedcaps{2*(i-1)+1}=md; + end + end % }}} + function listcaps(self) % {{{ + for i=1:length(self.icecaps), + disp(sprintf('%i: %s',i,self.icecaps{i}.miscellaneous.name)); + end + end % }}} + function viscousiterations(self) % {{{ + for i=1:length(self.icecaps), + ic=self.icecaps{i}; + mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps; + for j=2:length(ic.results.TransientSolution)-1, + mvi=max(mvi,ic.results.TransientSolution(j).StressbalanceConvergenceNumSteps); + end + disp(sprintf('%i, %s: %i',i,self.icecaps{i}.miscellaneous.name,mvi)); + end + end % }}} + function maxtimestep(self) % {{{ + for i=1:length(self.icecaps), + ic=self.icecaps{i}; + mvi=length(ic.results.TransientSolution); + timei=ic.results.TransientSolution(end).time; + disp(sprintf('%i, %s: %i/%g',i,self.icecaps{i}.miscellaneous.name,mvi,timei)); + end + mvi=length(self.earth.results.TransientSolution); + timei=self.earth.results.TransientSolution(end).time; + disp(sprintf('Earth: %i/%g',mvi,timei)); + end % }}} + function self=homogeneize(self,noearth) % {{{ + if nargin==1, + noearth=0; + end + mintimestep=Inf; + for i=1:length(self.icecaps), + ic=self.icecaps{i}; + mintimestep=min(mintimestep, length(ic.results.TransientSolution)); + end + if ~noearth, + mintimestep=min(mintimestep, length(self.earth.results.TransientSolution)); + end + + for i=1:length(self.icecaps), + ic=self.icecaps{i}; + ic.results.TransientSolution=ic.results.TransientSolution(1:mintimestep); + self.icecaps{i}=ic; + end + ic=self.earth; + if ~noearth, + ic.results.TransientSolution=ic.results.TransientSolution(1:mintimestep); + end + self.earth=ic; end % }}} end Index: /issm/trunk-jpl/src/m/classes/slr.js =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.js (revision 22954) +++ /issm/trunk-jpl/src/m/classes/slr.js (revision 22955) @@ -88,26 +88,26 @@ console.log(sprintf(' Sealevelrise solution parameters:')); - fielddisplay(this,'deltathickness','thickness change: ice height equivalent [m]'); - fielddisplay(this,'sealevel','current sea level (prior to computation) [m]'); - fielddisplay(this,'reltol','sea level rise relative convergence criterion, (NaN: not applied)'); - fielddisplay(this,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)'); - fielddisplay(this,'maxiter','maximum number of nonlinear iterations'); - fielddisplay(this,'love_h','load Love number for radial displacement'); - fielddisplay(this,'love_k','load Love number for gravitational potential perturbation'); - fielddisplay(this,'love_l','load Love number for horizontal displacements'); - fielddisplay(this,'tide_love_h','tidal love number (degree 2)'); - fielddisplay(this,'tide_love_k','tidal love number (degree 2)'); - fielddisplay(this,'fluid_love','secular fluid Love number'); - fielddisplay(this,'equatorial_moi','mean equatorial moment of inertia [kg m^2]'); - fielddisplay(this,'polar_moi','polar moment of inertia [kg m^2]'); - fielddisplay(this,'angular_velocity','mean rotational velocity of earth [per second]'); - fielddisplay(this,'rigid','rigid earth graviational potential perturbation'); - fielddisplay(this,'elastic','elastic earth graviational potential perturbation'); - fielddisplay(this,'rotation','rotational earth potential perturbation'); - fielddisplay(this,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); - fielddisplay(this,'steric_rate','rate of steric ocean expansion [mm/yr]'); - fielddisplay(this,'degacc',"accuracy (default .01 deg) for numerical discretization of the Green's functions"); - fielddisplay(this,'transitions','indices into parts of the mesh that will be icecaps'); - fielddisplay(this,'requested_outputs','additional outputs requested'); + fielddisplay(this,'deltathickness','thickness change (main loading of the slr solution core [m]'); + fielddisplay(this,'sealevel','current sea level (prior to computation) [m]'); + fielddisplay(this,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)'); + fielddisplay(this,'abstol','sea level rise absolute convergence criterion, NaN: not applied'); + fielddisplay(this,'maxiter','maximum number of nonlinear iterations'); + fielddisplay(this,'love_h','load Love number for radial displacement'); + fielddisplay(this,'love_k','load Love number for gravitational potential perturbation'); + fielddisplay(this,'love_l','load Love number for horizontal displacements'); + fielddisplay(this,'tide_love_h','tidal love number (degree 2)'); + fielddisplay(this,'tide_love_k','tidal love number (degree 2)'); + fielddisplay(this,'fluid_love','secular fluid Love number'); + fielddisplay(this,'equatorial_moi','mean equatorial moment of inertia [kg m^2]'); + fielddisplay(this,'polar_moi','polar moment of inertia [kg m^2]'); + fielddisplay(this,'angular_velocity','mean rotational velocity of earth [per second]'); + fielddisplay(this,'rigid','rigid earth graviational potential perturbation'); + fielddisplay(this,'elastic','elastic earth graviational potential perturbation'); + fielddisplay(this,'rotation','rotational earth potential perturbation'); + fielddisplay(this,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); + fielddisplay(this,'degacc',"accuracy (default .01 deg) for numerical discretization of the Green's functions"); + fielddisplay(this,'transitions','indices into parts of the mesh that will be icecaps'); + fielddisplay(this,'requested_outputs','additional outputs requested'); + fielddisplay(this,'steric_rate','rate of steric ocean expansion [mm/yr]'); } //}}} this.marshall=function(md,prefix,fid) { //{{{ Index: /issm/trunk-jpl/src/m/classes/slr.m =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.m (revision 22954) +++ /issm/trunk-jpl/src/m/classes/slr.m (revision 22955) @@ -8,4 +8,5 @@ deltathickness = NaN; sealevel = NaN; + spcthickness = NaN; maxiter = 0; reltol = 0; @@ -25,5 +26,11 @@ ocean_area_scaling = 0; steric_rate = 0; %rate of ocean expansion from steric effects. + geodetic_run_frequency = 1; %how many time steps we skip before we run the geodetic part of the solver during transient + geodetic = 0; %compute geodetic SLR? (in addition to steric?) degacc = 0; + loop_increment = 0; + horiz = 0; + Ngia = 0; + Ugia = 0; requested_outputs = {}; transitions = {}; @@ -46,10 +53,12 @@ %maximum of non-linear iterations. self.maxiter=5; + self.loop_increment=200; %computational flags: + self.geodetic=0; self.rigid=1; self.elastic=1; - self.rotation=0; self.ocean_area_scaling=0; + self.rotation=1; %tidal love numbers: @@ -73,4 +82,6 @@ self.steric_rate=0; + %how many time steps we skip before we run SLR solver during transient + self.geodetic_run_frequency=1; %output default: @@ -79,4 +90,7 @@ %transitions should be a cell array of vectors: self.transitions={}; + + %horizontal displacement? (not by default) + self.horiz=0; end % }}} @@ -86,4 +100,5 @@ md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1); md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1); @@ -98,7 +113,12 @@ md = checkfield(md,'fieldname','slr.abstol','size',[1 1]); md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1); + md = checkfield(md,'fieldname','slr.geodetic_run_frequency','size',[1 1],'>=',1); md = checkfield(md,'fieldname','slr.steric_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10); md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1); + md = checkfield(md,'fieldname','slr.loop_increment','NaN',1,'Inf',1,'>=',1); + md = checkfield(md,'fieldname','slr.horiz','NaN',1,'Inf',1,'values',[0 1]); + md = checkfield(md,'fieldname','slr.Ngia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','slr.Ugia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); %check that love numbers are provided at the same level of accuracy: @@ -112,6 +132,21 @@ [els,vertices]=find(maskpos>0); if length(els), - error('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!'); - end + warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!'); + end + + %check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not, + %a coupler to a planet model is provided. + if self.geodetic, + if md.transient.iscoupler, + %we are good; + else + if strcmpi(class(md.mesh),'mesh3dsurface'), + %we are good + else + error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!'); + end + end + end + end % }}} @@ -124,6 +159,7 @@ fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]'); fielddisplay(self,'sealevel','current sea level (prior to computation) [m]'); - fielddisplay(self,'reltol','sea level rise relative convergence criterion, (NaN: not applied)'); - fielddisplay(self,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)'); + fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]'); + fielddisplay(self,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)'); + fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied'); fielddisplay(self,'maxiter','maximum number of nonlinear iterations'); fielddisplay(self,'love_h','load Love number for radial displacement'); @@ -136,9 +172,14 @@ fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]'); fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]'); + fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); + fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)'); + fielddisplay(self,'Ngia','rate of viscous (GIA) geoid expansion (in mm/yr)'); + fielddisplay(self,'Ugia','rate of viscous (GIA) bedrock uplift (in mm/yr)'); + fielddisplay(self,'loop_increment','vector assembly (in the convolution) framentation'); + fielddisplay(self,'geodetic','compute geodetic SLR? ( in addition to steric?) default 0'); + fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)'); fielddisplay(self,'rigid','rigid earth graviational potential perturbation'); fielddisplay(self,'elastic','elastic earth graviational potential perturbation'); fielddisplay(self,'rotation','earth rotational potential perturbation'); - fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); - fielddisplay(self,'steric_rate','rate of steric ocean expansion [mm/yr]'); fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'); fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); @@ -149,5 +190,7 @@ %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2); WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); + %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1); WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','abstol','format','Double'); @@ -166,7 +209,13 @@ WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','Ngia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','Ugia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray'); + WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer'); %process requested outputs @@ -184,4 +233,5 @@ writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness); writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel); + writejs1Darray(fid,[modelname '.slr.spcthickness'],self.spcthickness); writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter); writejsdouble(fid,[modelname '.slr.reltol'],self.reltol); @@ -200,4 +250,5 @@ writejsdouble(fid,[modelname '.slr.rotation'],self.rotation); writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling); + writejsdouble(fid,[modelname '.slr.geodetic_run_frequency'],self.geodetic_run_frequency); writejs1Darray(fid,[modelname '.slr.steric_rate'],self.steric_rate); writejsdouble(fid,[modelname '.slr.degacc'],self.degacc); @@ -205,4 +256,7 @@ writejscellarray(fid,[modelname '.slr.transitions'],self.transitions); end % }}} + function self = extrude(self,md) % {{{ + self.sealevel=project3d(md,'vector',self.sealevel,'type','node'); + end % }}} end end Index: /issm/trunk-jpl/src/m/classes/slr.py =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.py (revision 22954) +++ /issm/trunk-jpl/src/m/classes/slr.py (revision 22955) @@ -60,5 +60,5 @@ string="%s\n%s"%(string,fielddisplay(self,'elastic','elastic earth graviational potential perturbation')) string="%s\n%s"%(string,fielddisplay(self,'rotation','earth rotational potential perturbation')) - string="%s\n%s"%(string,fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)')) + string="%s\n%s"%(string,fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]')) string="%s\n%s"%(string,fielddisplay(self,'steric_rate','rate of steric ocean expansion [mm/yr]')) string="%s\n%s"%(string,fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions')) @@ -71,6 +71,6 @@ #Convergence criterion: absolute, relative and residual - self.reltol=0.01 # 1 per cent - self.abstol=float('NaN') #default + self.reltol=float('NaN') #default + self.abstol=0.001 #1 mm of sea level rise #maximum of non-linear iterations. @@ -87,5 +87,5 @@ self.tide_love_k=0.3055; #degree 2 - #secular fluid love number: + #secular fluid love number: self.fluid_love=0.942; @@ -93,5 +93,5 @@ self.equatorial_moi=8.0077*10**37; # [kg m^2] self.polar_moi =8.0345*10**37; # [kg m^2] - + #mean rotational velocity of earth self.angular_velocity=7.2921*10**-5; # [s^-1]