Changeset 24556

Timestamp:

02/07/20 16:27:39 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on Josh's model

Location:

issm/trunk-jpl/src/c

Files:

• analyses/EnthalpyAnalysis.cpp (modified) (4 diffs)
• classes/Elements/Element.cpp (modified) (5 diffs)
• classes/Inputs2/PentaInput2.cpp (modified) (2 diffs)

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

@@ -507 +507 @@
                 element->AddInput2(WatercolumnEnum,watercolumns,finite_element);
         }
-        element->AddInput2(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,finite_element);
+        element->AddInput2(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,P1DGEnum);
 
         /*Clean up and return*/
@@ -825 +825 @@
         Input2* pressure_input=element->GetInput2(PressureEnum); _assert_(pressure_input);
         Input2* enthalpy_input=NULL;
-        if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput2(EnthalpyEnum); _assert_(enthalpy_input);}
+        if(dt>0.){
+                enthalpy_input = element->GetInput2(EnthalpyEnum); _assert_(enthalpy_input);
+        }
 
         /* Start  looping on the number of gaussian points: */
@@ -864 +866 @@
 
                 /* Build transient now */
-                if(reCast<bool,IssmDouble>(dt)){
+                if(dt>0.){
                         enthalpy_input->GetInputValue(&enthalpy, gauss);
                         scalar_transient=enthalpy*Jdet*gauss->weight;
@@ -891 +893 @@
                         element->ElementSizes(&hx,&hy,&hz);
                         kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
-                        vx_input->GetInputValue(&u,gauss);
-                        vy_input->GetInputValue(&v,gauss);
-                        vz_input->GetInputValue(&w,gauss);
-                        element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa/rho_ice);
-                        tau_parameter_hor=tau_parameter_anisotropic[0];
-                        tau_parameter_ver=tau_parameter_anisotropic[1];
-                        
-                        for(i=0;i<numnodes;i++) pe->values[i]+=scalar_def*(tau_parameter_hor*u*dbasis[0*numnodes+i]+tau_parameter_hor*v*dbasis[1*numnodes+i]+tau_parameter_ver*w*dbasis[2*numnodes+i]);
+                        vx_input->GetInputValue(&u,gauss);
+                        vy_input->GetInputValue(&v,gauss);
+                        vz_input->GetInputValue(&w,gauss);
+                        element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa/rho_ice);
+                        tau_parameter_hor=tau_parameter_anisotropic[0];
+                        tau_parameter_ver=tau_parameter_anisotropic[1];
+
+                        for(i=0;i<numnodes;i++) pe->values[i]+=scalar_def*(tau_parameter_hor*u*dbasis[0*numnodes+i]+tau_parameter_hor*v*dbasis[1*numnodes+i]+tau_parameter_ver*w*dbasis[2*numnodes+i]);
                 }
         }

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

@@ -18 +18 @@
 #include "../Inputs2/TransientInput2.h"
 #include "../Inputs2/ElementInput2.h"
+#include "../Inputs2/PentaInput2.h"
 #include "../Inputs2/DatasetInput2.h"
 #include "../Inputs2/ArrayInput2.h"
 /*}}}*/
+#define MAXVERTICES 6 /*Maximum number of vertices per element, currently Penta, to avoid dynamic mem allocation*/
 
 #ifdef _HAVE_SEMIC_
@@ -2804 +2806 @@
         const int NUM_VERTICES_MONTHS_PER_YEAR = NUM_VERTICES * 12;
 
-        int             i;
+        int             i,vertexlids[MAXVERTICES];
         IssmDouble* agd=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
         IssmDouble* melt=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
@@ -2822 +2824 @@
         IssmDouble mavg=1./12.; //factor for monthly average
 
+        /*Get vertex Lids for later*/
+        this->GetVerticesLidList(&vertexlids[0]);
+
         /*Get material parameters :*/
         rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
@@ -2915 +2920 @@
                                  * the temperatures as they are for the base of the penta and
                                  * yse yearlytemperatures for the top*/
-                                GetInputListOnVertices(&s[0],TemperatureEnum);
-                                yearlytemperatures[0] = s[0];
-                                yearlytemperatures[1] = s[1];
-                                yearlytemperatures[2] = s[2];
+                                PentaInput2* temp_input = xDynamicCast<PentaInput2*>(this->GetInput2(TemperatureEnum)); _assert_(temp_input);
+                                switch(temp_input->GetInputInterpolationType()){
+                                        case P1Enum:
+                                                temp_input->element_values[3] = yearlytemperatures[3];
+                                                temp_input->element_values[4] = yearlytemperatures[4];
+                                                temp_input->element_values[5] = yearlytemperatures[5];
+                                                temp_input->SetInput(P1Enum,NUM_VERTICES,&vertexlids[0],temp_input->element_values);
+                                                break;
+                                        default:
+                                                _error_("Interpolation "<<EnumToStringx(temp_input->GetInputInterpolationType())<<" not supported yet");
+                                }
                                 this->AddInput2(TemperatureEnum,&yearlytemperatures[0],P1Enum);
 
@@ -3247 +3259 @@
 } /*}}}*/
 void       Element::ResultToVector(Vector<IssmDouble>* vector,int output_enum){/*{{{*/
-
-        const int MAXVERTICES = 6;
 
         IssmDouble values[MAXVERTICES];

issm/trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp

@@ -85 +85 @@
 /*}}}*/
 void PentaInput2::DeepEcho(void){/*{{{*/
-        _printf_("PentaInput2 Echo:\n");
-        _printf_("   interpolation:      "<<EnumToStringx(this->interpolation)<<"\n");
-        _printf_("   Size:               "<<M<<"x"<<N<<"\n");
-        _printf_("   isserved:           "<<(isserved?"true":"false") << "\n");
-        _printf_("   isserved_collapsed: "<<isserved_collapsed << "\n");
-        if(isserved){
-                _printf_("   current values:      ");
-                _printf_("[ ");
-                for(int i=0;i<6;i++) _printf_(" "<<this->element_values[i]);
-                _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
-        }
-        printarray(this->values,this->M,this->N);
-        //_printf_(setw(15)<<"   PentaInput2 "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<" "<<(value?"true":"false") << "\n");
-}
-/*}}}*/
-void PentaInput2::Echo(void){/*{{{*/
         _printf_("PentaInput2 Echo:\n");
         _printf_("   interpolation:      "<<EnumToStringx(this->interpolation)<<"\n");
@@ -123 +107 @@
 }
 /*}}}*/
+void PentaInput2::Echo(void){/*{{{*/
+        _printf_(setw(15)<<"   PentaInput "<<setw(25)<<left<<EnumToStringx(-1));
+        if(isserved){
+                if(isserved_collapsed){
+                        _printf_("[ ");
+                        for(int i=0;i<3;i++) _printf_(" "<<this->element_values[i]);
+                        _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
+                }
+                else{
+                        _printf_("[ ");
+                        for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) _printf_(" "<<this->element_values[i]);
+                        _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
+                }
+        }
+}
+/*}}}*/
 int  PentaInput2::Id(void){/*{{{*/
         return -1;