Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 26226) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 26227) @@ -559,4 +559,5 @@ ./cores/sealevelchange_core.cpp \ ./analyses/SealevelchangeAnalysis.cpp\ + ./classes/GrdLoads.cpp\ ./classes/SealevelGeometry.cpp Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26226) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26227) @@ -31,7 +31,7 @@ class DatasetInput; class IoModel; -class SealevelMasks; class SealevelGeometry; class Gauss; +class GrdLoads; class ElementVector; template class Matrix; @@ -391,13 +391,13 @@ virtual void SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom)=0; - virtual void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom)=0; - virtual void SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, SealevelGeometry* slgeom)=0; + virtual void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0; + virtual void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0; virtual void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0; virtual void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0; virtual void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0; - virtual void SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom)=0; - virtual void SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0; - virtual void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0; - virtual void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom)=0; + virtual void SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom)=0; + virtual void SealevelchangeConvolution(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0; + virtual void SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0; + virtual void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26226) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26227) @@ -224,13 +224,13 @@ 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 SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, 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");}; void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");}; - void SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeConvolution(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26226) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26227) @@ -179,13 +179,13 @@ 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 SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, 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");}; void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");}; - void SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeConvolution(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26226) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26227) @@ -186,13 +186,13 @@ 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 SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, 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");}; void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");}; - void SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeConvolution(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");}; + void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26226) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26227) @@ -6806,5 +6806,5 @@ } /*}}}*/ -void Tria::SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/ +void Tria::SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/ /*Inputs:*/ @@ -6863,19 +6863,19 @@ if(slgeom->issubelement[SLGEOM_ICE][this->lid]){ int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid]; - subelementiceloads->SetValue(intj,Iavg,INS_VAL); + loads->vsubloads[SLGEOM_ICE]->SetValue(intj,Iavg,INS_VAL); Iavg=0; //avoid double counting centroid loads and subelement loads! } if(slgeom->issubelement[SLGEOM_WATER][this->lid]){ int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid]; - subelementhydroloads->SetValue(intj,Wavg,INS_VAL); + loads->vsubloads[SLGEOM_WATER]->SetValue(intj,Wavg,INS_VAL); Wavg=0; } if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; - subelementbploads->SetValue(intj,BPavg,INS_VAL); + loads->vsubloads[SLGEOM_OCEAN]->SetValue(intj,BPavg,INS_VAL); BPavg=0; } /*Plug remaining values into centroi load vector:*/ - loads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL); + loads->vloads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL); /*Keep track of barystatic contributions:*/ @@ -6884,11 +6884,9 @@ } /*}}}*/ -void Tria::SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/ +void Tria::SealevelchangeConvolution(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/ /*sal green function:*/ IssmDouble* G=NULL; - IssmDouble* GsubelIce=NULL; - IssmDouble* GsubelHydro=NULL; - IssmDouble* GsubelOcean=NULL; + IssmDouble* Gsub[SLGEOM_NUMLOADS]; IssmDouble SealevelGRD[NUMVERTICES]={0}; IssmDouble oceanaverage=0; @@ -6919,46 +6917,9 @@ if(sal){ this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); - this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&GsubelIce,&size); - this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&GsubelHydro,&size); - this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&GsubelOcean,&size); - - if(allsealevelloads){ - for(int i=0;inbar[SLGEOM_ICE]; - for (int e=0;enbar[SLGEOM_WATER]; - for (int e=0;enbar[SLGEOM_OCEAN]; - for (int e=0;enbar[SLGEOM_ICE]; - for (int e=0;enbar[SLGEOM_WATER]; - for (int e=0;enbar[SLGEOM_OCEAN]; - for (int e=0;einputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size); + + this->SealevelchangeGxL(&SealevelGRD[0],G, Gsub, loads, slgeom, nel); } @@ -6975,9 +6936,9 @@ if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; - subelementsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL); - } - else sealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL); - - if(!allsealevelloads){ + loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL); + } + else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL); + + if(!loads->sealevelloads){ oceanareas->SetValue(this->sid,oceanarea,INS_VAL); if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ @@ -6990,5 +6951,5 @@ return; } /*}}}*/ -void Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/ +void Tria::SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/ IssmDouble Sealevel[3]={0,0,0}; @@ -7003,16 +6964,8 @@ IssmDouble* GE=NULL; IssmDouble* GN=NULL; - IssmDouble* GsubelIce=NULL; - IssmDouble* GsubelHydro=NULL; - IssmDouble* GsubelOcean=NULL; - IssmDouble* GUsubelIce=NULL; - IssmDouble* GUsubelHydro=NULL; - IssmDouble* GUsubelOcean=NULL; - IssmDouble* GNsubelIce=NULL; - IssmDouble* GNsubelHydro=NULL; - IssmDouble* GNsubelOcean=NULL; - IssmDouble* GEsubelIce=NULL; - IssmDouble* GEsubelHydro=NULL; - IssmDouble* GEsubelOcean=NULL; + IssmDouble* Gsub[SLGEOM_NUMLOADS]; + IssmDouble* GUsub[SLGEOM_NUMLOADS]; + IssmDouble* GNsub[SLGEOM_NUMLOADS]; + IssmDouble* GEsub[SLGEOM_NUMLOADS]; int horiz; @@ -7039,79 +6992,35 @@ if(sal){ + this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); - this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&GsubelIce,&size); - this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&GsubelHydro,&size); - this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&GsubelOcean,&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size); + this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size); + if(elastic){ this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size); - this->inputs->GetArrayPtr(SealevelchangeGUsubelIceEnum,this->lid,&GUsubelIce,&size); - this->inputs->GetArrayPtr(SealevelchangeGUsubelHydroEnum,this->lid,&GUsubelHydro,&size); - this->inputs->GetArrayPtr(SealevelchangeGUsubelOceanEnum,this->lid,&GUsubelOcean,&size); + this->inputs->GetArrayPtr(SealevelchangeGUsubelIceEnum,this->lid,&GUsub[SLGEOM_ICE],&size); + this->inputs->GetArrayPtr(SealevelchangeGUsubelHydroEnum,this->lid,&GUsub[SLGEOM_WATER],&size); + this->inputs->GetArrayPtr(SealevelchangeGUsubelOceanEnum,this->lid,&GUsub[SLGEOM_OCEAN],&size); if(horiz){ + this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size); + this->inputs->GetArrayPtr(SealevelchangeGNsubelIceEnum,this->lid,&GNsub[SLGEOM_ICE],&size); + this->inputs->GetArrayPtr(SealevelchangeGNsubelHydroEnum,this->lid,&GNsub[SLGEOM_WATER],&size); + this->inputs->GetArrayPtr(SealevelchangeGNsubelOceanEnum,this->lid,&GNsub[SLGEOM_OCEAN],&size); + this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size); - this->inputs->GetArrayPtr(SealevelchangeGEsubelIceEnum,this->lid,&GEsubelIce,&size); - this->inputs->GetArrayPtr(SealevelchangeGEsubelOceanEnum,this->lid,&GEsubelOcean,&size); - this->inputs->GetArrayPtr(SealevelchangeGEsubelHydroEnum,this->lid,&GEsubelHydro,&size); - this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size); - this->inputs->GetArrayPtr(SealevelchangeGNsubelIceEnum,this->lid,&GNsubelIce,&size); - this->inputs->GetArrayPtr(SealevelchangeGNsubelHydroEnum,this->lid,&GNsubelHydro,&size); - this->inputs->GetArrayPtr(SealevelchangeGNsubelOceanEnum,this->lid,&GNsubelOcean,&size); - } - } - - for(int i=0;inbar[SLGEOM_ICE]; - for (int e=0;enbar[SLGEOM_WATER]; - for (int e=0;enbar[SLGEOM_OCEAN]; - for (int e=0;enbar[SLGEOM_ICE]; - for (int e=0;enbar[SLGEOM_WATER]; - for (int e=0;enbar[SLGEOM_OCEAN]; - for (int e=0;enbar[SLGEOM_ICE]; - for (int e=0;enbar[SLGEOM_WATER]; - for (int e=0;enbar[SLGEOM_OCEAN]; - for (int e=0;einputs->GetArrayPtr(SealevelchangeGEsubelIceEnum,this->lid,&GEsub[SLGEOM_ICE],&size); + this->inputs->GetArrayPtr(SealevelchangeGEsubelHydroEnum,this->lid,&GEsub[SLGEOM_WATER],&size); + this->inputs->GetArrayPtr(SealevelchangeGEsubelOceanEnum,this->lid,&GEsub[SLGEOM_OCEAN],&size); + } + } + + this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel); + + if(elastic){ + this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel); + if(horiz ){ + this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel); + this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel); } } @@ -7131,5 +7040,40 @@ } /*}}}*/ -void Tria::SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){/*{{{*/ +void Tria::SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel) { /*{{{*/ + + if(loads->sealevelloads){ + for(int i=0;isealevelloads[e]+loads->loads[e]); + } + for(int l=0;lnbar[l]; + for (int e=0;esubloads[l][e]); + } + if(l==SLGEOM_OCEAN){ + for (int e=0;esubsealevelloads[e]); + } + } + } + } + } + else{ //this is the initial convolution where only loads are provided + for(int i=0;iloads[e]); + } + for(int l=0;lnbar[l]; + for (int e=0;esubloads[l][e]); + } + } + } + } + +} /*}}}*/ +void Tria::SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){/*{{{*/ IssmDouble S=0; @@ -7148,6 +7092,6 @@ /*recover total load: */ - if(loads) S+=loads[this->Sid()]; - if(sealevelloads) S+=sealevelloads[this->Sid()]; + if(loads->loads) S+=loads->loads[this->Sid()]; + if(loads->sealevelloads) S+=loads->sealevelloads[this->Sid()]; /* Perturbation terms for moment of inertia (moi_list): @@ -7161,5 +7105,5 @@ return; }/*}}}*/ -void Tria::SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){/*{{{*/ +void Tria::SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* grdloads, SealevelGeometry* slgeom){/*{{{*/ IssmDouble SA=0; @@ -7181,16 +7125,6 @@ for(int i=0;iissubelement[i][this->lid]){ - switch(i){ - case SLGEOM_ICE: - loads=subelementiceloads; - break; - case SLGEOM_WATER: - loads=subelementhydroloads; - break; - case SLGEOM_OCEAN: - loads=subelementbploads; - sealevelloads=subelementsealevelloads; - break; - } + loads=grdloads->subloads[i]; + if(i==SLGEOM_OCEAN) sealevelloads=grdloads->subsealevelloads; intj=slgeom->subelementmapping[i][this->lid]; late=slgeom->latbarycentre[i][intj]/180*M_PI; @@ -7215,12 +7149,12 @@ return; }/*}}}*/ -void Tria::SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/ +void Tria::SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/ if(slgeom->isoceanin[this->lid]){ if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; - subelementloads->SetValue(intj,offset,ADD_VAL); - } - else loads->SetValue(this->sid,offset,ADD_VAL); + loads->vsubsealevelloads->SetValue(intj,offset,ADD_VAL); + } + else loads->vsealevelloads->SetValue(this->sid,offset,ADD_VAL); } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26226) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26227) @@ -173,10 +173,10 @@ void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae); void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae); - void SealevelchangeBarystaticLoads(Vector* loads, Vector* subelementiceloads, Vector* subelementhydroloads, Vector* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom); - void SealevelchangeConvolution(Vector* sealevelloads, Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom); - void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom); - void SealevelchangeShift(Vector* loads, Vector* subelementloads, IssmDouble offset, SealevelGeometry* slgeom); - void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom); - void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom); + void SealevelchangeBarystaticLoads(GrdLoads* loads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom); + void SealevelchangeConvolution(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom); + void SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom); + void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom); + void SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom); + void SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom); #endif /*}}}*/ @@ -214,9 +214,9 @@ Gauss* NewGauss(void); Gauss* NewGauss(int order); - Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order); - Gauss* NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order); - Gauss* NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order); - Gauss* NewGauss(IssmDouble fraction1,IssmDouble fraction2,int order); - Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert); + Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order); + Gauss* NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order); + Gauss* NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order); + Gauss* NewGauss(IssmDouble fraction1,IssmDouble fraction2,int order); + Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert); Gauss* NewGaussBase(int order); Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");}; @@ -237,9 +237,10 @@ void SetTemporaryElementType(int element_type_in){_error_("not implemented yet");}; void InputServe(Input* input_in); - Seg* SpawnSeg(int index1,int index2); + Seg* SpawnSeg(int index1,int index2); IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; void StabilizationParameterAnisotropic(IssmDouble* tau_parameter_ansiotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){_error_("not implemented yet");}; void UpdateConstraintsExtrudeFromBase(void); void UpdateConstraintsExtrudeFromTop(void); + void SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/GrdLoads.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/GrdLoads.cpp (revision 26227) +++ /issm/trunk-jpl/src/c/classes/GrdLoads.cpp (revision 26227) @@ -0,0 +1,71 @@ +/* + * \file GrdLoads.cpp + * \brief: Implementation of GrdLoads class + */ + +/*Headers: {{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./GrdLoads.h" +#include "./SealevelGeometry.h" +using namespace std; +/*}}}*/ + +/*Object constructors and destructor*/ +GrdLoads::GrdLoads(int nel,SealevelGeometry* slgeom){ /*{{{*/ + + + vloads=new Vector(nel); + for (int i=0;i(slgeom->nbar[i]); + + vsealevelloads=new Vector(nel); + vsealevelloads->Set(0);vsealevelloads->Assemble(); + + vsubsealevelloads=new Vector(slgeom->nbar[SLGEOM_OCEAN]); + + +}; /*}}}*/ +GrdLoads::~GrdLoads(){ /*{{{*/ + + delete vloads; + xDelete(loads); + for(int i=0;i(subloads[i]); + } + delete vsealevelloads; + xDelete(sealevelloads); + delete vsubsealevelloads; + xDelete(subsealevelloads); +}; /*}}}*/ + +void GrdLoads::BroadcastLoads(void){ /*{{{*/ + + /*Initialize barycentre vectors, now that we know their size: */ + vloads->Assemble(); + for (int i=0;iAssemble(); + } + + loads=vloads->ToMPISerial(); + for (int i=0;iToMPISerial(); + } + +} /*}}}*/ +void GrdLoads::AssembleSealevelLoads(void){ /*{{{*/ + + vsealevelloads->Assemble(); + vsubsealevelloads->Assemble(); + +} /*}}}*/ +void GrdLoads::BroadcastSealevelLoads(void){ /*{{{*/ + + sealevelloads=vsealevelloads->ToMPISerial(); + subsealevelloads=vsubsealevelloads->ToMPISerial(); + +} /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/GrdLoads.h =================================================================== --- /issm/trunk-jpl/src/c/classes/GrdLoads.h (revision 26227) +++ /issm/trunk-jpl/src/c/classes/GrdLoads.h (revision 26227) @@ -0,0 +1,32 @@ +/*!\file GrdLoads.h + * \brief: header file for GrdLoads + */ + +#ifndef _SEALEVELGRDLOADS_H_ +#define _SEALEVELGRDLOADS_H_ + +/*Headers:*/ +#include "./SealevelGeometry.h" +#include "../toolkits/toolkits.h" + +class GrdLoads{ + + public: + + Vector* vloads=NULL; + IssmDouble* loads=NULL; + Vector* vsubloads[SLGEOM_NUMLOADS]; + IssmDouble* subloads[SLGEOM_NUMLOADS]; + Vector* vsealevelloads=NULL; + IssmDouble* sealevelloads=NULL; + Vector* vsubsealevelloads=NULL; + IssmDouble* subsealevelloads=NULL; + + GrdLoads(int nel, SealevelGeometry* slgeom); + ~GrdLoads(); + + void AssembleSealevelLoads(void); + void BroadcastLoads(void); + void BroadcastSealevelLoads(void); +}; +#endif /* _SEALEVELGRDLOADS_H_ */ Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 26226) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 26227) @@ -19,4 +19,5 @@ #include "./Misfit.h" #include "./SealevelGeometry.h" +#include "./GrdLoads.h" #include "./BarystaticContributions.h" #include "./Nodalvalue.h" Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26226) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26227) @@ -23,8 +23,8 @@ void TransferForcing(FemModel* femmodel,int forcingenum); void TransferSealevel(FemModel* femmodel,int forcingenum); -void slcconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs); -IssmDouble SealevelloadsOceanAverage(Vector* sealevelloads,Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble oceanarea); -void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom); -void ConserveOceanMass(FemModel* femmodel,Vector* sealevelloads, Vector* subelementsealevelloads, IssmDouble offset, SealevelGeometry* slgeom); +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); +void ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom); void ivins_deformation_core(FemModel* femmodel); IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel); @@ -32,5 +32,5 @@ /*main cores:*/ -void sealevelchange_core(FemModel* femmodel){ /*{{{*/ +void sealevelchange_core(FemModel* femmodel){ /*{{{*/ SealevelGeometry* slgeom=NULL; @@ -83,5 +83,5 @@ } /*}}}*/ -void solidearthexternal_core(FemModel* femmodel){ /*{{{*/ +void solidearthexternal_core(FemModel* femmodel){ /*{{{*/ /*variables:*/ @@ -152,5 +152,5 @@ } /*}}}*/ -void couplerinput_core(FemModel* femmodel){ /*{{{*/ +void couplerinput_core(FemModel* femmodel){ /*{{{*/ /*Be very careful here, everything is well thought through, do not remove @@ -215,5 +215,5 @@ }; /*}}}*/ -void grd_core(FemModel* femmodel, SealevelGeometry* slgeom) { /*{{{*/ +void grd_core(FemModel* femmodel, SealevelGeometry* slgeom) { /*{{{*/ /*variables:{{{*/ @@ -223,16 +223,5 @@ IssmDouble rotationaxismotionvector[3]={0}; - Vector* loads=NULL; - IssmDouble* allloads=NULL; - Vector* subelementiceloads=NULL; - IssmDouble* allsubelementiceloads=NULL; - Vector* subelementhydroloads=NULL; - IssmDouble* allsubelementhydroloads=NULL; - Vector* subelementbploads=NULL; - IssmDouble* allsubelementbploads=NULL; - Vector* sealevelloads=NULL; - IssmDouble* allsealevelloads=NULL; - Vector* subelementsealevelloads=NULL; - IssmDouble* allsubelementsealevelloads=NULL; + GrdLoads* loads=NULL; Vector* oldsealevelloads=NULL; Vector* oceanareas=NULL; @@ -249,5 +238,4 @@ int step; IssmDouble time; - bool converged=false; int modelid,earthid; @@ -301,12 +289,5 @@ femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum); - loads=new Vector(nel); - subelementiceloads=new Vector(slgeom->nbar[SLGEOM_ICE]); - subelementhydroloads=new Vector(slgeom->nbar[SLGEOM_WATER]); - subelementbploads=new Vector(slgeom->nbar[SLGEOM_OCEAN]); - - sealevelloads=new Vector(nel); - sealevelloads->Set(0);sealevelloads->Assemble(); - subelementsealevelloads=new Vector(slgeom->nbar[SLGEOM_OCEAN]); + loads=new GrdLoads(nel,slgeom); subelementoceanareas=new Vector(slgeom->nbar[SLGEOM_OCEAN]); oceanareas=new Vector(nel); @@ -317,25 +298,17 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeBarystaticLoads(loads, subelementiceloads, subelementhydroloads, subelementbploads, barycontrib,slgeom); - } - loads->Assemble(); - subelementiceloads->Assemble(); - subelementhydroloads->Assemble(); - subelementbploads->Assemble(); - + element->SealevelchangeBarystaticLoads(loads, barycontrib,slgeom); + } //broadcast loads - allloads=loads->ToMPISerial(); - allsubelementiceloads=subelementiceloads->ToMPISerial(); - allsubelementhydroloads=subelementhydroloads->ToMPISerial(); - allsubelementbploads=subelementbploads->ToMPISerial(); + loads->BroadcastLoads(); //compute rotation axis motion: - RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsubelementiceloads,allsubelementhydroloads,allsubelementbploads, NULL, NULL,slgeom); + RotationAxisMotion(&rotationaxismotionvector[0],femmodel,loads,slgeom); /*skip computation of sea level if requested, which means sea level loads should be zeroed */ if(!computesealevel){ - allsealevelloads=xNewZeroInit(nel); - allsubelementsealevelloads=xNewZeroInit(slgeom->nbar[SLGEOM_OCEAN]); + loads->sealevelloads=xNewZeroInit(nel); + loads->subsealevelloads=xNewZeroInit(slgeom->nbar[SLGEOM_OCEAN]); goto deformation; } @@ -344,17 +317,16 @@ for(;;){ - oldsealevelloads=sealevelloads->Duplicate(); sealevelloads->Copy(oldsealevelloads); + oldsealevelloads=loads->vsealevelloads->Duplicate(); loads->vsealevelloads->Copy(oldsealevelloads); /*convolve load and sealevel loads on oceans:*/ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeConvolution(sealevelloads, subelementsealevelloads, oceanareas, subelementoceanareas, allsealevelloads, allloads,allsubelementiceloads, allsubelementhydroloads, allsubelementbploads, allsubelementsealevelloads, rotationaxismotionvector,slgeom); - } - - sealevelloads->Assemble(); - subelementsealevelloads->Assemble(); - + element->SealevelchangeConvolution(loads, oceanareas, subelementoceanareas, rotationaxismotionvector,slgeom); + } + + loads->AssembleSealevelLoads(); + /*compute ocean areas:*/ - if(!allsealevelloads){ //first time in the loop + if(!loads->sealevelloads){ //first time in the loop oceanareas->Assemble(); subelementoceanareas->Assemble(); @@ -364,22 +336,19 @@ //Conserve ocean mass: - oceanaverage=SealevelloadsOceanAverage(sealevelloads,subelementsealevelloads, oceanareas,subelementoceanareas, oceanarea); - - ConserveOceanMass(femmodel,sealevelloads,subelementsealevelloads,barycontrib->Total()/oceanarea - oceanaverage,slgeom); + oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, oceanarea); + ConserveOceanMass(femmodel,loads,barycontrib->Total()/oceanarea - oceanaverage,slgeom); //broadcast sea level loads - allsealevelloads=sealevelloads->ToMPISerial(); - allsubelementsealevelloads=subelementsealevelloads->ToMPISerial(); + loads->BroadcastSealevelLoads(); //compute rotation axis motion: - RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsubelementiceloads,allsubelementhydroloads,allsubelementbploads, allsealevelloads,allsubelementsealevelloads,slgeom); + RotationAxisMotion(&rotationaxismotionvector[0],femmodel,loads, slgeom); //convergence? - slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs); - if (converged)break; + if(slcconvergence(loads->vsealevelloads,oldsealevelloads,eps_rel,eps_abs))break; //early return? if(iterations>=max_nonlinear_iterations)break; - else iterations++; + iterations++; } @@ -391,5 +360,5 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeDeformationConvolution(allsealevelloads, allsubelementsealevelloads,allloads, allsubelementiceloads, allsubelementhydroloads, allsubelementbploads,rotationaxismotionvector,slgeom); + element->SealevelchangeDeformationConvolution(loads, rotationaxismotionvector,slgeom); } @@ -414,5 +383,5 @@ } /*}}}*/ -void dynstr_core(FemModel* femmodel){ /*{{{*/ +void dynstr_core(FemModel* femmodel){ /*{{{*/ /*variables:*/ @@ -480,5 +449,5 @@ } /*}}}*/ -void coupleroutput_core(FemModel* femmodel){ /*{{{*/ +void coupleroutput_core(FemModel* femmodel){ /*{{{*/ /*parameters:*/ @@ -500,5 +469,5 @@ } }; /*}}}*/ -void ivins_deformation_core(FemModel* femmodel){ /*{{{*/ +void ivins_deformation_core(FemModel* femmodel){ /*{{{*/ int gsize; @@ -551,7 +520,5 @@ } /*}}}*/ - -//Geometry: -void sealevelchange_initialgeometry(FemModel* femmodel) { /*{{{*/ +void sealevelchange_initialgeometry(FemModel* femmodel) { /*{{{*/ /*Geometry core where we compute geometrical kernels and weights:*/ @@ -668,5 +635,167 @@ }/*}}}*/ -/*Support routines:*/ +/*subroutines:*/ +bool slcconvergence(Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ + + bool converged=true; + IssmDouble ndS,nS; + Vector *dRSLg = NULL; + + //compute norm(du) and norm(u) if requested + dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0); + ndS=dRSLg->Norm(NORM_TWO); + + if (xIsNan(ndS)){ + _error_("convergence criterion is NaN (RSL_old=" << RSLg_old->Norm(NORM_TWO) << " RSL=" << RSLg->Norm(NORM_TWO) << ")"); + } + + if(!xIsNan(eps_rel)){ + nS=RSLg_old->Norm(NORM_TWO); + if (xIsNan(nS)) _error_("convergence criterion is NaN! (check the initial RSL)"); + } + + //clean up + delete dRSLg; + + //print + if(!xIsNan(eps_rel)){ + if((ndS/nS) " << eps_rel*100 << " %\n"); + converged=false; + } + } + if(!xIsNan(eps_abs)){ + if(ndS " << eps_abs << " \n"); + converged=false; + } + } + + /*assign output*/ + return converged; + +} /*}}}*/ +IssmDouble SealevelloadsOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* suboceanareas, IssmDouble oceanarea){ /*{{{*/ + + IssmDouble sealevelloadsaverage; + IssmDouble subsealevelloadsaverage; + + Vector* vsealevelloadsvolume=loads->vsealevelloads->Duplicate(); + Vector* vsubsealevelloadsvolume=loads->vsubsealevelloads->Duplicate(); + + vsealevelloadsvolume->PointwiseMult(loads->vsealevelloads,oceanareas); + vsubsealevelloadsvolume->PointwiseMult(loads->vsubsealevelloads,suboceanareas); + + vsealevelloadsvolume->Sum(&sealevelloadsaverage); + vsubsealevelloadsvolume->Sum(&subsealevelloadsaverage); + delete vsealevelloadsvolume; + delete vsubsealevelloadsvolume; + + //return (sealevelloadsaverage+subsealevelloadsaverage)/oceanarea; + return (sealevelloadsaverage)/oceanarea; +} /*}}}*/ +void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom){ /*{{{*/ + + IssmDouble moi_list[3]={0,0,0}; + IssmDouble moi_list_sub[3]={0,0,0}; + IssmDouble moi_list_cpu[3]={0,0,0}; + IssmDouble* tide_love_h = NULL; + IssmDouble* tide_love_k = NULL; + IssmDouble* load_love_k = NULL; + IssmDouble tide_love_k2secular; + IssmDouble moi_e, moi_p; + IssmDouble m1, m2, m3; + bool rotation=false; + + /*early return?:*/ + femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); + if(!rotation)return; + + /*retrieve parameters: */ + femmodel->parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum); + femmodel->parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum); + femmodel->parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum); + femmodel->parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum); + femmodel->parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum); + femmodel->parameters->FindParam(&moi_p,RotationalPolarMoiEnum); + + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + + element->SealevelchangeMomentOfInertiaCentroid(&moi_list[0],loads,slgeom); + + element->SealevelchangeMomentOfInertiaSubElement(&moi_list_sub[0],loads, slgeom); + + moi_list_cpu[0] += moi_list[0]+moi_list_sub[0]; + moi_list_cpu[1] += moi_list[1]+moi_list_sub[1]; + moi_list_cpu[2] += moi_list[2]+moi_list_sub[2]; + } + + + ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*compute perturbation terms for angular velocity vector: */ + m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0]; + m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1]; + m3 = -(1+load_love_k[2])/moi_p * moi_list[2]; // term associated with fluid number (3-order-of-magnitude smaller) is negelected + + /*Assign output pointers:*/ + m[0]=m1; + m[1]=m2; + m[2]=m3; +} /*}}}*/ +void ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/ + + /*Shift sealevel loads by ocean average, only on ocean! :*/ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeShift(loads, offset,slgeom); + } + loads->AssembleSealevelLoads(); + +} /*}}}*/ +IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel){ /*{{{*/ + + int* indices=xNew(nel); + for(int i=0;i* vloadcopy=new Vector(nel); + IssmDouble* loadcopy=xNew(nel); + + vloadcopy->SetValues(nel,indices,load,INS_VAL); + vloadcopy->Assemble(); + + + if(subload){ + for (int i=0;ielements->Size();i++){ + if (slgeom->issubelement[loadtype][i]){ + int se= slgeom->subelementmapping[loadtype][i]; + IssmDouble subloadi=subload[se]; + Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + vloadcopy->SetValue(element->Sid(),subloadi,ADD_VAL); + } + } + } + vloadcopy->Assemble(); + loadcopy=vloadcopy->ToMPISerial(); + + return loadcopy; + +} /*}}}*/ + +/*Coupling routines:*/ void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/ @@ -914,164 +1043,2 @@ } /*}}}*/ -void slcconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/ - - bool converged=true; - IssmDouble ndS,nS; - Vector *dRSLg = NULL; - - //compute norm(du) and norm(u) if requested - dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0); - ndS=dRSLg->Norm(NORM_TWO); - - if (xIsNan(ndS)){ - _error_("convergence criterion is NaN (RSL_old=" << RSLg_old->Norm(NORM_TWO) << " RSL=" << RSLg->Norm(NORM_TWO) << ")"); - } - - if(!xIsNan(eps_rel)){ - nS=RSLg_old->Norm(NORM_TWO); - if (xIsNan(nS)) _error_("convergence criterion is NaN! (check the initial RSL)"); - } - - //clean up - delete dRSLg; - - //print - if(!xIsNan(eps_rel)){ - if((ndS/nS) " << eps_rel*100 << " %\n"); - converged=false; - } - } - if(!xIsNan(eps_abs)){ - if(ndS " << eps_abs << " \n"); - converged=false; - } - } - - /*assign output*/ - *pconverged=converged; - -} /*}}}*/ -IssmDouble SealevelloadsOceanAverage(Vector* sealevelloads,Vector* subelementsealevelloads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble oceanarea){ /*{{{*/ - - IssmDouble sealevelloadsaverage; - IssmDouble subelementsealevelloadsaverage; - - Vector* sealevelloadsvolume=sealevelloads->Duplicate(); - Vector* subelementsealevelloadsvolume=subelementsealevelloads->Duplicate(); - - sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas); - subelementsealevelloadsvolume->PointwiseMult(subelementsealevelloads,subelementoceanareas); - - sealevelloadsvolume->Sum(&sealevelloadsaverage); - subelementsealevelloadsvolume->Sum(&subelementsealevelloadsaverage); - delete sealevelloadsvolume; - delete subelementsealevelloadsvolume; - - //return (sealevelloadsaverage+subelementsealevelloadsaverage)/oceanarea; - return (sealevelloadsaverage)/oceanarea; -} /*}}}*/ -void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){ /*{{{*/ - - IssmDouble moi_list[3]={0,0,0}; - IssmDouble moi_list_sub[3]={0,0,0}; - IssmDouble moi_list_cpu[3]={0,0,0}; - IssmDouble* tide_love_h = NULL; - IssmDouble* tide_love_k = NULL; - IssmDouble* load_love_k = NULL; - IssmDouble tide_love_k2secular; - IssmDouble moi_e, moi_p; - IssmDouble m1, m2, m3; - bool rotation=false; - - /*early return?:*/ - femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); - if(!rotation)return; - - /*retrieve parameters: */ - femmodel->parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum); - femmodel->parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum); - femmodel->parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum); - femmodel->parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum); - femmodel->parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum); - femmodel->parameters->FindParam(&moi_p,RotationalPolarMoiEnum); - - for(Object* & object : femmodel->elements->objects){ - Element* element = xDynamicCast(object); - - element->SealevelchangeMomentOfInertiaCentroid(&moi_list[0],loads,sealevelloads,slgeom); - - element->SealevelchangeMomentOfInertiaSubElement(&moi_list_sub[0],subelementiceloads, subelementhydroloads, subelementbploads, subelementsealevelloads,slgeom); - - moi_list_cpu[0] += moi_list[0]+moi_list_sub[0]; - moi_list_cpu[1] += moi_list[1]+moi_list_sub[1]; - moi_list_cpu[2] += moi_list[2]+moi_list_sub[2]; - } - - - ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*compute perturbation terms for angular velocity vector: */ - m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0]; - m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1]; - m3 = -(1+load_love_k[2])/moi_p * moi_list[2]; // term associated with fluid number (3-order-of-magnitude smaller) is negelected - - /*Assign output pointers:*/ - m[0]=m1; - m[1]=m2; - m[2]=m3; -} /*}}}*/ -void ConserveOceanMass(FemModel* femmodel,Vector* sealevelloads, Vector* subelementsealevelloads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/ - - /*Shift sealevel loads by ocean average, only on ocean! :*/ - for(Object* & object : femmodel->elements->objects){ - Element* element = xDynamicCast(object); - element->SealevelchangeShift(sealevelloads,subelementsealevelloads,offset,slgeom); - } - sealevelloads->Assemble(); - subelementsealevelloads->Assemble(); - -} /*}}}*/ -IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel){ /*{{{*/ - - int* indices=xNew(nel); - for(int i=0;i* vloadcopy=new Vector(nel); - IssmDouble* loadcopy=xNew(nel); - - vloadcopy->SetValues(nel,indices,load,INS_VAL); - vloadcopy->Assemble(); - - - if(subload){ - for (int i=0;ielements->Size();i++){ - if (slgeom->issubelement[loadtype][i]){ - int se= slgeom->subelementmapping[loadtype][i]; - IssmDouble subloadi=subload[se]; - Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); - vloadcopy->SetValue(element->Sid(),subloadi,ADD_VAL); - } - } - } - vloadcopy->Assemble(); - loadcopy=vloadcopy->ToMPISerial(); - - return loadcopy; - -} /*}}}*/ -