Index: /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 26295) +++ /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 26296) @@ -61,5 +61,5 @@ /*some constant parameters: */ parameters->AddObject(iomodel->CopyConstantObject("md.esa.hemisphere",EsaHemisphereEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rigid",SolidearthSettingsRigidEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.selfattraction",SolidearthSettingsSelfAttractionEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum)); Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26295) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26296) @@ -92,6 +92,6 @@ int isexternal=0; - IssmDouble* G_rigid = NULL; - IssmDouble* G_rigid_local = NULL; + IssmDouble* G_gravi = NULL; + IssmDouble* G_gravi_local = NULL; IssmDouble* G_viscoelastic = NULL; IssmDouble* G_viscoelastic_interpolated = NULL; @@ -108,5 +108,5 @@ IssmDouble planetradius=0; IssmDouble planetarea=0; - bool rigid=false; + bool selfattraction=false; bool elastic=false; bool viscous=false; @@ -145,5 +145,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.abstol",SolidearthSettingsAbstolEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.maxiter",SolidearthSettingsMaxiterEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rigid",SolidearthSettingsRigidEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.selfattraction",SolidearthSettingsSelfAttractionEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum)); @@ -153,6 +153,7 @@ parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.grdocean",SolidearthSettingsGrdOceanEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.ocean_area_scaling",SolidearthSettingsOceanAreaScalingEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.sealevelloading",SolidearthSettingsSealevelLoadingEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.isgrd",SolidearthSettingsGRDEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.compute_bp_grd",SolidearthSettingsComputeBpGrdEnum)); @@ -228,7 +229,7 @@ parameters->FindParam(&grdmodel,GrdModelEnum); - if(grdmodel!=IvinsEnum){ + if(grdmodel==ElasticEnum){ /*Deal with elasticity {{{*/ - iomodel->FetchData(&rigid,"md.solidearth.settings.rigid"); + iomodel->FetchData(&selfattraction,"md.solidearth.settings.selfattraction"); iomodel->FetchData(&elastic,"md.solidearth.settings.elastic"); iomodel->FetchData(&viscous,"md.solidearth.settings.viscous"); @@ -236,5 +237,5 @@ iomodel->FetchData(&horiz,"md.solidearth.settings.horiz"); - if(rigid){ + if(selfattraction){ /*compute green functions for a range of angles*/ iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); @@ -271,17 +272,17 @@ // // Providing "t" will cause ensurecontiguous to be called. if(viscous){ - IssmDouble* stacktimes=NULL; + IssmDouble* viscoustimes=NULL; ntimesteps=precomputednt; nt=reCast((final_time-start_time)/timeacc)+1; - parameters->AddObject(new IntParam(StackNumStepsEnum,nt)); - /*Initialize stack times:*/ - stacktimes=xNew(nt); + parameters->AddObject(new IntParam(SealevelchangeViscousNumStepsEnum,nt)); + /*Initialize viscous stack times:*/ + viscoustimes=xNew(nt); for(int t=0;tAddObject(new DoubleVecParam(StackTimesEnum,stacktimes,nt)); - parameters->AddObject(new IntParam(StackIndexEnum,0)); - xDelete(stacktimes); + parameters->AddObject(new DoubleVecParam(SealevelchangeViscousTimesEnum,viscoustimes,nt)); + parameters->AddObject(new IntParam(SealevelchangeViscousIndexEnum,0)); + xDelete(viscoustimes); } else { @@ -300,9 +301,9 @@ #endif } - if(rigid){ + if(selfattraction){ #ifdef _HAVE_AD_ - G_rigid=xNew(M,"t"); + G_gravi=xNew(M,"t"); #else - G_rigid=xNew(M); + G_gravi=xNew(M); #endif } @@ -310,5 +311,5 @@ if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); - if(rigid){ + if(selfattraction){ /*compute combined legendre + love number (elastic green function:*/ @@ -327,17 +328,17 @@ #endif } - if(rigid){ + if(selfattraction){ #ifdef _HAVE_AD_ - G_rigid_local=xNew(m,"t"); + G_gravi_local=xNew(m,"t"); #else - G_rigid_local=xNew(m); + G_gravi_local=xNew(m); #endif } - if(rigid){ + if(selfattraction){ for(int i=lower_row;i(i)*degacc * M_PI / 180.0; - G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); + G_gravi_local[i-lower_row]= .5/sin(alpha/2.0); } } @@ -348,6 +349,6 @@ for(int t=0;tAddObject(new DoubleVecParam(SealevelchangeGRigidEnum,G_rigid,M)); + parameters->AddObject(new DoubleVecParam(SealevelchangeGSelfAttractionEnum,G_gravi,M)); if(viscous || elastic){ parameters->AddObject(new DoubleVecParam(SealevelchangeGViscoElasticEnum,G_viscoelastic_interpolated,M*nt)); @@ -507,6 +508,6 @@ /*free resources: */ - xDelete(G_rigid); - xDelete(G_rigid_local); + xDelete(G_gravi); + xDelete(G_gravi_local); if(elastic){ xDelete(love_h); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26295) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26296) @@ -391,5 +391,5 @@ virtual void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0; - virtual void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom)=0; + virtual void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom)=0; virtual void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0; virtual void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0; @@ -400,5 +400,5 @@ virtual void SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0; virtual void SealevelchangeOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom)=0; - virtual void SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0; + virtual void SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0; virtual void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0; virtual void SealevelchangeUpdateViscousFields()=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26295) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26296) @@ -223,5 +223,5 @@ #ifdef _HAVE_SEALEVELCHANGE_ void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; - void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");}; @@ -232,5 +232,5 @@ void SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeUpdateViscousFields(){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26295) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26296) @@ -178,5 +178,5 @@ #ifdef _HAVE_SEALEVELCHANGE_ void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; - void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");}; @@ -187,5 +187,5 @@ void SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeUpdateViscousFields(){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26295) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26296) @@ -185,5 +185,5 @@ #ifdef _HAVE_SEALEVELCHANGE_ void GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; - void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");}; @@ -194,5 +194,5 @@ void SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeUpdateViscousFields(){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26295) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26296) @@ -1937,5 +1937,5 @@ bool istrapeze1, istrapeze2; - IssmDouble phi1,phi2, f11,f12,f21,f22; + IssmDouble phi1,phi2, d,e,f,g,h1,h2; int point1, point2, i0,i1,i2,j0,j1,j2; IssmDouble weights1[3],weights2[3]; @@ -1951,5 +1951,4 @@ IssmDouble xyz3[3][3]={0}; IssmDouble xyz[8][3]={0}; - IssmDouble f1o; IssmDouble w[8][NUMVERTICES]={0}; IssmDouble areasub=0; @@ -2019,12 +2018,12 @@ //If just one levelset is all negative, just take the partitioning of the other, no interaction between them if (levelset1[0]<=0 && levelset1[1]<=0 && levelset1[2]<=0) { - this->GetFractionGeometry(loadweights,&phi2,&point2,&f21,&f22,&istrapeze2,levelset2); - this->GetBarycenterFromLevelset(platbar,plongbar, phi2, f21, f22, late, longe, point2,istrapeze2,planetradius); + this->GetFractionGeometry(loadweights,&phi2,&point2,&f,&g,&istrapeze2,levelset2); + this->GetBarycenterFromLevelset(platbar,plongbar, phi2, f, g, late, longe, point2,istrapeze2,planetradius); *ploadarea=area*phi2; return; } if (levelset2[0]<=0 && levelset2[1]<=0 && levelset2[2]<=0) { - this->GetFractionGeometry(loadweights,&phi1,&point1,&f11,&f12,&istrapeze1,levelset1); - this->GetBarycenterFromLevelset(platbar,plongbar, phi1, f11, f12, late, longe, point1,istrapeze1,planetradius); + this->GetFractionGeometry(loadweights,&phi1,&point1,&d,&e,&istrapeze1,levelset1); + this->GetBarycenterFromLevelset(platbar,plongbar, phi1, d, e, late, longe, point1,istrapeze1,planetradius); *ploadarea=area*phi1; return; @@ -2032,6 +2031,6 @@ - this->GetFractionGeometry(&weights1[0],&phi1,&point1,&f11,&f12,&istrapeze1,levelset1); - this->GetFractionGeometry(&weights2[0],&phi2,&point2,&f21,&f22,&istrapeze2,levelset2); + this->GetFractionGeometry(&weights1[0],&phi1,&point1,&d,&e,&istrapeze1,levelset1); + this->GetFractionGeometry(&weights2[0],&phi2,&point2,&f,&g,&istrapeze2,levelset2); @@ -2039,5 +2038,5 @@ if (istrapeze1==istrapeze2 && point1==point2 && phi1==phi2){ //the two levelsets are redundant: levelset1 = positivescalar * levelset2 - this->GetBarycenterFromLevelset(platbar,plongbar, phi1, f11, f12, late, longe, point1,istrapeze1,planetradius); + this->GetBarycenterFromLevelset(platbar,plongbar, phi1, d, e, late, longe, point1,istrapeze1,planetradius); *ploadarea=area*phi1; for (int i=0;i point_f12} where the levelsets intersect (only used when f1o>=0 and f1o<=1) - if(point2==i0) f1o= f22*(f11-f21)/(f11*f22-f12*f21); - else if(point2==i1) f1o=f21*(1.0-f22-f11)/(f21*(f12-f11)-f12*f22); - else f1o= (f22*(1.0-f21-f11)+f21*f11)/(f22*(f12-f11) +f21*f11); + //Let our element be triangle ABC with: + i0=point1; //A + i1=(point1+1)%3; //B + i2=(point1+2)%3; //C + + j0=point2; //Can be A, B or C + j1=(point2+1)%3; //anticlockwise point from j0 + j2=(point2+2)%3; //clockwise point from j0 + + /* Below we define the relative fractional lengths of ABC where the zero-level contours of the two level sets intersect with ABC and each other. For example D is the intersection of level set 1 with side AB with fractional length d=[AD]/[AB]. + + levelset1 intersects ABC on D and E: + A------D---B + <--d---> + + A-----E----C + <--e--> + + levelset2 intersects ABC on F and G: + j0---F------j1 + <--f-> + + j0-------G--j2 + <---g----> + + levelset1 and 2 intersect on H (when that intersection exists inside the element) + D----H------E + <-h1-> + + F-----H----G + <--h2-> + */ + + if (point2==i0){ + h1= g*(d-f)/(d*g-e*f); + h2= e/g * h1; + } + else if (point2==i1){ + h1=f*(1.0-g-d)/(f*(e-d)-e*g); + h2= 1.0-e/f * h1; + } + else if (point2==i2){ + h1= (g*(1.0-f-d)+f*d)/(g*(e-d) +f*d); + h2= (d*(f+e-1))/(g*(e-d) +f*d); + } //interpolant weights of each point. Any field F[0,1,2] provided at the original vertices [0,1,2] will be equal on point k to sum_i (F[i] * w[k][i]) - w[0][0]=1; - w[1][1]=1; - w[2][2]=1; - w[3][i0]=1.0-f11; w[3][i1]=f11; - w[4][i0]=1.0-f12; w[4][i2]=f12; - w[5][j0]=1.0-f21; w[5][j1]=f21; - w[6][j0]=1.0-f22; w[6][j2]=f22; - for (int j=0;j<3;j++) w[7][j]=w[3][j]*(1.0-f1o) + w[4][j]*f1o; //we interpolate the intersection point between point_f11 and point_f12 at fraction f1o + w[0][0]=1; //A + w[1][1]=1; //B + w[2][2]=1; //C + w[3][i0]=1.0-d; w[3][i1]=d; //D + w[4][i0]=1.0-e; w[4][i2]=e; //E + w[5][j0]=1.0-f; w[5][j1]=f; //F + w[6][j0]=1.0-g; w[6][j2]=g; //G + for (int j=0;j<3;j++) w[7][j]=w[3][j]*(1.0-h1) + w[4][j]*h1; //H: we interpolate the intersection point H between D and E at fraction h1 for (int k=0;k<8;k++){ @@ -2088,68 +2120,92 @@ //point2 can be either i0,i1 or i2. We start the search with i1 and i2 as they have less computational cost in ifs if(point2==i2){ /*{{{*/ - if (f12>1.0-f21){ /*{{{*/ + if (e>1.0-f){ /*{{{*/ if (!istrapeze1 && !istrapeze2){ - tria1[0]=5; tria1[1]= 7; tria1[2]= 4; + tria1[0]=5; tria1[1]= 7; tria1[2]= 4; // FHE + area1=h2*(e+f-1.0)*g; } else if (!istrapeze1 && istrapeze2){ - tria1[0]=i0; tria1[1]= 3; tria1[2]= 7; - tria2[0]=i0; tria2[1]= 7; tria2[2]= 5; + tria1[0]=i0; tria1[1]= 3; tria1[2]= 5; //ADF + area1=d*(1.0-f); + tria2[0]=3; tria2[1]= 7; tria2[2]= 5; //DHF + area2=d*h1*(e+f-1.0); } else if (istrapeze1 && !istrapeze2){ - tria1[0]=7; tria1[1]= 6; tria1[2]= 4; - tria2[0]=4; tria2[1]= 6; tria2[2]= i2; + tria1[0]=7; tria1[1]= 6; tria1[2]= 4; //HGE + area1=g*(1.0-h2)*(e+f-1.0); + tria2[0]=4; tria2[1]= 6; tria2[2]= i2; //EGC + area2=g*(1.0-e); } else { //istrapeze1 && istrapeze2 - tria1[0]=3; tria1[1]= i1; tria1[2]= 7; - tria2[0]=7; tria2[1]= i1; tria2[2]= 6; + tria1[0]=3; tria1[1]= i1; tria1[2]= 6; //DBG + area1=(1.0-d)*(1.0-g); + tria2[0]=3; tria2[1]= 6; tria2[2]= 7; //DGH + area2=g*((1.0-f)*(1.0-h2)+h2*e)+d*(1.0-e-g); } /*}}}*/ } - else if (f12<=1.0-f21){ /*{{{*/ + else if (e<=1.0-f){ /*{{{*/ if (!istrapeze1 && !istrapeze2){ } else if (!istrapeze1 && istrapeze2){ - tria1[0]=i0; tria1[1]= 3; tria1[2]= 4; + tria1[0]=i0; tria1[1]= 3; tria1[2]= 4; //ADE + area1=d*e; } else if (istrapeze1 && !istrapeze2){ - tria1[0]=5; tria1[1]= 6; tria1[2]= i2; + tria1[0]=5; tria1[1]= 6; tria1[2]= i2; //FGC + area1=f*g; } else { //istrapeze1 && istrapeze2 - tria1[0]=3; tria1[1]= i1; tria1[2]= 4; - tria2[0]=4; tria2[1]= i1; tria2[2]= 5; - tria3[0]=5; tria3[1]= i1; tria3[2]= 6; + tria1[0]=3; tria1[1]= i1; tria1[2]= 5; //DBF + area1=(1.0-d)*(1.0-f); + tria2[0]=4; tria2[1]= 3; tria2[2]= 5; //EDF + area2=d*(1.0-e-f); + tria3[0]=5; tria3[1]= i1; tria3[2]= 6; //FBG + area3=f*(1.0-g); } /*}}}*/ } }/*}}}*/ else if(point2==i1){ /*{{{*/ - if (f11>1.0-f22){ /*{{{*/ + if (d>1.0-g){ /*{{{*/ if (!istrapeze1 && !istrapeze2){ - tria1[0]=6; tria1[1]= 3; tria1[2]= 7; + tria1[0]=6; tria1[1]= 3; tria1[2]= 7; //GDH + area1=(1.0-h2)*(d+g-1.0)*f; } else if (!istrapeze1 && istrapeze2){ - tria1[0]=i0; tria1[1]= 6; tria1[2]= 7; - tria2[0]=i0; tria2[1]= 7; tria2[2]= 4; + tria1[0]=i0; tria1[1]= 6; tria1[2]= 4; //AGE + area1=(1.0-g)*e; + tria2[0]=6; tria2[1]= 7; tria2[2]= 4; //GHE + area2=e*(1.0-h1)*(g+d-1.0); } else if (istrapeze1 && !istrapeze2){ - tria1[0]=3; tria1[1]= i1; tria1[2]= 7; - tria2[0]=7; tria2[1]= i1; tria2[2]= 5; + tria1[0]=i1; tria1[1]= 5; tria1[2]= 3; //BFD + area1=(1.0-d)*f; + tria2[0]=3; tria2[1]= 5; tria2[2]= 7; //DFH + area2=f*h2*(d+g-1.0); } else { //istrapeze1 && istrapeze2 - tria1[0]=7; tria1[1]= 5; tria1[2]= 4; - tria2[0]=4; tria2[1]= 5; tria2[2]= i2; + tria1[0]=7; tria1[1]= 5; tria1[2]= 4; //HFE + area1=e*((1.0-d)*(1.0-h1)+h1*g)+f*(1.0-g-e); + tria2[0]=4; tria2[1]= 5; tria2[2]= i2;//EFC + area2=(1.0-e)*(1.0-f); } /*}}}*/ } - else if (f11<=1.0-f22){ /*{{{*/ + else if (d<=1.0-g){ /*{{{*/ if (!istrapeze1 && !istrapeze2){ } else if (!istrapeze1 && istrapeze2){ - tria1[0]=i0; tria1[1]= 3; tria1[2]= 4; + tria1[0]=i0; tria1[1]= 3; tria1[2]= 4; //ADE + area1=d*e; } else if (istrapeze1 && !istrapeze2){ - tria1[0]=6; tria1[1]= i1; tria1[2]= 5; + tria1[0]=6; tria1[1]= i1; tria1[2]= 5; //GBF + area1=f*g; } else { //istrapeze1 && istrapeze2 - tria1[0]=3; tria1[1]= 6; tria1[2]= 4; - tria2[0]=4; tria2[1]= 6; tria2[2]= 5; - tria3[0]=4; tria3[1]= 5; tria3[2]= i2; + tria1[0]=3; tria1[1]= 6; tria1[2]= 4; //DGE + area1=e*(1.0-d-g); + tria2[0]=6; tria2[1]= i2; tria2[2]= 4; //GCE + area2=(1.0-g)*(1.0-e); + tria3[0]=6; tria3[1]= 5; tria3[2]= i2; //GFC + area3=g*(1.0-f); } /*}}}*/ } @@ -2157,57 +2213,81 @@ }/*}}}*/ else{ /*{{{*/ - if (f11<=f21 && f12>=f22){ /*{{{*/ + if (d<=f && e>=g){ /*{{{*/ if (!istrapeze1 && !istrapeze2){ - tria1[0]=i0; tria1[1]= 3; tria1[2]= 7; - tria2[0]=i0; tria2[1]= 7; tria2[2]= 6; + tria1[0]=i0; tria1[1]= 3; tria1[2]= 7; //ADH + area1=h1*d*e; + tria2[0]=i0; tria2[1]= 7; tria2[2]= 6; //AHG + area2=(1.0-h2)*f*g; } else if (!istrapeze1 && istrapeze2){ - tria1[0]=6; tria1[1]= 7; tria1[2]= 4; + tria1[0]=6; tria1[1]= 7; tria1[2]= 4; //GHE + area1=(e-g)*d*(1.0-h1); } else if (istrapeze1 && !istrapeze2){ - tria1[0]=3; tria1[1]= 5; tria1[2]= 7; + tria1[0]=3; tria1[1]= 5; tria1[2]= 7; //DFH + area1=(f-d)*g*h2; } else { //istrapeze1 && istrapeze2 - tria1[0]=7; tria1[1]= 5; tria1[2]= i1; - tria2[0]=7; tria2[1]= i1; tria2[2]= 4; - tria3[0]=4; tria3[1]= i1; tria3[2]= i2; + tria1[0]=5; tria1[1]= i1; tria1[2]= 7; //FBH + area1=g*h2*(1.0-f); + tria2[0]=7; tria2[1]= i1; tria2[2]= i2;//HBC + area2=1.0+d*(h1-e*h1-1.0)+g*h2*(d-1.0); + tria3[0]=7; tria3[1]= i2; tria3[2]= 4; //HCE + area3=d*(1.0-h1)*(1.0-e); } /*}}}*/ } - else if (f11>=f21 && f12<=f22){ /*{{{*/ + else if (d>=f && e<=g){ /*{{{*/ if (!istrapeze1 && !istrapeze2){ - tria1[0]=i0; tria1[1]= 5; tria1[2]= 7; - tria2[0]=i0; tria2[1]= 4; tria2[2]= 7; + tria1[0]=i0; tria1[1]= 5; tria1[2]= 7; //AFH + area1=h2*f*g; + tria2[0]=i0; tria2[1]= 7; tria2[2]= 4; //AHE + area2=(1.0-h1)*d*e; }else if (!istrapeze1 && istrapeze2){ - tria1[0]=5; tria1[1]= 3; tria1[2]= 7; + tria1[0]=5; tria1[1]= 3; tria1[2]= 7; //FDH + area1=(d-f)*e*h1; }else if (istrapeze1 && !istrapeze2){ - tria1[0]=4; tria1[1]= 7; tria1[2]= 6; + tria1[0]=4; tria1[1]= 7; tria1[2]= 6; //EHG + area1=(g-e)*f*(1.0-h2); }else { //istrapeze1 && istrapeze2 - tria1[0]=3; tria1[1]= i1; tria1[2]= 7; - tria2[0]=7; tria2[1]= i1; tria2[2]= 6; - tria3[0]=6; tria3[1]= i1; tria3[2]= i2; + tria1[0]=3; tria1[1]= i1; tria1[2]= 7; //DBH + area1=e*h1*(1.0-d); + tria2[0]=7; tria2[1]= i1; tria2[2]= 6; //HCG + area2=f*(1.0-h2)*(1.0-g); + tria3[0]=6; tria3[1]= i1; tria3[2]= i2; //HBC + area3=1.0+f*(h2-g*h2-1.0)+e*h1*(f-1.0); } /*}}}*/ } - else if (f11<=f21 && f12<=f22){ /*{{{*/ + else if (d<=f && e<=g){ /*{{{*/ if (!istrapeze1 && !istrapeze2){ - tria1[0]=i0; tria1[1]= 3; tria1[2]= 4; + tria1[0]=i0; tria1[1]= 3; tria1[2]= 4;//ADE + area1=d*e; }else if (!istrapeze1 && istrapeze2){ }else if (istrapeze1 && !istrapeze2){ - tria1[0]=3; tria1[1]= 5; tria1[2]= 4; - tria2[0]=4; tria2[1]= 5; tria2[2]= 6; + tria1[0]=3; tria1[1]= 5; tria1[2]= 6; //DFG + area1=g*(f-d); + tria2[0]=3; tria2[1]= 6; tria2[2]= 4; //DGE + area2=d*(g-e); }else { //istrapeze1 && istrapeze2 - tria1[0]=6; tria1[1]= 5; tria1[2]= i1; - tria2[0]=6; tria2[1]= i1; tria2[2]= i2; + tria1[0]=5; tria1[1]= i2; tria1[2]= 6; //FCG + area1=f*(1.0-g); + tria2[0]=5; tria2[1]= i1; tria2[2]= i2;//FBC + area2=1.0-f; } /*}}}*/ } - else if (f11>=f21 && f12>=f22){ /*{{{*/ + else if (d>=f && e>=g){ /*{{{*/ if (!istrapeze1 && !istrapeze2){ - tria1[0]=i0; tria1[1]= 5; tria1[2]= 6; + tria1[0]=i0; tria1[1]= 5; tria1[2]= 6; //AFG + area1=f*g; }else if (!istrapeze1 && istrapeze2){ - tria1[0]=5; tria1[1]= 3; tria1[2]= 6; - tria2[0]=6; tria2[1]= 3; tria2[2]= 4; + tria1[0]=5; tria1[1]= 4; tria1[2]= 6; //FEG + area1=f*(e-g); + tria2[0]=5; tria2[1]= 3; tria2[2]= 4; //FDE + area2=e*(d-f); }else if (istrapeze1 && !istrapeze2){ }else { //istrapeze1 && istrapeze2 - tria1[0]=4; tria1[1]= 3; tria1[2]= i1; - tria2[0]=i1; tria2[1]= i2; tria2[2]= 4; + tria1[0]=3; tria1[1]= i1; tria1[2]= i2; //DBC + area1=1.0-d; + tria2[0]=3; tria2[1]= i2; tria2[2]= 4;//DCE + area2=d*(1.0-e); } /*}}}*/ } @@ -2221,5 +2301,5 @@ } } - area1= GetTriangleAreaSpherical(xyz1); + area1*=area; //dimensionalize the fractional area from [0:1] to [0:area] } if(tria2[0]>-1){ @@ -2230,5 +2310,5 @@ } } - area2= GetTriangleAreaSpherical(xyz2); + area2*=area; } if(tria3[0]>-1){ @@ -2239,5 +2319,5 @@ } } - area3= GetTriangleAreaSpherical(xyz3); + area3*=area; } @@ -6100,9 +6180,9 @@ IssmDouble xyz_list[NUMVERTICES][3]; - /*stacks:*/ - IssmDouble* stackRSL = NULL; - IssmDouble* stackU = NULL; - IssmDouble* stackN = NULL; - IssmDouble* stackE = NULL; + /*viscous stacks:*/ + IssmDouble* viscousRSL = NULL; + IssmDouble* viscousU = NULL; + IssmDouble* viscousN = NULL; + IssmDouble* viscousE = NULL; #ifdef _HAVE_RESTRICT_ @@ -6112,5 +6192,5 @@ IssmDouble* __restrict__ GE=NULL; IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL; - IssmDouble* __restrict__ G_rigid_precomputed=NULL; + IssmDouble* __restrict__ G_gravi_precomputed=NULL; IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL; IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL; @@ -6121,5 +6201,5 @@ IssmDouble* GE=NULL; IssmDouble* G_viscoelastic_precomputed=NULL; - IssmDouble* G_rigid_precomputed=NULL; + IssmDouble* G_gravi_precomputed=NULL; IssmDouble* U_viscoelastic_precomputed=NULL; IssmDouble* H_viscoelastic_precomputed=NULL; @@ -6128,8 +6208,9 @@ /*viscoelastic green function:*/ int index; - int M; + int M; + IssmDouble doubleindex,lincoef; /*Computational flags:*/ - bool computerigid = false; + bool computeselfattraction = false; bool computeelastic = false; bool computerotation = false; @@ -6140,4 +6221,5 @@ IssmDouble start_time,final_time; int nt,precomputednt; + int grd, grdmodel; /*Rotational:*/ @@ -6155,7 +6237,8 @@ IssmDouble polenudge; /*}}}*/ + /*Recover parameters:{{{ */ rho_earth=FindParam(MaterialsEarthDensityEnum); - this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&computeselfattraction,SolidearthSettingsSelfAttractionEnum); this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); this->parameters->FindParam(&computerotation,SolidearthSettingsRotationEnum); @@ -6166,4 +6249,10 @@ this->parameters->FindParam(&NewtonG,ConstantsNewtonGravityEnum); this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + this->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); + this->parameters->FindParam(&grdmodel,GrdModelEnum); + + /*early return:*/ + if (!grd || grdmodel!=ElasticEnum) return; //Love numbers won't be found in this case, return before loading them + if(!computeselfattraction)return; if(computerotation){ @@ -6177,10 +6266,7 @@ /*}}}*/ - /*early return:*/ - if(!computerigid)return; - /*Recover precomputed green function kernels:{{{*/ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGSelfAttractionEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_gravi_precomputed,&M); if(computeelastic){ @@ -6215,5 +6301,5 @@ } else{ - nt=1; //in elastic, or if we run only rigid, we need only one step + nt=1; //in elastic, or if we run only selfattraction, we need only one step } @@ -6245,6 +6331,15 @@ delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda; alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); - index=reCast( alpha/M_PI*reCast(M-1) ); + doubleindex=alpha/M_PI*reCast(M-1); //maps 0(doubleindex); //truncates doubleindex to integer part _assert_(index>=0 && index(3*nt); - stackU=xNewZeroInit(3*nt); - - this->inputs->SetArrayInput(StackRSLEnum,this->lid,stackRSL,3*nt); - this->inputs->SetArrayInput(StackUEnum,this->lid,stackU,3*nt); - this->parameters->SetParam(0,StackIndexEnum); + viscousRSL=xNewZeroInit(3*nt); + viscousU=xNewZeroInit(3*nt); + + this->inputs->SetArrayInput(SealevelchangeViscousRSLEnum,this->lid,viscousRSL,3*nt); + this->inputs->SetArrayInput(SealevelchangeViscousUEnum,this->lid,viscousU,3*nt); + this->parameters->SetParam(0,SealevelchangeViscousIndexEnum); if(horiz){ - stackN=xNewZeroInit(3*nt); - stackE=xNewZeroInit(3*nt); - this->inputs->SetArrayInput(StackNEnum,this->lid,stackRSL,3*nt); - this->inputs->SetArrayInput(StackEEnum,this->lid,stackU,3*nt); + viscousN=xNewZeroInit(3*nt); + viscousE=xNewZeroInit(3*nt); + this->inputs->SetArrayInput(SealevelchangeViscousNEnum,this->lid,viscousRSL,3*nt); + this->inputs->SetArrayInput(SealevelchangeViscousEEnum,this->lid,viscousU,3*nt); } } @@ -6472,5 +6572,5 @@ slgeom->late[this->lid]=late; - /*compute fractional areas and load weights for ocean:*/ + /*compute areas and load weights for ocean and flag elements only partially in the ocean:*/ if(isoceanonly){ slgeom->LoadArea[SLGEOM_OCEAN][this->lid]=area; @@ -6582,5 +6682,5 @@ /*Deal with ice loads if we are on grounded ice:*/ if(isice && !isoceanonly && computeice){ - if(isiceonly){ + if(isiceonly && !isocean){ slgeom->LoadArea[SLGEOM_ICE][this->lid]=area; for(int i=0;iLoadWeigths[SLGEOM_ICE][i][this->lid]=1.0/3.0; @@ -6635,6 +6735,6 @@ IssmDouble loadweightsocean[3]; //to keep memory of these loads, no need to recompute for bottom pressure. IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble latbar=0; - IssmDouble longbar=0; + IssmDouble latbar=slgeom->late[this->lid]; + IssmDouble longbar=slgeom->longe[this->lid]; IssmDouble constant; IssmDouble nanconstant=NAN; @@ -6645,10 +6745,10 @@ #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticIceLatbarEnum,&nanconstant,P0Enum); - this->AddInput(SealevelBarystaticIceLongbarEnum,&nanconstant,P0Enum); - this->AddInput(SealevelBarystaticHydroLatbarEnum,&nanconstant,P0Enum); - this->AddInput(SealevelBarystaticHydroLongbarEnum,&nanconstant,P0Enum); - this->AddInput(SealevelBarystaticOceanLatbarEnum,&nanconstant,P0Enum); - this->AddInput(SealevelBarystaticOceanLongbarEnum,&nanconstant,P0Enum); + this->AddInput(SealevelBarystaticIceLatbarEnum,&latbar,P0Enum); + this->AddInput(SealevelBarystaticIceLongbarEnum,&longbar,P0Enum); + this->AddInput(SealevelBarystaticHydroLatbarEnum,&latbar,P0Enum); + this->AddInput(SealevelBarystaticHydroLongbarEnum,&longbar,P0Enum); + this->AddInput(SealevelBarystaticOceanLatbarEnum,&latbar,P0Enum); + this->AddInput(SealevelBarystaticOceanLongbarEnum,&longbar,P0Enum); #endif @@ -6725,5 +6825,5 @@ } /*}}}*/ -void Tria::SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){ /*{{{*/ +void Tria::SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){ /*{{{*/ /*Declarations:{{{*/ @@ -6740,5 +6840,5 @@ #ifdef _HAVE_RESTRICT_ IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL; - IssmDouble* __restrict__ G_rigid_precomputed=NULL; + IssmDouble* __restrict__ G_gravi_precomputed=NULL; IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL; IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL; @@ -6750,5 +6850,5 @@ #else IssmDouble* G_viscoelastic_precomputed=NULL; - IssmDouble* G_rigid_precomputed=NULL; + IssmDouble* G_gravi_precomputed=NULL; IssmDouble* U_viscoelastic_precomputed=NULL; IssmDouble* H_viscoelastic_precomputed=NULL; @@ -6759,13 +6859,15 @@ #endif - /*elastic green function:*/ + /*viscoelastic green function:*/ int index; - int M; + int M; + IssmDouble doubleindex,lincoef; /*Computational flags:*/ - bool computerigid = false; + bool computeselfattraction = false; bool computeelastic = false; bool computeviscous = false; int horiz; + int grd, grdmodel; bool istime=true; @@ -6777,5 +6879,5 @@ /*Recover parameters:{{{ */ rho_earth=FindParam(MaterialsEarthDensityEnum); - this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&computeselfattraction,SolidearthSettingsSelfAttractionEnum); this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum); @@ -6784,12 +6886,15 @@ this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + this->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); + this->parameters->FindParam(&grdmodel,GrdModelEnum); /*}}}*/ /*early return:*/ - if(!computerigid)return; + if (!grd || grdmodel!=ElasticEnum) return; //Love numbers won't be found in this case, return before loading them + if(!computeselfattraction)return; /*Recover precomputed green function kernels:{{{*/ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGSelfAttractionEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_gravi_precomputed,&M); if(computeelastic){ @@ -6826,5 +6931,5 @@ } else{ - nt=1; //in elastic, or if we run only rigid, we need only one step + nt=1; //in elastic, or if we run only selfattraction, we need only one step } Gsubel=xNew(SLGEOM_NUMLOADS); @@ -6884,6 +6989,12 @@ delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda; alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0))); - index=reCast( alpha/M_PI*reCast(M-1) ); - + doubleindex=alpha/M_PI*reCast(M-1); //maps 0(doubleindex); //truncates doubleindex to integer part + + lincoef=doubleindex-index; //where between index and index+1 is alpha + if (index==M-1){ //avoids out of bound case + index-=1; + lincoef=1; + } for(int t=0;tparameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - this->parameters->FindParam(&stacknumsteps,StackNumStepsEnum); - this->parameters->FindParam(&stacktimes,NULL,StackTimesEnum); - this->parameters->FindParam(&stackindex,StackIndexEnum); + this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum); + this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum); + this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum); this->parameters->FindParam(¤ttime,TimeEnum); - this->inputs->GetArrayPtr(StackRSLEnum,this->lid,&stackRSL,&dummy); - this->inputs->GetArrayPtr(StackUEnum,this->lid,&stackU,&dummy); + this->inputs->GetArrayPtr(SealevelchangeViscousRSLEnum,this->lid,&viscousRSL,&dummy); + this->inputs->GetArrayPtr(SealevelchangeViscousUEnum,this->lid,&viscousU,&dummy); if(horiz){ - this->inputs->GetArrayPtr(StackNEnum,this->lid,&stackN,&dummy); - this->inputs->GetArrayPtr(StackEEnum,this->lid,&stackE,&dummy); + this->inputs->GetArrayPtr(SealevelchangeViscousNEnum,this->lid,&viscousN,&dummy); + this->inputs->GetArrayPtr(SealevelchangeViscousEEnum,this->lid,&viscousE,&dummy); } @@ -6986,10 +7101,10 @@ int offset=1; lincoeff=0; - newindex=stacknumsteps-2; - - for(int t=stackindex;tcurrenttime){ + newindex=viscousnumsteps-2; + + for(int t=viscousindex;tcurrenttime){ newindex=t-1; - lincoeff=(currenttime-stacktimes[newindex])/(stacktimes[t]-stacktimes[newindex]); + lincoeff=(currenttime-viscoustimes[newindex])/(viscoustimes[t]-viscoustimes[newindex]); foundtime=true; offset=0; @@ -6999,21 +7114,21 @@ if(!foundtime) lincoeff=1; - stacktimes[newindex]=currenttime; + viscoustimes[newindex]=currenttime; for(int i=0;iparameters->SetParam(stackindex,StackIndexEnum); - this->parameters->SetParam(stacktimes,stacknumsteps,StackTimesEnum); + this->parameters->SetParam(viscousindex,SealevelchangeViscousIndexEnum); + this->parameters->SetParam(viscoustimes,viscousnumsteps,SealevelchangeViscousTimesEnum); /*free allocations:*/ - xDelete(stacktimes); + xDelete(viscoustimes); } @@ -7099,13 +7214,9 @@ } /*}}}*/ -void Tria::SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/ +void Tria::SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* polarmotionvector,SealevelGeometry* slgeom){ /*{{{*/ /*sal green function:*/ IssmDouble* G=NULL; IssmDouble* Gsub[SLGEOM_NUMLOADS]; - IssmDouble SealevelGRD[NUMVERTICES]={0}; - IssmDouble oceanaverage=0; - IssmDouble oceanarea=0; - IssmDouble rho_water; bool computefuture=false; @@ -7120,9 +7231,8 @@ IssmDouble Grotm3[3]; - this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&sal,SolidearthSettingsSelfAttractionEnum); this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum); this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); this->parameters->FindParam(&nel,MeshNumberofelementsEnum); - this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); if(sal){ @@ -7132,5 +7242,5 @@ this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size); - this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true,StackRSLEnum,computefuture=false); + this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false); } @@ -7142,5 +7252,5 @@ for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ - sealevelpercpu[this->vertices[i]->lid]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2]; + sealevelpercpu[this->vertices[i]->lid]+=Grotm1[i]*polarmotionvector[0]+Grotm2[i]*polarmotionvector[1]+Grotm3[i]*polarmotionvector[2]; } } @@ -7151,6 +7261,4 @@ /*sal green function:*/ - IssmDouble* G=NULL; - IssmDouble* Gsub[SLGEOM_NUMLOADS]; IssmDouble oceanaverage=0; IssmDouble oceanarea=0; @@ -7188,11 +7296,11 @@ return; } /*}}}*/ -void Tria::SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/ - - IssmDouble Sealevel[3]={0,0,0}; - IssmDouble SealevelRSL[3]={0,0,0}; - IssmDouble SealevelU[3]={0,0,0}; - IssmDouble SealevelN[3]={0,0,0}; - IssmDouble SealevelE[3]={0,0,0}; +void Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* polarmotionvector,SealevelGeometry* slgeom){ /*{{{*/ + + IssmDouble SealevelGrd[3]={0,0,0}; + IssmDouble RSLGrd[3]={0,0,0}; + IssmDouble UGrd[3]={0,0,0}; + IssmDouble NGrd[3]={0,0,0}; + IssmDouble EGrd[3]={0,0,0}; int nel,nbar; bool sal = false; @@ -7224,11 +7332,12 @@ bool elastic=false; bool percpu=false; + bool planethasocean=false; this->parameters->FindParam(&nel,MeshNumberofelementsEnum); - this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&sal,SolidearthSettingsSelfAttractionEnum); this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); this->parameters->FindParam(&elastic,SolidearthSettingsElasticEnum); this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - + this->parameters->FindParam(&planethasocean,SolidearthSettingsGrdOceanEnum); @@ -7257,11 +7366,11 @@ } } - this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel,percpu=false,StackRSLEnum,computefuture=true); + this->SealevelchangeGxL(&RSLGrd[0],G, Gsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true); if(elastic){ - this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel,percpu=false,StackUEnum,computefuture=true); + this->SealevelchangeGxL(&UGrd[0],GU, GUsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true); if(horiz ){ - this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel,percpu=false,StackNEnum,computefuture=true); - this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel,percpu=false,StackEEnum,computefuture=true); + this->SealevelchangeGxL(&NGrd[0],GN, GNsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true); + this->SealevelchangeGxL(&EGrd[0],GE, GEsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true); } } @@ -7275,5 +7384,5 @@ for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ - SealevelRSL[i]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2]; + RSLGrd[i]+=Grotm1[i]*polarmotionvector[0]+Grotm2[i]*polarmotionvector[1]+Grotm3[i]*polarmotionvector[2]; } } @@ -7286,5 +7395,5 @@ for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ - SealevelU[i]+=GUrotm1[i]*rotationvector[0]+GUrotm2[i]*rotationvector[1]+GUrotm3[i]*rotationvector[2]; + UGrd[i]+=GUrotm1[i]*polarmotionvector[0]+GUrotm2[i]*polarmotionvector[1]+GUrotm3[i]*polarmotionvector[2]; } } @@ -7299,6 +7408,6 @@ for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ - SealevelN[i]+=GNrotm1[i]*rotationvector[0]+GNrotm2[i]*rotationvector[1]+GNrotm3[i]*rotationvector[2]; - SealevelE[i]+=GErotm1[i]*rotationvector[0]+GErotm2[i]*rotationvector[1]+GErotm3[i]*rotationvector[2]; + NGrd[i]+=GNrotm1[i]*polarmotionvector[0]+GNrotm2[i]*polarmotionvector[1]+GNrotm3[i]*polarmotionvector[2]; + EGrd[i]+=GErotm1[i]*polarmotionvector[0]+GErotm2[i]*polarmotionvector[1]+GErotm3[i]*polarmotionvector[2]; } } @@ -7307,37 +7416,45 @@ } + if (planethasocean){ //We must also output the RSL on vertices to compute the ocean mass conservation + for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ + sealevelpercpu[this->vertices[i]->lid]=RSLGrd[i]; + } + } + } + /*Create geoid: */ - for(int i=0;iAddInput(SealevelGRDEnum,Sealevel,P1Enum); - this->AddInput(BedGRDEnum,SealevelU,P1Enum); + this->AddInput(SealevelGRDEnum,SealevelGrd,P1Enum); + this->AddInput(BedGRDEnum,UGrd,P1Enum); if(horiz){ - this->AddInput(BedNorthGRDEnum,SealevelN,P1Enum); - this->AddInput(BedEastGRDEnum,SealevelE,P1Enum); + this->AddInput(BedNorthGRDEnum,NGrd,P1Enum); + this->AddInput(BedEastGRDEnum,EGrd,P1Enum); } } /*}}}*/ -void Tria::SealevelchangeGxL(IssmDouble* sealevelout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu, int stackenum, bool computefuture) { /*{{{*/ - - IssmDouble* sealevel=NULL; +void Tria::SealevelchangeGxL(IssmDouble* grdfieldout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/ + + IssmDouble* grdfield=NULL; int i,e,l,t,nbar; - bool viscous=false; - IssmDouble* stack=NULL; + bool computeviscous=false; + IssmDouble* viscousfield=NULL; int nt=1; //important - int stackindex=0; //important - int stacknumsteps=1; //important - - this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum); - if(viscous){ - this->parameters->FindParam(&stacknumsteps,StackNumStepsEnum); - if(computefuture) nt=stacknumsteps; + int viscousindex=0; //important + int viscousnumsteps=1; //important + + this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum); + if(computeviscous){ + this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum); + if(computefuture) nt=viscousnumsteps; else nt=1; //allocate - sealevel=xNewZeroInit(3*nt); - } - else sealevel=xNewZeroInit(3*nt); + grdfield=xNewZeroInit(3*nt); + } + else grdfield=xNewZeroInit(3*nt); if(loads->sealevelloads){ @@ -7347,6 +7464,6 @@ for (e=0;esealevelloads[e]+loads->loads[e]); + int index=i*nel*viscousnumsteps+e*viscousnumsteps+t; + grdfield[i*nt+t]+=G[index]*(loads->sealevelloads[e]+loads->loads[e]); } } @@ -7355,6 +7472,6 @@ for (e=0;esubloads[l][e]); + int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t; + grdfield[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]); } } @@ -7362,6 +7479,6 @@ for (e=0;esubsealevelloads[e]); + int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t; + grdfield[i*nt+t]+=Gsub[l][index]*(loads->subsealevelloads[e]); } } @@ -7376,6 +7493,6 @@ for (e=0;eloads[e]); + int index=i*nel*viscousnumsteps+e*viscousnumsteps+t; + grdfield[i*nt+t]+=G[index]*(loads->loads[e]); } } @@ -7384,6 +7501,6 @@ for (e=0;esubloads[l][e]); + int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t; + grdfield[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]); } } @@ -7392,25 +7509,25 @@ } - if(viscous){ - IssmDouble* sealevelinterp=NULL; - IssmDouble* stacktimes=NULL; + if(computeviscous){ + IssmDouble* grdfieldinterp=NULL; + IssmDouble* viscoustimes=NULL; IssmDouble final_time; IssmDouble lincoeff; IssmDouble timeacc; - this->parameters->FindParam(&stackindex,StackIndexEnum); - this->parameters->FindParam(&stacktimes,NULL,StackTimesEnum); + this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum); + this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum); this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum); this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum); - this->inputs->GetArrayPtr(stackenum,this->lid,&stack,NULL); + this->inputs->GetArrayPtr(viscousenum,this->lid,&viscousfield,NULL); if(computefuture){ - sealevelinterp=xNew(3*nt); - if(stacktimes[stackindex](3*nt); + if(viscoustimes[viscousindex]lids[this->vertices[i]->lid]!=this->lid)continue; - sealevelinterp[i*nt+t]= sealevel[i*nt+0]; + grdfieldinterp[i*nt+t]= grdfield[i*nt+0]; } } @@ -7418,5 +7535,5 @@ for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; - sealevelinterp[i*nt+t]= (1-lincoeff)*sealevel[i*nt+(t-stackindex-1)]+lincoeff*sealevel[i*nt+(t-stackindex)]; + grdfieldinterp[i*nt+t]= (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]+lincoeff*grdfield[i*nt+(t-viscousindex)]; } } @@ -7425,24 +7542,25 @@ } - /*update sealevel at present time using stack at present time: */ + /*update grdfield at present time using viscous stack at present time: */ for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; - sealevel[i*nt+0]+=stack[i*stacknumsteps+stackindex]; - } - - if(computefuture){ /*update stack with future deformation from present load: */ - - for(int t=nt-1;t>=stackindex;t--){ + grdfield[i*nt+0]+=viscousfield[i*viscousnumsteps+viscousindex]; + } + + if(computefuture){ /*update viscous stack with future deformation from present load: */ + + for(int t=nt-1;t>=viscousindex;t--){ for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; - stack[i*stacknumsteps+t]+=sealevelinterp[i*nt+t]-sealevelinterp[i*nt+stackindex]-stack[i*stacknumsteps+stackindex]; + //offset viscousfield to remove all deformation that has already been added + viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*nt+t]-grdfieldinterp[i*nt+viscousindex]-viscousfield[i*viscousnumsteps+viscousindex]; } } - /*Re-add stack now that we updated:*/ - this->inputs->SetArrayInput(stackenum,this->lid,stack,3*stacknumsteps); + /*Re-add viscous stack now that we updated:*/ + this->inputs->SetArrayInput(viscousenum,this->lid,viscousfield,3*viscousnumsteps); } /*Free allocatoins:*/ - xDelete(stacktimes); + xDelete(viscoustimes); } @@ -7451,10 +7569,10 @@ for(i=0;ilids[this->vertices[i]->lid]==this->lid){ - sealevelout[this->vertices[i]->lid]=sealevel[i*nt+0]; + grdfieldout[this->vertices[i]->lid]=grdfield[i*nt+0]; } } } else{ - for(i=0;i* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae); - void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom); + void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom); void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae); void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae); @@ -176,5 +176,5 @@ void SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom); void SealevelchangeOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom); - void SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom); + void SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom); void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom); void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom); Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26295) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26296) @@ -24,6 +24,6 @@ void TransferSealevel(FemModel* femmodel,int forcingenum); bool slcconvergence(Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs); -IssmDouble SealevelloadsOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble oceanarea); -void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom); +IssmDouble SealevelloadsOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble totaloceanarea); +void PolarMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom); void ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom); void ivins_deformation_core(FemModel* femmodel); @@ -221,10 +221,10 @@ BarystaticContributions* barycontrib=NULL; GenericParam* barycontribparam=NULL; - IssmDouble rotationaxismotionvector[3]={0}; + IssmDouble polarmotionvector[3]={0}; GrdLoads* loads=NULL; Vector* oldsealevelloads=NULL; Vector* oceanareas=NULL; - IssmDouble oceanarea; + IssmDouble totaloceanarea; Vector* subelementoceanareas=NULL; IssmDouble oceanaverage; @@ -244,6 +244,8 @@ int grd=0; int grdmodel; - int computesealevel=0; + int sealevelloading=0; bool viscous=false; + bool rotation=false; + bool planethasocean=false; IssmDouble* sealevelpercpu=NULL; @@ -255,10 +257,14 @@ /*retrieve parameters:{{{*/ femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); + femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); - femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum); + femmodel->parameters->FindParam(&step,StepEnum); + femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->parameters->FindParam(&sealevelloading,SolidearthSettingsSealevelLoadingEnum); femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum); - femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); femmodel->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum); + femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); + femmodel->parameters->FindParam(&planethasocean,SolidearthSettingsGrdOceanEnum); /*}}}*/ @@ -314,11 +320,12 @@ loads->BroadcastLoads(); - //compute rotation axis motion: - RotationAxisMotion(&rotationaxismotionvector[0],femmodel,loads,slgeom); - - /*skip computation of sea level if requested, which means sea level loads should be zeroed */ - if(!computesealevel){ + //compute polar motion: + PolarMotion(&polarmotionvector[0],femmodel,loads,slgeom); + + /*skip computation of sea level equation if requested, which means sea level loads should be zeroed */ + if(!sealevelloading){ loads->sealevelloads=xNewZeroInit(nel); loads->subsealevelloads=xNewZeroInit(slgeom->nbar[SLGEOM_OCEAN]); + goto deformation; } @@ -332,5 +339,5 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeConvolution(sealevelpercpu, loads , rotationaxismotionvector,slgeom); + element->SealevelchangeConvolution(sealevelpercpu, loads , polarmotionvector,slgeom); } @@ -347,17 +354,17 @@ oceanareas->Assemble(); subelementoceanareas->Assemble(); - oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.); - if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 + oceanareas->Sum(&totaloceanarea); _assert_(totaloceanarea>0.); + if(scaleoceanarea) totaloceanarea=3.619e+14; // use true ocean area, m^2 } //Conserve ocean mass: - oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, oceanarea); - ConserveOceanMass(femmodel,loads,barycontrib->Total()/oceanarea - oceanaverage,slgeom); + oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, totaloceanarea); + ConserveOceanMass(femmodel,loads,barycontrib->Total()/totaloceanarea - oceanaverage,slgeom); //broadcast sea level loads loads->BroadcastSealevelLoads(); - //compute rotation axis motion: - RotationAxisMotion(&rotationaxismotionvector[0],femmodel,loads, slgeom); + //compute polar motion: + PolarMotion(&polarmotionvector[0],femmodel,loads, slgeom); //convergence? @@ -376,19 +383,46 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeDeformationConvolution(loads, rotationaxismotionvector,slgeom); + element->SealevelchangeDeformationConvolution(sealevelpercpu, loads, polarmotionvector,slgeom); } if(VerboseSolution()) _printf0_(" updating GRD fields\n"); - /*Update bedrock motion and geoid:*/ - if(computesealevel){ - femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water); //given that we converged, no need to recompute ocean average + if(planethasocean){ + if (!sealevelloading){ //we haven't done so before, so we need to compute the ocean average and area + loads->sealevelloads=NULL; //needed to trigger the calculation of areas + /*retrieve sea level average and ocean area:*/ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeOceanAverage(loads, oceanareas, subelementoceanareas, sealevelpercpu, slgeom); + } + loads->sealevelloads=xNewZeroInit(nel); + loads->AssembleSealevelLoads(); + /*compute ocean areas:*/ + oceanareas->Assemble(); + subelementoceanareas->Assemble(); + oceanareas->Sum(&totaloceanarea); _assert_(totaloceanarea>0.); + if(scaleoceanarea) totaloceanarea=3.619e+14; // use true ocean area, m^2 + + //Conserve ocean mass + //Note that here we create sea-level loads but they will not generate GRD as we have already run all the convolutions + oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, totaloceanarea); + ConserveOceanMass(femmodel,loads,barycontrib->Total()/totaloceanarea - oceanaverage,slgeom); + } + + femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/totaloceanarea- oceanaverage/rho_water); //cumulate barystatic contributions and save to results: barycontrib->Cumulate(femmodel->parameters); - barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea); + barycontrib->Save(femmodel->results,femmodel->parameters,totaloceanarea); barycontrib->Reset(); } + if (rotation) { + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,polarmotionvector[0],step,time)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,polarmotionvector[1],step,time)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,polarmotionvector[2],step,time)); + } + + /*Update bedrock motion and geoid:*/ femmodel->inputs->AXPY(1,SealevelGRDEnum,SealevelEnum); femmodel->inputs->AXPY(1,BedGRDEnum,BedEnum); @@ -409,5 +443,4 @@ /*parameters: */ int step; - int computesealevel=0; bool isocean=false; IssmDouble time; @@ -417,8 +450,4 @@ IssmDouble cumgmslc=0; IssmDouble gmtslc=0; - - /*early return if we are not computing sea level, but rather deformation: */ - femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum); - if (!computesealevel)return; /*early return if we have no ocean transport:*/ @@ -610,5 +639,5 @@ femmodel->parameters->FindParam(&areae,&nel,AreaeEnum); - /*initialize SealevelMasks structure: */ + /*initialize SealevelloadMasks structure: */ slgeom=new SealevelGeometry(femmodel->elements->Size(),femmodel->vertices->Size()); @@ -639,5 +668,5 @@ for(Object* & object : femmodel->elements->objects){ Element* element=xDynamicCast(object); - element->SealevelchangeGeometryFractionKernel(slgeom); + element->SealevelchangeGeometrySubElementKernel(slgeom); } @@ -699,5 +728,5 @@ } /*}}}*/ -IssmDouble SealevelloadsOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* suboceanareas, IssmDouble oceanarea){ /*{{{*/ +IssmDouble SealevelloadsOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* suboceanareas, IssmDouble totaloceanarea){ /*{{{*/ IssmDouble sealevelloadsaverage; @@ -715,7 +744,7 @@ delete vsubsealevelloadsvolume; - return (sealevelloadsaverage+subsealevelloadsaverage)/oceanarea; + return (sealevelloadsaverage+subsealevelloadsaverage)/totaloceanarea; } /*}}}*/ -void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom){ /*{{{*/ +void PolarMotion(IssmDouble* polarmotionvector, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom){ /*{{{*/ IssmDouble moi_list[3]={0,0,0}; @@ -746,5 +775,4 @@ element->SealevelchangeMomentOfInertiaCentroid(&moi_list[0],loads,slgeom); - element->SealevelchangeMomentOfInertiaSubElement(&moi_list_sub[0],loads, slgeom); @@ -754,4 +782,5 @@ } + for (int i=0;i<3;i++) moi_list[i]=0; ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); @@ -770,7 +799,7 @@ /*Assign output pointers:*/ - m[0]=m1; - m[1]=m2; - m[2]=m3; + polarmotionvector[0]=m1; + polarmotionvector[1]=m2; + polarmotionvector[2]=m3; } /*}}}*/ void ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26295) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26296) @@ -360,4 +360,7 @@ SolidearthSettingsViscousEnum, SealevelchangeGeometryDoneEnum, + SealevelchangeViscousNumStepsEnum, + SealevelchangeViscousTimesEnum, + SealevelchangeViscousIndexEnum, RotationalEquatorialMoiEnum, TidalLoveHEnum, @@ -370,7 +373,7 @@ LoveTimeFreqEnum, LoveIsTimeEnum, - SealevelchangeGRigidEnum, + SealevelchangeGSelfAttractionEnum, SealevelchangeGViscoElasticEnum, - SolidearthSettingsComputesealevelchangeEnum, + SolidearthSettingsSealevelLoadingEnum, SolidearthSettingsGRDEnum, SolidearthSettingsRunFrequencyEnum, @@ -379,9 +382,10 @@ SolidearthSettingsHorizEnum, SolidearthSettingsMaxiterEnum, + SolidearthSettingsGrdOceanEnum, SolidearthSettingsOceanAreaScalingEnum, RotationalPolarMoiEnum, SolidearthSettingsReltolEnum, SealevelchangeRequestedOutputsEnum, - SolidearthSettingsRigidEnum, + SolidearthSettingsSelfAttractionEnum, SolidearthSettingsRotationEnum, SolidearthSettingsMaxSHCoeffEnum, @@ -397,7 +401,4 @@ SettingsSolverResidueThresholdEnum, SettingsWaitonlockEnum, - StackNumStepsEnum, - StackTimesEnum, - StackIndexEnum, SmbAIceEnum, SmbAIdxEnum, @@ -832,4 +833,8 @@ SealevelchangeGEsubelHydroEnum, SealevelchangeGNsubelHydroEnum, + SealevelchangeViscousRSLEnum, + SealevelchangeViscousUEnum, + SealevelchangeViscousNEnum, + SealevelchangeViscousEEnum, SedimentHeadEnum, SedimentHeadOldEnum, @@ -953,8 +958,4 @@ SolidearthExternalDisplacementUpRateEnum, SolidearthExternalGeoidRateEnum, - StackRSLEnum, - StackUEnum, - StackNEnum, - StackEEnum, StrainRateeffectiveEnum, StrainRateparallelEnum, @@ -1430,4 +1431,5 @@ SealevelInertiaTensorYZEnum, SealevelInertiaTensorZZEnum, + SealevelchangePolarMotionEnum, SealevelNmotionEnum, SealevelUmotionEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26295) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26296) @@ -378,7 +378,7 @@ case LoveTimeFreqEnum : return "LoveTimeFreq"; case LoveIsTimeEnum : return "LoveIsTime"; - case SealevelchangeGRigidEnum : return "SealevelchangeGRigid"; + case SealevelchangeGSelfAttractionEnum : return "SealevelchangeGSelfAttraction"; case SealevelchangeGViscoElasticEnum : return "SealevelchangeGViscoElastic"; - case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange"; + case SolidearthSettingsSealevelLoadingEnum : return "SolidearthSettingsSealevelLoading"; case SolidearthSettingsGRDEnum : return "SolidearthSettingsGRD"; case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency"; @@ -391,5 +391,5 @@ case SolidearthSettingsReltolEnum : return "SolidearthSettingsReltol"; case SealevelchangeRequestedOutputsEnum : return "SealevelchangeRequestedOutputs"; - case SolidearthSettingsRigidEnum : return "SolidearthSettingsRigid"; + case SolidearthSettingsSelfAttractionEnum : return "SolidearthSettingsSelfAttraction"; case SolidearthSettingsRotationEnum : return "SolidearthSettingsRotation"; case SealevelchangeRunCountEnum : return "SealevelchangeRunCount"; @@ -404,7 +404,7 @@ case SettingsSolverResidueThresholdEnum : return "SettingsSolverResidueThreshold"; case SettingsWaitonlockEnum : return "SettingsWaitonlock"; - case StackNumStepsEnum : return "StackNumSteps"; - case StackTimesEnum : return "StackTimes"; - case StackIndexEnum : return "StackIndex"; + case SealevelchangeViscousNumStepsEnum : return "SealevelchangeViscousNumSteps"; + case SealevelchangeViscousTimesEnum : return "SealevelchangeViscousTimes"; + case SealevelchangeViscousIndexEnum : return "SealevelchangeViscousIndex"; case SmbAIceEnum : return "SmbAIce"; case SmbAIdxEnum : return "SmbAIdx"; @@ -957,8 +957,8 @@ case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate"; case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate"; - case StackRSLEnum : return "StackRSL"; - case StackUEnum : return "StackU"; - case StackNEnum : return "StackN"; - case StackEEnum : return "StackE"; + case SealevelchangeViscousRSLEnum : return "SealevelchangeViscousRSL"; + case SealevelchangeViscousUEnum : return "SealevelchangeViscousU"; + case SealevelchangeViscousNEnum : return "SealevelchangeViscousN"; + case SealevelchangeViscousEEnum : return "SealevelchangeViscousE"; case StrainRateeffectiveEnum : return "StrainRateeffective"; case StrainRateparallelEnum : return "StrainRateparallel"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26295) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26296) @@ -387,7 +387,7 @@ if(stage==4){ if (strcmp(name,"LoveIsTime")==0) return LoveIsTimeEnum; - else if (strcmp(name,"SealevelchangeGRigid")==0) return SealevelchangeGRigidEnum; + else if (strcmp(name,"SealevelchangeGSelfAttraction")==0) return SealevelchangeGSelfAttractionEnum; else if (strcmp(name,"SealevelchangeGViscoElastic")==0) return SealevelchangeGViscoElasticEnum; - else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum; + else if (strcmp(name,"SolidearthSettingsSealevelLoading")==0) return SolidearthSettingsSealevelLoadingEnum; else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum; else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; @@ -400,5 +400,5 @@ else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum; else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum; - else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum; + else if (strcmp(name,"SolidearthSettingsSelfAttraction")==0) return SolidearthSettingsSelfAttractionEnum; else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum; else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum; @@ -413,7 +413,7 @@ else if (strcmp(name,"SettingsSolverResidueThreshold")==0) return SettingsSolverResidueThresholdEnum; else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum; - else if (strcmp(name,"StackNumSteps")==0) return StackNumStepsEnum; - else if (strcmp(name,"StackTimes")==0) return StackTimesEnum; - else if (strcmp(name,"StackIndex")==0) return StackIndexEnum; + else if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum; + else if (strcmp(name,"SealevelchangeViscousTimes")==0) return SealevelchangeViscousTimesEnum; + else if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum; else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum; else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; @@ -978,8 +978,8 @@ else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum; else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum; - else if (strcmp(name,"StackRSL")==0) return StackRSLEnum; - else if (strcmp(name,"StackU")==0) return StackUEnum; - else if (strcmp(name,"StackN")==0) return StackNEnum; - else if (strcmp(name,"StackE")==0) return StackEEnum; + else if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum; + else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum; + else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum; + else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum; else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 26295) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 26296) @@ -9,11 +9,12 @@ abstol = 0; maxiter = 0; - rigid = 0; - elastic = 0; - viscous = 0; - rotation = 0; + selfattraction = 1; + elastic = 1; + viscous = 1; + rotation = 1; + grdocean = 1; ocean_area_scaling = 0; runfrequency = 1; %how many time steps we skip before we run grd_core - computesealevelchange = 1; %will sea-level be coputed? + sealevelloading = 1; %will sea-level loads be computed? isgrd = 0; %will GRD patterns be computed? compute_bp_grd = 0; %will GRD patterns for bottom pressures be computed? @@ -21,8 +22,29 @@ timeacc = 0; %time step accurary required to compute Green tables horiz = 0; %compute horizontal deformation - grdmodel = 0; %grd model (0 by default, 1 for elastic, 2 for Ivins) + grdmodel = 0; %grd model (0 by default, 1 for (visco)-elastic, 2 for Ivins) cross_section_shape = 0; %cross section only used when grd model is Ivins end + methods (Static) + function self = loadobj(self) % {{{ + % This function is directly called by matlab when a model object is + % loaded. If the input is a struct it is an old version of this class and + % old fields must be recovered (make sure they are in the deprecated + % model properties) + + if isstruct(self) + % 2021, Jun 4 + if isfield(self,'rigid') + self.selfattraction = self.rigid; + end + if isfield(self,'computesealevelchange') + self.sealevelloading = self.computesealevelchange; + end + self = structtoobj(solidearthsettings(),self); + + end + end% }}} + end methods + function self = solidearthsettings(varargin) % {{{ switch nargin @@ -43,12 +65,13 @@ %computational flags: - self.rigid=1; + self.selfattraction=1; self.elastic=1; self.viscous=1; self.rotation=1; + self.grdocean=1; self.ocean_area_scaling=0; self.compute_bp_grd=0; self.isgrd=0; - self.computesealevelchange=1; + self.sealevelloading=1; %numerical discretization accuracy @@ -83,8 +106,6 @@ md = checkfield(md,'fieldname','solidearth.settings.grdmodel','>=',0,'<=',2); md = checkfield(md,'fieldname','solidearth.settings.cross_section_shape','numel',[1],'values',[1,2]); - - %checks on computational flags - if self.elastic==1 & self.rigid==0, - error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set'); + if self.elastic==1 & self.selfattraction==0, + error('solidearthsettings checkconsistency error message: need selfattraction on if elastic flag is set'); end if self.viscous==1 & self.elastic==0, @@ -103,4 +124,7 @@ end end + if self.sealevelloading==1 & self.grdocean==0 + error('solidearthsettings checkconsistency error message: need grdocean on if sealevelloading flag is set'); + end end @@ -116,16 +140,17 @@ fielddisplay(self,'abstol','sea level change absolute convergence criterion, NaN: not applied'); fielddisplay(self,'maxiter','maximum number of nonlinear iterations'); + fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module [default: 1]'); fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); - fielddisplay(self,'computesealevelchange','compute sealevel change (default 1)'); + fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default 1)'); fielddisplay(self,'isgrd','compute GRD patterns (default 1)'); fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default 1)'); fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)'); - fielddisplay(self,'rigid','rigid earth graviational potential perturbation'); - fielddisplay(self,'elastic','elastic earth graviational potential perturbation'); - fielddisplay(self,'viscous','viscous earth graviational potential perturbation'); - fielddisplay(self,'rotation','earth rotational potential perturbation'); + fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field'); + fielddisplay(self,'elastic','enables elastic deformation from surface loading'); + fielddisplay(self,'viscous','enables viscous deformation from surface loading'); + fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields'); fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'); fielddisplay(self,'timeacc','time accuracy (default 1 yr) for numerical discretization of the Green''s functions'); - fielddisplay(self,'grdmodel','type of deformation model, 1 for elastic, 2 for visco-elastic from Ivins'); + fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)'); fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore'); end % }}} @@ -134,8 +159,9 @@ WriteData(fid,prefix,'object',self,'fieldname','abstol','name','md.solidearth.settings.abstol','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','maxiter','name','md.solidearth.settings.maxiter','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','rigid','name','md.solidearth.settings.rigid','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','selfattraction','name','md.solidearth.settings.selfattraction','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','elastic','name','md.solidearth.settings.elastic','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','viscous','name','md.solidearth.settings.viscous','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','rotation','name','md.solidearth.settings.rotation','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','grdocean','name','md.solidearth.settings.grdocean','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','name','md.solidearth.settings.ocean_area_scaling','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','runfrequency','name','md.solidearth.settings.runfrequency','format','Integer'); @@ -143,5 +169,5 @@ WriteData(fid,prefix,'object',self,'fieldname','timeacc','name','md.solidearth.settings.timeacc','format','Double','scale',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','horiz','name','md.solidearth.settings.horiz','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','computesealevelchange','name','md.solidearth.settings.computesealevelchange','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','sealevelloading','name','md.solidearth.settings.sealevelloading','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','isgrd','name','md.solidearth.settings.isgrd','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','compute_bp_grd','name','md.solidearth.settings.compute_bp_grd','format','Integer'); @@ -154,8 +180,9 @@ writejsdouble(fid,[modelname '.solidearth.settings.reltol'],self.reltol); writejsdouble(fid,[modelname '.solidearth.settings.abstol'],self.abstol); - writejsdouble(fid,[modelname '.solidearth.settings.rigid'],self.rigid); + writejsdouble(fid,[modelname '.solidearth.settings.selfattraction'],self.selfattraction); writejsdouble(fid,[modelname '.solidearth.settings.elastic'],self.elastic); writejsdouble(fid,[modelname '.solidearth.settings.viscous'],self.viscous); writejsdouble(fid,[modelname '.solidearth.settings.rotation'],self.rotation); + writejsdouble(fid,[modelname '.solidearth.settings.grdocean'],self.rotation); writejsdouble(fid,[modelname '.solidearth.settings.ocean_area_scaling'],self.ocean_area_scaling); writejsdouble(fid,[modelname '.solidearth.settings.run_frequency'],self.run_frequency); Index: /issm/trunk-jpl/test/NightlyRun/test2001.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2001.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2001.m (revision 26296) @@ -17,5 +17,6 @@ md.initialization.sealevel=zeros(md.mesh.numberofvertices,1); md.solidearth.settings.cross_section_shape=1; % for square-edged x-section -md.solidearth.settings.computesealevelchange=0; %do not compute sea level, only deformation +md.solidearth.settings.grdocean=0; %do not compute sea level, only deformation +md.solidearth.settings.sealevelloading=0; %do not compute sea level, only deformation md.solidearth.requested_outputs={'Sealevel','BedGRD'}; Index: /issm/trunk-jpl/test/NightlyRun/test2002.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 26296) @@ -3,7 +3,8 @@ %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %700 km resolution mesh +load ../Data/SlcTestMesh.mat; +md.mesh=SlcMesh; %700 km resolution mesh -%Geometry for the bed, arbitrary thickness of 1000: +%Geometry for the bed, arbitrary thickness of 100: md.geometry.bed=zeros(md.mesh.numberofvertices,1); md.geometry.base=md.geometry.bed; @@ -70,5 +71,5 @@ md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; -md.solidearth.settings.computesealevelchange=1; +md.solidearth.settings.sealevelloading=1; md.solidearth.settings.isgrd=1; md.solidearth.settings.ocean_area_scaling=0; @@ -87,5 +88,5 @@ %eustatic run: -md.solidearth.settings.rigid=0; +md.solidearth.settings.selfattraction=0; md.solidearth.settings.elastic=0; md.solidearth.settings.rotation=0; @@ -96,15 +97,15 @@ Beustatic=md.results.TransientSolution.Bed; -%eustatic + rigid run: -md.solidearth.settings.rigid=1; +%eustatic + selfattraction run: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=0; md.solidearth.settings.rotation=0; md.solidearth.settings.viscous=0; md=solve(md,'tr'); -Srigid=md.results.TransientSolution.Sealevel; -Brigid=md.results.TransientSolution.Bed; +Sselfattraction=md.results.TransientSolution.Sealevel; +Bselfattraction=md.results.TransientSolution.Bed; -%eustatic + rigid + elastic run: -md.solidearth.settings.rigid=1; +%eustatic + selfattraction + elastic run: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=0; @@ -114,6 +115,6 @@ Belastic=md.results.TransientSolution.Bed; -%eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid=1; +%eustatic + selfattraction + elastic + rotation run: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; @@ -124,5 +125,5 @@ %Fields and tolerances to track changes -field_names={'Seustatic','Srigid','Selastic','Srotation','Beustatic','Brigid','Belastic','Brotation'}; +field_names={'Seustatic','Sselfattraction','Selastic','Srotation','Beustatic','Bselfattraction','Belastic','Brotation'}; field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; -field_values={Seustatic,Srigid,Selastic,Srotation,Beustatic,Brigid,Belastic,Brotation}; +field_values={Seustatic,Sselfattraction,Selastic,Srotation,Beustatic,Bselfattraction,Belastic,Brotation}; Index: /issm/trunk-jpl/test/NightlyRun/test2003.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2003.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2003.m (revision 26296) @@ -3,5 +3,6 @@ %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %1000 km resolution mesh +load ../Data/SlcTestMesh.mat; +md.mesh=SlcMesh; %700 km resolution mesh %Geometry for the bed, arbitrary thickness of 1000: @@ -62,5 +63,6 @@ md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; -md.solidearth.settings.computesealevelchange=0; +md.solidearth.settings.sealevelloading=0; +md.solidearth.settings.grdocean=0; md.solidearth.settings.isgrd=1; md.solidearth.settings.ocean_area_scaling=0; @@ -68,5 +70,4 @@ md.solidearth.settings.horiz=1; md.solidearth.requested_outputs={'Sealevel','Bed', 'BedEast', 'BedNorth'}; - %Physics: @@ -81,6 +82,6 @@ md.timestepping.final_time=1; -%eustatic + rigid + elastic run: -md.solidearth.settings.rigid=1; +%eustatic + selfattraction + elastic run: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=0; @@ -94,6 +95,6 @@ BNnoRotation=md.results.TransientSolution.BedNorth; -%eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid=1; +%eustatic + selfattraction + elastic + rotation run: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; @@ -108,5 +109,5 @@ %Fields and tolerances to track changes -field_names ={'noRotation','Rotation'}; +field_names ={'Sealevel', 'Uplift', 'NorthDisplacement', 'EastDisplacement'}; field_tolerances={1e-13,1e-13,1e-13,1e-13}; field_values={SRotation-SnoRotation,BURotation-BUnoRotation,BNRotation-BNnoRotation,BERotation-BEnoRotation}; Index: /issm/trunk-jpl/test/NightlyRun/test2004.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 26296) @@ -395,5 +395,5 @@ md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; -md.solidearth.settings.computesealevelchange=1; +md.solidearth.settings.sealevelloading=1; md.solidearth.settings.isgrd=1; md.solidearth.settings.ocean_area_scaling=0; @@ -423,5 +423,5 @@ %eustatic run: -md.solidearth.settings.rigid=0; +md.solidearth.settings.selfattraction=0; md.solidearth.settings.elastic=0; md.solidearth.settings.rotation=0; @@ -432,13 +432,13 @@ Seustatic=md.results.TransientSolution.Sealevel; -%eustatic + rigid run: -md.solidearth.settings.rigid=1; +%eustatic + selfattraction run: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=0; md.solidearth.settings.rotation=0; md=solve(md,'Transient'); -Srigid=md.results.TransientSolution.Sealevel; - -%eustatic + rigid + elastic run: -md.solidearth.settings.rigid=1; +Sselfattraction=md.results.TransientSolution.Sealevel; + +%eustatic + selfattraction + elastic run: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=0; @@ -446,6 +446,6 @@ Selastic=md.results.TransientSolution.Sealevel; -%eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid=1; +%eustatic + selfattraction + elastic + rotation run: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; @@ -457,3 +457,3 @@ field_names ={'Eustatic','Rigid','Elastic','Rotation'}; field_tolerances={1e-13,1e-13,1e-13,1e-13}; -field_values={Seustatic,Srigid,Selastic,Srotation}; +field_values={Seustatic,Sselfattraction,Selastic,Srotation}; Index: /issm/trunk-jpl/test/NightlyRun/test2005.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2005.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2005.m (revision 26296) @@ -1,51 +1,63 @@ -%Test Name: EarthSlc +%Test Name: Earthslc_time %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh +load ../Data/SlcTestMesh.mat; +md.mesh=SlcMesh; %700 km resolution mesh -%parameterize solidearth solution: +%Geometry for the bed, arbitrary thickness of 100: +md.geometry.bed=zeros(md.mesh.numberofvertices,1); +md.geometry.base=md.geometry.bed; +md.geometry.thickness=100*ones(md.mesh.numberofvertices,1); +md.geometry.surface=md.geometry.bed+md.geometry.thickness; + %solidearth loading: {{{ -md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); -md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); -md.dsl.global_average_thermosteric_sea_level_change=[1;0]; -md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); -md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +md.masstransport.spcthickness=[md.geometry.thickness;0]; +md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); %antarctica -late=sum(md.mesh.lat(md.mesh.elements),2)/3; -longe=sum(md.mesh.long(md.mesh.elements),2)/3; +xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3; +ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3; +ze=md.mesh.z(md.mesh.elements)*[1;1;1]/3; +re=sqrt(xe.^2+ye.^2+ze.^2); + +late=asind(ze./re); +longe=atan2d(ye,xe); pos=find(late < -80); -md.solidearth.surfaceload.icethicknesschange(pos)=-100; +md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; +posant=pos; %greenland -pos=find(late>70 & late<80 & longe>-60 & longe<-30); -md.solidearth.surfaceload.icethicknesschange(pos)=-100; +pos=find(late>60 & late<90 & longe>-75 & longe<-15); +md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; +posgre=pos; %elastic loading from love numbers: -md.solidearth.lovenumbers=lovenumbers('maxdeg',100); +md.solidearth.lovenumbers=lovenumbers('maxdeg',1000); %}}} %mask: {{{ mask=gmtmask(md.mesh.lat,md.mesh.long); +oceanmask=-ones(md.mesh.numberofvertices,1); +pos=find(mask==0); oceanmask(pos)=1; + icemask=ones(md.mesh.numberofvertices,1); -pos=find(mask==0); -icemask(pos)=-1; -pos=find(sum(mask(md.mesh.elements),2)<3); -icemask(md.mesh.elements(pos,:))=-1; +icemask(md.mesh.elements(posant))=-1; +icemask(md.mesh.elements(posgre))=-1; + md.mask.ice_levelset=icemask; -md.mask.ocean_levelset=-icemask; - -%make sure that the elements that have loads are fully grounded: -pos=find(md.solidearth.surfaceload.icethicknesschange); -md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; - -%make sure wherever there is an ice load, that the mask is set to ice: -%pos=find(md.solidearth.surfaceload.icethicknesschange); % TODO: Do we need to do this twice? -md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +md.mask.ocean_levelset=oceanmask; % }}} -md.solidearth.settings.ocean_area_scaling=0; +%time stepping: +md.timestepping.start_time=0; +md.timestepping.time_step=1; +md.timestepping.final_time=10; -%Geometry for the bed, arbitrary: -md.geometry.bed=-ones(md.mesh.numberofvertices,1); +%masstransport: +md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.initialization.vx=zeros(md.mesh.numberofvertices,1); +md.initialization.vy=zeros(md.mesh.numberofvertices,1); +md.initialization.sealevel=zeros(md.mesh.numberofvertices,1); +md.initialization.str=0; %Materials: @@ -56,32 +68,32 @@ %Solution parameters +md.cluster.np=3; md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; -md.solidearth.settings.computesealevelchange=1; +md.solidearth.settings.sealevelloading=1; +md.solidearth.settings.isgrd=1; +md.solidearth.settings.ocean_area_scaling=0; +md.solidearth.settings.grdmodel=1; -% max number of iteration reverted back to 10 (i.e., the original default value) -md.solidearth.settings.maxiter=10; - -%eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid=1; +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; +md.solidearth.settings.viscous=0; -%transient settings: -md.timestepping.start_time=0; -md.timestepping.final_time=10; -md.timestepping.time_step=1; -md.transient.isslc=1; -md.transient.issmb=0; -md.transient.ismasstransport=0; +%Physics: +md.transient.issmb=0; md.transient.isstressbalance=0; md.transient.isthermal=0; -dh=md.solidearth.surfaceload.icethicknesschange; -deltathickness=zeros(md.mesh.numberofelements+1,10); -for i=1:10, - deltathickness(1:end-1,i)=dh*i; +md.transient.ismasstransport=1; +md.transient.isslc=1; +md.solidearth.requested_outputs={'Sealevel'}; + +dh=md.masstransport.spcthickness; +deltathickness=zeros(md.mesh.numberofvertices+1,10); +for i=0:10, + deltathickness(1:end-1,i+1)=md.geometry.thickness+dh(1:end-1)*i; end -deltathickness(end,:)=0:1:9; -md.solidearth.surfaceload.icethicknesschange=deltathickness; +deltathickness(end,:)=0:1:10; +md.masstransport.spcthickness=deltathickness; %hack: md.geometry.surface=zeros(md.mesh.numberofvertices,1); @@ -89,4 +101,5 @@ md.geometry.base=-ones(md.mesh.numberofvertices,1); md.geometry.bed=md.geometry.base; + %run transient solution: Index: /issm/trunk-jpl/test/NightlyRun/test2006.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2006.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2006.m (revision 26296) @@ -3,22 +3,31 @@ %mesh earth: md=model; -md.cluster=generic('name',oshostname(),'np',5); -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh - -%parameterize solidearth solution: +load ../Data/SlcTestMesh.mat; +md.mesh=SlcMesh; %700 km resolution mesh + +%Geometry for the bed, arbitrary thickness of 100: +md.geometry.bed=zeros(md.mesh.numberofvertices,1); +md.geometry.base=md.geometry.bed; +md.geometry.thickness=100*ones(md.mesh.numberofvertices,1); +md.geometry.surface=md.geometry.bed+md.geometry.thickness; + %solidearth loading: {{{ -md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); -md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); -md.dsl.global_average_thermosteric_sea_level_change=[1;0]; -md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); -md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +md.masstransport.spcthickness=[md.geometry.thickness;0]; +md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); %antarctica -late=sum(md.mesh.lat(md.mesh.elements),2)/3; -longe=sum(md.mesh.long(md.mesh.elements),2)/3; +xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3; +ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3; +ze=md.mesh.z(md.mesh.elements)*[1;1;1]/3; +re=sqrt(xe.^2+ye.^2+ze.^2); + +late=asind(ze./re); +longe=atan2d(ye,xe); pos=find(late < -80); -md.solidearth.surfaceload.icethicknesschange(pos)=-100; +md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; +posant=pos; %greenland -pos=find(late>70 & late<80 & longe>-60 & longe<-30); -md.solidearth.surfaceload.icethicknesschange(pos)=-100; +pos=find(late>60 & late<90 & longe>-75 & longe<-15); +md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; +posgre=pos; %elastic loading from love numbers: @@ -28,25 +37,27 @@ %mask: {{{ mask=gmtmask(md.mesh.lat,md.mesh.long); +oceanmask=-ones(md.mesh.numberofvertices,1); +pos=find(mask==0); oceanmask(pos)=1; + icemask=ones(md.mesh.numberofvertices,1); -pos=find(mask==0); -icemask(pos)=-1; -pos=find(sum(mask(md.mesh.elements),2)<3); -icemask(md.mesh.elements(pos,:))=-1; +icemask(md.mesh.elements(posant))=-1; +icemask(md.mesh.elements(posgre))=-1; + md.mask.ice_levelset=icemask; -md.mask.ocean_levelset=-icemask; - -%make sure that the elements that have loads are fully grounded: -pos=find(md.solidearth.surfaceload.icethicknesschange); -md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; - -%make sure wherever there is an ice load, that the mask is set to ice: -pos=find(md.solidearth.surfaceload.icethicknesschange); -md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +md.mask.ocean_levelset=oceanmask; % }}} -md.solidearth.settings.ocean_area_scaling=0; - -%Geometry for the bed, arbitrary: -md.geometry.bed=-ones(md.mesh.numberofvertices,1); +%time stepping: +md.timestepping.start_time=0; +md.timestepping.time_step=1; +md.timestepping.final_time=10; + +%masstransport: +md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.initialization.vx=zeros(md.mesh.numberofvertices,1); +md.initialization.vy=zeros(md.mesh.numberofvertices,1); +md.initialization.sealevel=zeros(md.mesh.numberofvertices,1); +md.initialization.str=0; %Materials: @@ -57,32 +68,32 @@ %Solution parameters +md.cluster.np=3; md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; -md.solidearth.settings.computesealevelchange=1; - -% max number of iteration reverted back to 10 (i.e., the original default value) -md.solidearth.settings.maxiter=10; - -%eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid=1; +md.solidearth.settings.sealevelloading=1; +md.solidearth.settings.isgrd=1; +md.solidearth.settings.ocean_area_scaling=0; +md.solidearth.settings.grdmodel=1; + +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; - -%transient settings: -md.timestepping.start_time=0; -md.timestepping.final_time=10; -md.timestepping.time_step=1; -md.transient.isslc=1; -md.transient.issmb=0; -md.transient.ismasstransport=0; +md.solidearth.settings.viscous=0; + +%Physics: +md.transient.issmb=0; md.transient.isstressbalance=0; md.transient.isthermal=0; -dh=md.solidearth.surfaceload.icethicknesschange; -deltathickness=zeros(md.mesh.numberofelements+1,10); -for i=1:10, - deltathickness(1:end-1,i)=dh*i; +md.transient.ismasstransport=1; +md.transient.isslc=1; +md.solidearth.requested_outputs={'Sealevel'}; + +dh=md.masstransport.spcthickness; +deltathickness=zeros(md.mesh.numberofvertices+1,10); +for i=0:10, + deltathickness(1:end-1,i+1)=md.geometry.thickness+dh(1:end-1)*i; end -deltathickness(end,:)=0:1:9; -md.solidearth.surfaceload.icethicknesschange=deltathickness; +deltathickness(end,:)=0:1:10; +md.masstransport.spcthickness=deltathickness; %hack: Index: /issm/trunk-jpl/test/NightlyRun/test2007.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2007.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2007.m (revision 26296) @@ -1,9 +1,10 @@ -%Test Name: offline solid earth solution +%Test Name: External_OfflineSolidearthSolution %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh +load ../Data/SlcTestMesh.mat; +md.mesh=SlcMesh; %700 km resolution mesh -%Geometry for the bed, arbitrary thickness of 1000: +%Geometry for the bed, arbitrary md.geometry.bed=-ones(md.mesh.numberofvertices,1); md.geometry.base=md.geometry.bed; @@ -46,4 +47,5 @@ %Solution parameters md.cluster.np=3; +md.solidearth.settings.isgrd=0; md.solidearth.settings.horiz=1; Index: /issm/trunk-jpl/test/NightlyRun/test2008.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2008.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2008.m (revision 26296) @@ -1,7 +1,8 @@ -%Test Name: External:AdditionalSolidearthSolution +%Test Name: External_AdditionalSolidearthSolution %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh +load ../Data/SlcTestMesh.mat; +md.mesh=SlcMesh; %700 km resolution mesh %Geometry for the bed, arbitrary thickness of 100: @@ -70,5 +71,5 @@ md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; -md.solidearth.settings.computesealevelchange=1; +md.solidearth.settings.sealevelloading=1; md.solidearth.settings.isgrd=1; md.solidearth.settings.ocean_area_scaling=0; @@ -86,6 +87,6 @@ md.solidearth.settings.maxiter=10; -%eustatic + rigid + elastic: -md.solidearth.settings.rigid=1; +%eustatic + selfattraction + elastic: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=0; Index: /issm/trunk-jpl/test/NightlyRun/test2010.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2010.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2010.m (revision 26296) @@ -3,25 +3,34 @@ %mesh earth: md=model; -rad_e = 6.371012*10^3; % mean radius of Earth, km -md.mesh=gmshplanet('radius',rad_e,'resolution',1000.0); % km resolution +load ../Data/SlcTestMesh.mat; +md.mesh=SlcMesh; %700 km resolution mesh + +%Geometry for the bed, arbitrary thickness of 1000: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); +md.geometry.base=md.geometry.bed; +md.geometry.thickness=100*ones(md.mesh.numberofvertices,1); +md.geometry.surface=md.geometry.bed+md.geometry.thickness; + %parameterize slc solution: -%slc loading: {{{ -late=sum(md.mesh.lat(md.mesh.elements),2)/3; -longe=sum(md.mesh.long(md.mesh.elements),2)/3; +%solidearth loading: {{{ +md.masstransport.spcthickness=[md.geometry.thickness;0]; +md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); -md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); -pos=find(late <-75 & longe >0); -md.solidearth.surfaceload.icethicknesschange(pos(6:7))=-1; -md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); -md.dsl.global_average_thermosteric_sea_level_change=[0;0]; -md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); -md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3; +ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3; +ze=md.mesh.z(md.mesh.elements)*[1;1;1]/3; +re=sqrt(xe.^2+ye.^2+ze.^2); -md.solidearth.settings.ocean_area_scaling = 1; +late=asind(ze./re); +longe=atan2d(ye,xe); +%greenland +pos=find(late>60 & late<90 & longe>-75 & longe<-15); +md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; +posice=pos; %elastic loading from love numbers: -md.solidearth.lovenumbers=lovenumbers('maxdeg',1000); +md.solidearth.lovenumbers=lovenumbers('maxdeg',100); %}}} @@ -29,68 +38,89 @@ mask=gmtmask(md.mesh.lat,md.mesh.long); icemask=ones(md.mesh.numberofvertices,1); -pos=find(mask==0); -icemask(pos)=-1; -pos=find(sum(mask(md.mesh.elements),2)<3); -icemask(md.mesh.elements(pos,:))=-1; +icemask(md.mesh.elements(posice,:))=-0.5; + +oceanmask=-ones(md.mesh.numberofvertices,1); +pos=find(mask==0); oceanmask(pos)=1; icemask(~pos)=1; + md.mask.ice_levelset=icemask; -md.mask.ocean_levelset=-icemask; +md.mask.ocean_levelset=oceanmask; -%make sure that the elements that have loads are fully grounded: -pos=find(md.solidearth.surfaceload.icethicknesschange); -md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; +% use model representation of ocen area (not the true area) +md.solidearth.settings.ocean_area_scaling = 0; -%make sure wherever there is an ice load, that the mask is set to ice: -md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); -pos=find(md.solidearth.surfaceload.icethicknesschange); -md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; -% }}} +%materials +md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); +md.initialization.sealevel=zeros(md.mesh.numberofvertices,1); +md.initialization.str=0; -%geometry {{{ -di=md.materials.rho_ice/md.materials.rho_water; -md.geometry.thickness=ones(md.mesh.numberofvertices,1); -md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); -md.geometry.base=md.geometry.surface-md.geometry.thickness; -md.geometry.bed=md.geometry.base; -% }}} -%materials {{{ -md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); -md.materials.rheology_B=paterson(md.initialization.temperature); -md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); -% }}} -%Miscellaneous {{{ +md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.initialization.vx=zeros(md.mesh.numberofvertices,1); +md.initialization.vy=zeros(md.mesh.numberofvertices,1); + +%Miscellaneous md.miscellaneous.name='test2010'; -% }}} -%Solution parameters {{{ + +%Solution parameters md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; -md.solidearth.settings.computesealevelchange=1; -% }}} +md.solidearth.settings.sealevelloading=0; +md.solidearth.settings.grdocean=1; +md.solidearth.settings.isgrd=1; +md.solidearth.settings.ocean_area_scaling=0; +md.solidearth.settings.grdmodel=1; +md.solidearth.settings.horiz=1; +md.solidearth.requested_outputs={'Sealevel','SealevelBarystaticIceArea','SealevelBarystaticIceLoad','SealevelBarystaticIceMask','SealevelBarystaticIceLatbar' 'SealevelBarystaticIceLongbar'}; -%eustatic + rigid + elastic run: -md.solidearth.settings.rigid=1; +%Physics: +md.transient.issmb=0; +md.transient.isstressbalance=0; +md.transient.isthermal=0; +md.transient.ismasstransport=1; +md.transient.isslc=1; + +md.timestepping.start_time=0; +md.timestepping.time_step=1; +md.timestepping.final_time=1; + +%eustatic + selfattraction + elastic + rotation run: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; +md.solidearth.settings.viscous=0; md.cluster=generic('name',oshostname(),'np',3); +%md.verbose=verbose('111111111'); +md=solve(md,'Transient'); +moi_p=md.solidearth.rotational.polarmoi; +moi_e=md.solidearth.rotational.equatorialmoi; +tide_love_k2=md.solidearth.lovenumbers.tk(3); +load_love_k2=md.solidearth.lovenumbers.k(3); +tide_love_k2secular=md.solidearth.lovenumbers.tk2secular; % uncomment following 2 lines for -md=solve(md,'Sealevelchange'); -eus=md.results.SealevelchangeSolution.Bslc; -slc=md.results.SealevelchangeSolution.Sealevel; -moixz=md.results.SealevelchangeSolution.SealevelInertiaTensorXZ; -moiyz=md.results.SealevelchangeSolution.SealevelInertiaTensorYZ; -moizz=md.results.SealevelchangeSolution.SealevelInertiaTensorZZ; +eus=md.results.TransientSolution.Bslc; +slc=md.results.TransientSolution.Sealevel; +moixz=md.results.TransientSolution.SealevelInertiaTensorXZ / (1/(1-tide_love_k2/tide_love_k2secular) * (1+load_love_k2)/(moi_p-moi_e) ); +moiyz=md.results.TransientSolution.SealevelInertiaTensorYZ / (1/(1-tide_love_k2/tide_love_k2secular) * (1+load_love_k2)/(moi_p-moi_e) ); +moizz=md.results.TransientSolution.SealevelInertiaTensorZZ / ( -(1+load_love_k2)/moi_p); + +areaice=md.results.TransientSolution.SealevelBarystaticIceArea; +loadice=md.results.TransientSolution.SealevelBarystaticIceLoad; % analytical moi => just checking FOR ICE only!!! {{{ % ...have to mute ** slc induced MOI in Tria.cpp ** prior to the comparison -%rad_e = rad_e*1e3; % now in meters -%areas=GetAreasSphericalTria(md.mesh.elements,md.mesh.lat,md.mesh.long,rad_e); -%lat=late*pi/180; lon=longe*pi/180; -%moi_xz = sum(-md.materials.rho_freshwater.*md.solidearth.surfaceload.icethicknesschange.*areas.*rad_e^2.*sin(lat).*cos(lat).*cos(lon)); -%moi_yz = sum(-md.materials.rho_freshwater.*md.solidearth.surfaceload.icethicknesschange.*areas.*rad_e^2.*sin(lat).*cos(lat).*sin(lon)); +rad_e = md.solidearth.planetradius; + +lat=md.results.TransientSolution.SealevelBarystaticIceLatbar*pi/180; +lon=md.results.TransientSolution.SealevelBarystaticIceLongbar*pi/180; +moi_xz = sum(-loadice.*areaice.*rad_e^2.*sin(lat).*cos(lat).*cos(lon)); +moi_yz = sum(-loadice.*areaice.*rad_e^2.*sin(lat).*cos(lat).*sin(lon)); +moi_zz = sum(-loadice.*areaice.*rad_e^2.*(1.0-sin(lat).^2)); +theoretical_value_check=[moixz/moi_xz moiyz/moi_yz moizz/moizz] % }}} %Fields and tolerances to track changes field_names ={'eus','slc','moixz','moiyz','moizz'}; -field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13}; +field_tolerances={1e-13,1e-13,1e-13,1e-13}; field_values={eus,slc,moixz,moiyz,moizz}; Index: /issm/trunk-jpl/test/NightlyRun/test2011.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2011.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2011.m (revision 26296) @@ -1,7 +1,8 @@ -%Test Name: EarthSlc +%Test Name: SlcBarystatic %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %700 km resolution mesh +load ../Data/SlcTestMesh.mat; +md.mesh=SlcMesh; %700 km resolution mesh %Geometry for the bed, arbitrary thickness of 1000: @@ -70,5 +71,5 @@ %Miscellaneous -md.miscellaneous.name='test2002'; +md.miscellaneous.name='test2011'; %Solution parameters @@ -76,5 +77,5 @@ md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; -md.solidearth.settings.computesealevelchange=1; +md.solidearth.settings.sealevelloading=1; md.solidearth.settings.isgrd=1; md.solidearth.settings.ocean_area_scaling=0; @@ -94,7 +95,8 @@ %eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid=1; +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; -md.solidearth.settings.rotation=1; +md.solidearth.settings.rotation=0; +md.solidearth.settings.viscous=0; md=solve(md,'tr'); @@ -106,5 +108,5 @@ weights=md.results.TransientSolution(1).SealevelBarystaticIceWeights; dH=md.results.TransientSolution(1).DeltaIceThickness(md.mesh.elements); -dHavg=sum(dH.*weights); +dHavg=sum(dH.*weights,2); oceanarea=sum(md.results.TransientSolution(1).SealevelBarystaticOceanArea); bslc2=-sum(dHavg.*area)*md.materials.rho_ice/md.materials.rho_water/oceanarea; @@ -113,3 +115,3 @@ field_names={'BarystaticIce','BarystaticIce2','BarystaticIceDiff'}; field_tolerances={1e-13,1e-13,1e-13}; -field_values={bslc,bslc2,bslc2-bslc2}; +field_values={bslc,bslc2,bslc2-bslc}; Index: /issm/trunk-jpl/test/NightlyRun/test2090.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2090.m (revision 26295) +++ /issm/trunk-jpl/test/NightlyRun/test2090.m (revision 26296) @@ -3,5 +3,6 @@ %mesh earth: md=model; -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh +load ../Data/SlcTestMesh.mat; +md.mesh=SlcMesh; %700 km resolution mesh %Geometry for the bed, arbitrary thickness of 1000: @@ -73,5 +74,5 @@ time1=md.timestepping.start_time:md.timestepping.time_step:md.timestepping.final_time; md.masstransport.spcthickness=repmat(md.masstransport.spcthickness, [1 length(time1)]); -md.masstransport.spcthickness(1:end-1,2:end)=0; +md.masstransport.spcthickness(1:end-1,3:end)=0; md.masstransport.spcthickness(end,:)=time1; @@ -92,13 +93,13 @@ %Solution parameters -md.cluster.np=1; +md.cluster.np=3; md.solidearth.settings.reltol=NaN; md.solidearth.settings.abstol=1e-3; -md.solidearth.settings.computesealevelchange=1; +md.solidearth.settings.sealevelloading=1; md.solidearth.settings.isgrd=1; md.solidearth.settings.ocean_area_scaling=0; md.solidearth.settings.grdmodel=1; md.solidearth.settings.timeacc=md.timestepping.time_step; -md.solidearth.settings.degacc=.01; +md.solidearth.settings.degacc=.1; %Physics: @@ -113,6 +114,6 @@ md.solidearth.settings.maxiter=10; -%eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid=1; +%eustatic + selfattraction + elastic + rotation run: +md.solidearth.settings.selfattraction=1; md.solidearth.settings.elastic=1; md.solidearth.settings.viscous=1; @@ -127,4 +128,7 @@ end +hold on +plot(md.timestepping.start_time:md.timestepping.time_step:(md.timestepping.final_time-md.timestepping.time_step), [B(916,:)]); + %Fields and tolerances to track changes field_names={'Sealevel','Bed'};