Changeset 27017

Timestamp:

05/20/22 09:52:09 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moving min thickness as postprocessing. Will move other calving laws as well if successful

Location:

issm/trunk-jpl/src/c

Files:

• analyses/LevelsetAnalysis.cpp (modified) (2 diffs)
• analyses/LevelsetAnalysis.h (modified) (1 diff)
• cores/movingfront_core.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

@@ -641 +641 @@
         }
 }/*}}}*/
+void           LevelsetAnalysis::PostProcess(FemModel* femmodel){/*{{{*/
+
+        /*This function is only used by "discrete calving laws" for which we change
+         * the value of the levelset after the advection step (level set equation
+         * solve) based on the law*/
+
+        /*Intermediaries*/
+        int  calvinglaw;
+        IssmDouble newlevelset[6];
+        femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
+
+        if(calvinglaw==CalvingMinthicknessEnum){
+
+                IssmDouble mig_max = femmodel->parameters->FindParam(MigrationMaxEnum);
+                IssmDouble dt      = femmodel->parameters->FindParam(TimesteppingTimeStepEnum);
+
+                /*Get current distance to terminus*/
+                InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
+                femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
+
+                /*Intermediaries*/
+                IssmDouble thickness,bed,sealevel,distance,levelset;
+                IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum);
+
+                /*Loop over all elements of this partition*/
+                for(Object* & object : femmodel->elements->objects){
+                        Element* element  = xDynamicCast<Element*>(object);
+
+                        /*no need to postprocess an ice free element*/
+                        if(!element->IsIceInElement()) continue;
+
+                        int      numnodes     = element->GetNumberOfNodes(); _assert_(numnodes<7);
+                        Gauss*   gauss        = element->NewGauss();
+                        Input *H_input        = element->GetInput(ThicknessEnum);              _assert_(H_input);
+                        Input *b_input        = element->GetInput(BedEnum);                    _assert_(b_input);
+                        Input *sl_input       = element->GetInput(SealevelEnum);               _assert_(sl_input);
+                        Input *dis_input      = element->GetInput(DistanceToCalvingfrontEnum); _assert_(dis_input);
+                        Input *levelset_input = element->GetInput(MaskIceLevelsetEnum);        _assert_(levelset_input);
+
+                        /*Potentially constrain nodes of this element*/
+                        for(int in=0;in<numnodes;in++){
+                                gauss->GaussNode(element->GetElementType(),in);
+
+                                levelset_input->GetInputValue(&levelset,gauss);
+                                H_input->GetInputValue(&thickness,gauss);
+                                b_input->GetInputValue(&bed,gauss);
+                                sl_input->GetInputValue(&sealevel,gauss);
+                                dis_input->GetInputValue(&distance,gauss);
+
+                                if(thickness<min_thickness && bed<sealevel && fabs(distance)<mig_max*dt && levelset<0){
+                                        newlevelset[in] = +400.; //Arbitrary > 0 number
+                                }
+                                else{
+                                        newlevelset[in] = levelset;
+                                }
+                        }
+                        element->AddInput(MaskIceLevelsetEnum,&newlevelset[0],element->GetElementType());
+                        delete gauss;
+                }
+        }
+}/*}}}*/
 void           LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 
@@ -653 +714 @@
    femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
 
-        if(calvinglaw==CalvingMinthicknessEnum || calvinglaw==CalvingVonmisesEnum){
+        if(calvinglaw==CalvingVonmisesEnum){
 
                 /*Intermediaries*/

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h

@@ -32 +32 @@
                 void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
                 void           UpdateConstraints(FemModel* femmodel);
+
+                /*Specific methods to LevelsetAnalysis*/
+                void  PostProcess(FemModel* femmodel);
 };
 #endif

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

@@ -114 +114 @@
 
         /* solve level set equation */
-        analysis = new LevelsetAnalysis();
-        analysis->Core(femmodel);
-        delete analysis;
+        LevelsetAnalysis lsanalysis;
+        lsanalysis.Core(femmodel);
+        lsanalysis.PostProcess(femmodel);
 
         /*Kill ice berg to avoid free body motion*/