Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26229) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26230) @@ -397,5 +397,6 @@ virtual void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=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 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 SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26229) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26230) @@ -230,5 +230,6 @@ void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_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 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 SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26229) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26230) @@ -185,5 +185,6 @@ void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_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 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 SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26229) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26230) @@ -192,5 +192,6 @@ void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_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 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 SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26229) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26230) @@ -6305,6 +6305,4 @@ #endif /*}}}*/ - - return; @@ -6351,4 +6349,10 @@ slgeom->isoceanin[this->lid]=isocean; //keep track for later. area=areae[this->sid]; + + /*Compute element ids, used to speed up computations in convolution phase:{{{*/ + for(int i=0;ilids[this->vertices[i]->lid]=this->lid; + } + /*}}}*/ /*set barycentre for all elements, to be updated for fractional loads in the next routine: */ @@ -6884,5 +6888,5 @@ } /*}}}*/ -void Tria::SealevelchangeConvolution(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/ +void Tria::SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/ /*sal green function:*/ @@ -6896,4 +6900,5 @@ bool sal = false; bool rotation= false; + bool percpu= false; int size; int nel,nbar; @@ -6907,4 +6912,13 @@ this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); + if(sal){ + this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&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); + + this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true); + } + if(rotation){ Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum); @@ -6912,19 +6926,26 @@ Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum); - for(int i=0;iinputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&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); - - this->SealevelchangeGxL(&SealevelGRD[0],G, Gsub, loads, slgeom, nel); - } + 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]; + } + } + } + return; +} /*}}}*/ +void Tria::SealevelchangeOceanAverage(GrdLoads* loads, Vector* oceanareas, Vector* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){ /*{{{*/ + + /*sal green function:*/ + IssmDouble* G=NULL; + IssmDouble* Gsub[SLGEOM_NUMLOADS]; + IssmDouble oceanaverage=0; + IssmDouble oceanarea=0; + IssmDouble rho_water; + + this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); /*retrieve ocean average and area:*/ for(int i=0;iLoadWeigths[SLGEOM_OCEAN][i][this->lid]; + oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]; } #ifdef _ISSM_DEBUG_ @@ -6940,4 +6961,5 @@ else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL); + /*add ocean area into a global oceanareas vector:*/ if(!loads->sealevelloads){ oceanareas->SetValue(this->sid,oceanarea,INS_VAL); @@ -6947,5 +6969,4 @@ } } - return; @@ -6976,4 +6997,5 @@ bool rotation= false; bool elastic=false; + bool percpu=false; this->parameters->FindParam(&nel,MeshNumberofelementsEnum); @@ -6983,11 +7005,5 @@ this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); - if(rotation){ - Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum); - Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum); - Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum); - - for(int i=0;iSealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel); + this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel,percpu=false); if(elastic){ - this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel); + this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel,percpu=false); if(horiz ){ - this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel); - this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel); + this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel,percpu=false); + this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel,percpu=false); + } + } + } + + if(rotation){ + Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum); + Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum); + Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum); + + 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]; } } @@ -7040,8 +7068,11 @@ } /*}}}*/ -void Tria::SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel) { /*{{{*/ +void Tria::SealevelchangeGxL(IssmDouble* sealevelout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu) { /*{{{*/ + + IssmDouble sealevel[3]={0}; if(loads->sealevelloads){ for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; for (int e=0;esealevelloads[e]+loads->loads[e]); @@ -7062,4 +7093,5 @@ else{ //this is the initial convolution where only loads are provided for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; for (int e=0;eloads[e]); @@ -7072,4 +7104,14 @@ } } + } + + /*store values computed on vertices, but don't repeat the computation if another element already computed it!:*/ + if(percpu){ + for(int i=0;ilids[this->vertices[i]->lid]==this->lid)sealevelout[this->vertices[i]->lid]=sealevel[i]; + } + } + else{ + for(int i=0;i* oceanareas, Vector* subelementoceanareas, IssmDouble* rotationvector,SealevelGeometry* slgeom); + 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 SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom); @@ -242,5 +243,5 @@ void UpdateConstraintsExtrudeFromBase(void); void UpdateConstraintsExtrudeFromTop(void); - void SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel); + void SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool optimize); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp (revision 26229) +++ /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp (revision 26230) @@ -17,5 +17,5 @@ /*Object constructors and destructor*/ -SealevelGeometry::SealevelGeometry(int localnelin){ /*{{{*/ +SealevelGeometry::SealevelGeometry(int localnelin,int localnodsin){ /*{{{*/ localnel=localnelin; for(int i=0;i(localnel); isoceanin=xNew(localnel); + lids=xNew(localnodsin); }; /*}}}*/ @@ -59,4 +60,5 @@ xDelete(longe); xDelete(isoceanin); + xDelete(lids); }; /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/SealevelGeometry.h =================================================================== --- /issm/trunk-jpl/src/c/classes/SealevelGeometry.h (revision 26229) +++ /issm/trunk-jpl/src/c/classes/SealevelGeometry.h (revision 26230) @@ -35,6 +35,7 @@ int nsubel[SLGEOM_NUMLOADS]; int nbar[SLGEOM_NUMLOADS]; + int* lids; - SealevelGeometry(int localnel); + SealevelGeometry(int localnel,int localnods); ~SealevelGeometry(); void InitializeMappingsAndBarycentres(void); Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26229) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26230) @@ -245,4 +245,5 @@ int grdmodel; int computesealevel=0; + IssmDouble* sealevelpercpu=NULL; /*}}}*/ @@ -292,4 +293,5 @@ subelementoceanareas=new Vector(slgeom->nbar[SLGEOM_OCEAN]); oceanareas=new Vector(nel); + sealevelpercpu=xNewZeroInit(femmodel->vertices->Size()); if(VerboseSolution()) _printf0_(" starting GRD convolutions\n"); @@ -322,5 +324,11 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeConvolution(loads, oceanareas, subelementoceanareas, rotationaxismotionvector,slgeom); + element->SealevelchangeConvolution(sealevelpercpu, loads , rotationaxismotionvector,slgeom); + } + + /*retrieve sea level average and ocean area:*/ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeOceanAverage(loads, oceanareas, subelementoceanareas, sealevelpercpu, slgeom); } @@ -535,5 +543,5 @@ femmodel->parameters->FindParam(&grdmodel,GrdModelEnum); nel=femmodel->elements->NumberOfElements(); - + /*early return?:*/ if(grdmodel==IvinsEnum) return; @@ -594,5 +602,5 @@ /*initialize SealevelMasks structure: */ - slgeom=new SealevelGeometry(femmodel->elements->Size()); + slgeom=new SealevelGeometry(femmodel->elements->Size(),femmodel->vertices->Size()); /*Verbose: */