Changeset 26152

Timestamp:

03/26/21 10:48:14 (4 years ago)

Author:

Eric.Larour

Message:

CHG: SealevelMasks is not working on nel, but rather on local number of elements!
sealevelchange_core and Tria.cpp: need a routine to shift only on ocean for
mass conservation.

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Penta.h (modified) (1 diff)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/Elements/Tetra.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (1 diff)
• classes/Elements/Tria.h (modified) (1 diff)
• classes/SealevelMasks.h (modified) (1 diff)
• cores/sealevelchange_core.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -398 +398 @@
                 virtual void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector,  SealevelMasks* masks)=0;
                 virtual void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0;
+                virtual void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks)=0;
                 #endif
 

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -235 +235 @@
                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
                 void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
+                void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
                 #endif
 

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -190 +190 @@
                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
                 void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
+                void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
 
 #endif

issm/trunk-jpl/src/c/classes/Elements/Tetra.h

@@ -196 +196 @@
                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
                 void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
+                void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
 
 #endif

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -6965 +6965 @@
         return;
 }/*}}}*/
+void       Tria::SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){ /*{{{*/
+        
+        if(masks->isoceanin[this->lid]){
+                loads->SetValue(this->sid,offset,ADD_VAL);
+        }
+
+} /*}}}*/
 #endif
 

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -180 +180 @@
                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks);
                 void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks);
+                void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks);
                 #endif
                 /*}}}*/

issm/trunk-jpl/src/c/classes/SealevelMasks.h

@@ -18 +18 @@
                 
                 /*SealevelMasks constructors, destructors :*/
-                SealevelMasks(int nel){
+                SealevelMasks(int localnel){
                         /*allocate fields:*/
-                        this->isiceonly=xNew<bool>(nel);
-                        this->isfullyfloating=xNew<bool>(nel);
-                        this->notfullygrounded=xNew<bool>(nel);
-                        this->isoceanin=xNew<bool>(nel);
+                        this->isiceonly=xNew<bool>(localnel);
+                        this->isfullyfloating=xNew<bool>(localnel);
+                        this->notfullygrounded=xNew<bool>(localnel);
+                        this->isoceanin=xNew<bool>(localnel);
                 };
                 ~SealevelMasks(){

issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp

@@ -12 +12 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+
 /*support routines local definitions:{{{*/
 void TransferForcing(FemModel* femmodel,int forcingenum);
@@ -19 +20 @@
 SealevelMasks* sealevel_masks(FemModel* femmodel);
 void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealeveloads);
+void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, IssmDouble offset, SealevelMasks* masks);
 /*}}}*/
 
@@ -503 +505 @@
         if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
 
-        //substract ocean average and barystatic contributionfrom sea level loads: 
-        sealevelloads->Shift(barycontrib->Total()/oceanarea - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
-        
+        //conserve ocean mass: 
+        oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
+        ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks);
+
         //broadcast sea level loads 
         allsealevelloads=sealevelloads->ToMPISerial();
+
+        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,DummyEnum,allsealevelloads,nel,1,1,1));
 
         //compute rotation axis motion:
@@ -523 +528 @@
                 }
                 sealevelloads->Assemble();
-        
-                //substract ocean average and barystatic contribution
+                        
+                //Conserve ocean mass: 
                 oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
-                sealevelloads->Shift(barycontrib->Total()/oceanarea- oceanaverage);
-                
+                ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks);
+
                 //broadcast sea level loads 
                 allsealevelloads=sealevelloads->ToMPISerial();
@@ -1318 +1323 @@
         m[2]=m3;
 } /*}}}*/
+void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, IssmDouble offset, SealevelMasks* masks){ /*{{{*/
+
+        /*Shift sealevel loads by ocean average, only on ocean! :*/
+        for(Object* & object : femmodel->elements->objects){
+                Element* element = xDynamicCast<Element*>(object);
+                element->SealevelchangeShift(sealevelloads,offset,masks);
+        }
+        sealevelloads->Assemble();
+
+} /*}}}*/