Changeset 16699

Timestamp:

11/11/13 09:07:42 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added case in thermalanalysis

Location:

issm/trunk-jpl/src/c

Files:

• analyses/ThermalAnalysis.cpp (modified) (1 diff)
• classes/Elements/Element.h (modified) (2 diffs)
• classes/Elements/Penta.cpp (modified) (2 diffs)
• classes/Elements/Penta.h (modified) (2 diffs)
• classes/Elements/Seg.h (modified) (2 diffs)
• classes/Elements/Tria.cpp (modified) (2 diffs)
• classes/Elements/Tria.h (modified) (2 diffs)

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

@@ -116 +116 @@
 }/*}}}*/
 void ThermalAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-        _error_("not implemented yet");
+
+        bool        converged;
+        int         i,rheology_law;
+        IssmDouble  B_average,s_average,T_average=0.;
+        int        *doflist   = NULL;
+        IssmDouble *xyz_list  = NULL;
+        bool        hack      = false;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Fetch dof list and allocate solution vector*/
+        element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+        IssmDouble* values    = xNew<IssmDouble>(numnodes);
+
+        /*Use the dof list to index into the solution vector: */
+        for(i=0;i<numnodes;i++){
+                values[i]=solution[doflist[i]];
+
+                /*Check solution*/
+                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
+                //if(values[i]<0)      _printf_("temperature < 0°K found in solution vector\n");
+                //if(values[i]>275)    _printf_("temperature > 275°K found in solution vector (Paterson's rheology associated is negative)\n");
+        }
+
+        /*Force temperature between [Tpmp-50 Tpmp] to disable penalties*/
+        if(hack){
+                IssmDouble* pressure = xNew<IssmDouble>(numnodes);
+                element->GetInputListOnVertices(pressure,PressureEnum);
+                for(i=0;i<numnodes;i++){
+                        if(values[i]>element->TMeltingPoint(pressure[i]))     values[i]=element->TMeltingPoint(pressure[i]);
+                        if(values[i]<element->TMeltingPoint(pressure[i])-50.) values[i]=element->TMeltingPoint(pressure[i])-50.;
+                }
+                xDelete<IssmDouble>(pressure);
+        }
+
+        /*Get all inputs and parameters*/
+
+        element->GetInputValue(&converged,ConvergedEnum);
+        if(converged){
+                element->AddInput(TemperatureEnum,values,P1Enum);
+
+                /*Update Rheology only if converged (we must make sure that the temperature is below melting point
+                 * otherwise the rheology could be negative*/
+                element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+                switch(rheology_law){
+                        case NoneEnum:
+                                /*Do nothing: B is not temperature dependent*/
+                                break;
+                        case PatersonEnum:
+                                for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
+                                B_average=Paterson(T_average);
+                                element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+                                break;
+                        case ArrheniusEnum:{
+                                Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
+                                surface_input->GetInputAverage(&s_average);
+                                element->GetVerticesCoordinates(&xyz_list);
+                                for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
+                                //B_average=Arrhenius(T_average,
+                                                        //s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0),
+                                                        //element->GetMaticeParameter(MaterialsRheologyNEnum));
+                                element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+                                break;
+                                }
+                        default:
+                                _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
+
+                }
+        }
+        else{
+                element->AddInput(TemperaturePicardEnum,values,P1Enum);
+        }
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(values);
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<int>(doflist);
 }/*}}}*/

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

@@ -40 +40 @@
                 virtual void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
                 virtual void   AddInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
+                virtual void   AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
                 virtual void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>*  Kfs)=0;
                 virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
@@ -95 +96 @@
                 virtual void   Delta18oParameterization(void)=0;
                 virtual void   SmbGradients()=0;
+                virtual IssmDouble TMeltingPoint(IssmDouble pressure)=0;
                 virtual void   ResetCoordinateSystem()=0;
                 virtual int    VelocityInterpolation()=0;

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

@@ -124 +124 @@
 void  Penta::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){
 
-        /*Call inputs method*/
+        _assert_(this->inputs);
         this->inputs->AddInput(new PentaInput(input_enum,values,interpolation_enum));
 }
 /*}}}*/
+/*FUNCTION Penta::AddMaterialInput{{{*/
+void  Penta::AddMaterialInput(int input_enum,IssmDouble* values, int interpolation_enum){
+
+        _assert_(this->material);
+        this->material->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
+}
+/*}}}*/
+
 /*FUNCTION Penta::BedNormal {{{*/
 void Penta::BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]){
@@ -3151 +3159 @@
           /*Update inputs*/
           this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&smb[0],P1Enum));
+}
+/*}}}*/
+/*FUNCTION Penta::TMeltingPoint{{{*/
+IssmDouble Penta::TMeltingPoint(IssmDouble pressure){
+
+        _assert_(matpar);
+        return this->matpar->TMeltingPoint(pressure);
+
 }
 /*}}}*/

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

@@ -112 +112 @@
                 void   ResetCoordinateSystem(void);
                 void   SmbGradients();
+                IssmDouble  TMeltingPoint(IssmDouble pressure);
                 IssmDouble SurfaceArea(void);
                 void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
@@ -182 +183 @@
                 /*Penta specific routines:{{{*/
                 void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
+                void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);
                 void             BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]);
                 ElementMatrix* CreateBasalMassMatrix(void);

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

@@ -67 +67 @@
                 /*Element virtual functions definitions: {{{*/
                 void        AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
+                void        AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
                 void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
                 void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
@@ -137 +138 @@
                 void        ResetCoordinateSystem(void){_error_("not implemented yet");};
                 void          SmbGradients(){_error_("not implemented yet");};
+                IssmDouble  TMeltingPoint(IssmDouble pressure){_error_("not implemented yet");};
                 IssmDouble  SurfaceArea(void){_error_("not implemented yet");};
                 void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement){_error_("not implemented yet");};

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

@@ -170 +170 @@
 
         /*Call inputs method*/
+        _assert_(this->inputs);
         this->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
+}
+/*}}}*/
+/*FUNCTION Tria::AddMaterialInput{{{*/
+void  Tria::AddMaterialInput(int input_enum,IssmDouble* values, int interpolation_enum){
+
+        _assert_(this->material);
+        this->material->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
 }
 /*}}}*/
@@ -2641 +2649 @@
           /*Update inputs*/
           this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&smb[0],P1Enum));
+}
+/*}}}*/
+/*FUNCTION Tria::TMeltingPoint{{{*/
+IssmDouble Tria::TMeltingPoint(IssmDouble pressure){
+
+        _assert_(matpar);
+        return this->matpar->TMeltingPoint(pressure);
+
 }
 /*}}}*/

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

@@ -118 +118 @@
                 void        ResetCoordinateSystem(void);
                 void          SmbGradients();
+                IssmDouble  TMeltingPoint(IssmDouble pressure);
                 int         VelocityInterpolation();
                 int         PressureInterpolation();
@@ -193 +194 @@
                 /*Tria specific routines:{{{*/
                 void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
+                void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);
                 ElementMatrix* CreateKMatrix(void);
                 ElementMatrix* CreateKMatrixBalancethickness(void);