Changeset 26255

Timestamp:

05/10/21 12:30:20 (4 years ago)

Author:

Cheng Gong

Message:

CHG: add marshall and a test for the levelset misfit output definition

Location:

issm/trunk-jpl

Files:

• src/c/analyses/LevelsetAnalysis.cpp (modified) (10 diffs)
• src/c/classes/Cflevelsetmisfit.cpp (modified) (1 diff)
• src/c/datastructures/DataSet.cpp (modified) (1 diff)
• test/Archives/Archive3203.arch (added)
• test/NightlyRun/test3203.m (added)

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

@@ -10 +10 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+#include "../classes/Inputs/DatasetInput.h"
+
 
 void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
@@ -49 +51 @@
 
         iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
+//      iomodel->FetchDataToInput(inputs,elements,"md.mask.ref_levelset",MaskIceRefLevelsetEnum);
         iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
         iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
@@ -67 +70 @@
                         iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
                         break;
+                case CalvingTestEnum:
+                        iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum);
+                        iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum);
+                        iomodel->FetchDataToInput(inputs,elements,"md.calving.Qrates",CalvingNonlinearQRateEnum);
+                        iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
+                        break;
                 case CalvingMinthicknessEnum:
                         iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
@@ -103 +112 @@
 void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 
+        IssmDouble* coefficients=NULL;
+        int nco;
+
         parameters->AddObject(iomodel->CopyConstantObject("md.levelset.stabilization",LevelsetStabilizationEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.levelset.reinit_frequency",LevelsetReinitFrequencyEnum));
@@ -117 +129 @@
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum));
                         break;
+                case CalvingTestEnum:
+                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.dt_order",CalvingDtimeorderEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.nonlinear_law",CalvingNonlinearLawEnum));
+                        
+                        /* Add the coefficients for the new nonlinear law */
+                        iomodel->FetchData(&coefficients, &nco, NULL, "md.calving.coefficients");
+              parameters->AddObject(new DoubleMatParam(CalvingNonlinearCoeffEnum, coefficients, nco, 1));
+                        xDelete<IssmDouble>(coefficients);
                 case CalvingMinthicknessEnum:
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum));
@@ -215 +236 @@
         IssmDouble*    dbasis   = xNew<IssmDouble>(2*numnodes);
         IssmDouble*    Bprime = NULL;
+        IssmDouble              tau;
         if(stabilization==2){
                 Bprime   = xNew<IssmDouble>(dim*numnodes);
@@ -340 +362 @@
 
         /*Intermediaries */
-        int domaintype;
+        int domaintype, stabilization;
         IssmDouble  Jdet,dt;
         IssmDouble  lsf;
         IssmDouble* xyz_list = NULL;
+   IssmDouble  vx,vy,vel,h,tau;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -351 +374 @@
         ElementVector* pe = basalelement->NewElementVector();
         basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
+   basalelement->FindParam(&stabilization,LevelsetStabilizationEnum);
 
         if(dt!=0.){
@@ -359 +383 @@
                 basalelement->GetVerticesCoordinates(&xyz_list);
                 Input* levelset_input     = basalelement->GetInput(MaskIceLevelsetEnum);                    _assert_(levelset_input);
+                
+        //      h=element->CharacteristicLength();
 
                 /* Start  looping on the number of gaussian points: */
@@ -473 +499 @@
         }
 
+        if(calvinglaw==CalvingTestEnum){
+
+                /*Get minimum thickness threshold*/
+                femmodel->parameters->FindParam(&min_thickness,CalvingMinthicknessEnum);
+
+                /*Loop over all elements of this partition*/
+                for(Object* & object : femmodel->elements->objects){
+                        Element* element  = xDynamicCast<Element*>(object);
+
+                        int      numnodes = element->GetNumberOfNodes();
+                        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);
+
+                        /*Potentially constrain nodes of this element*/
+                        for(int in=0;in<numnodes;in++){
+                                gauss->GaussNode(element->GetElementType(),in);
+                                Node* node=element->GetNode(in);
+                                if(!node->IsActive()) continue;
+
+                                H_input->GetInputValue(&thickness,gauss);
+                                b_input->GetInputValue(&bed,gauss);
+                                sl_input->GetInputValue(&sealevel,gauss);
+                                if(thickness<min_thickness && bed<sealevel){
+                                        node->ApplyConstraint(0,+1.);
+                                }
+                                else {
+                                        /* no ice, set no spc */
+                                        node->DofInFSet(0);
+                                }
+                        }
+                        delete gauss;
+                }
+        }
         if(calvinglaw==CalvingHabEnum){
 

issm/trunk-jpl/src/c/classes/Cflevelsetmisfit.cpp ¶

@@ -89 +89 @@
 
         marshallhandle->call(this->definitionenum);
-        marshallhandle->call(this->local);
         marshallhandle->call(this->model_enum);
         marshallhandle->call(this->name);
         marshallhandle->call(this->observation_enum);
-        marshallhandle->call(this->timeinterpolation);
         marshallhandle->call(this->weights_enum);
         marshallhandle->call(this->datatime);

issm/trunk-jpl/src/c/datastructures/DataSet.cpp ¶

@@ -251 +251 @@
                                 this->AddObject(res);
                         }
+                        else if(obj_enum==CflevelsetmisfitEnum){
+                                Cflevelsetmisfit* Cflevelset=new Cflevelsetmisfit();
+                                Cflevelset->Marshall(marshallhandle);
+                                this->AddObject(Cflevelset);
+                        }
                         else _error_("could not recognize enum type: " << obj_enum << ": " << EnumToStringx(obj_enum) ); 
                 }