Changeset 24053

Timestamp:

06/28/19 05:47:22 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added mouling relaxation

Location:

issm/trunk-jpl/src/c/classes

Files:

• ExternalResults/GenericExternalResult.h (modified) (1 diff)
• Loads/Moulin.cpp (modified) (3 diffs)
• Loads/Moulin.h (modified) (1 diff)

issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h

@@ -87 +87 @@
                 } /*}}}*/
                 GenericExternalResult(int in_id, int in_enum_type,ResultType in_values, int in_M,int in_N,int in_step,IssmDouble in_time){/*{{{*/
-                        _error_("template GenericExternalResult(int in_id, int in_enum_type,double* in_values, int in_M,int in_N,int in_step,IssmDouble in_time) not implemented for this ResultType\n");
+                        id        = in_id;
+                        step      = in_step;
+                        time      = in_time;
+                        M         = in_M;
+                        N         = in_N;
+
+                        /*Copy result in values*/
+                        if(M*N){
+                                value=xNew<IssmDouble>(M*N);
+                                xMemCpy<IssmDouble>(value,in_values,M*N);
+                        }
+                        else value=NULL;
+
+                        /*Convert enum to name*/
+                        EnumToStringx(&this->result_name,in_enum_type);
                 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp

@@ -147 +147 @@
 void  Moulin::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
 
-        /*No loads applied, do nothing: */
-        return;
+        /*recover some parameters*/
+        ElementMatrix* Ke=NULL;
+        int analysis_type;
+        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+        switch(analysis_type){
+                case HydrologyGlaDSAnalysisEnum:
+                        Ke = this->CreateKMatrixHydrologyGlaDS();
+                        break;
+                case HydrologyShaktiAnalysisEnum:
+                        /*do nothing: */
+                        return;
+                case HydrologyDCInefficientAnalysisEnum:
+                        /*do nothing: */
+                        return;
+                case HydrologyDCEfficientAnalysisEnum:
+                        /*do nothing: */
+                        return;
+                default:
+                        _error_("Don't know why we should be here");
+
+        }
+        /*Add to global matrix*/
+        if(Ke){
+                Ke->AddToGlobal(Kff,Kfs);
+                delete Ke;
+        }
 
 }
@@ -313 +338 @@
 /*}}}*/
 
-ElementVector* Moulin::CreatePVectorHydrologyGlaDS(void){/*{{{*/
+ElementMatrix* Moulin::CreateKMatrixHydrologyGlaDS(void){/*{{{*/
 
         /*If this node is not the master node (belongs to another partition of the
@@ -319 +344 @@
         if(node->IsClone()) return NULL;
 
-        IssmDouble moulin_load;
-
         /*Initialize Element matrix*/
+        ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters);
+
+        /*Get all inputs and parameters*/
+        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
+        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+        IssmDouble g         = element->FindParam(ConstantsGEnum);
+        IssmDouble Am        = 0.; //For now...
+
+        /*Load vector*/
+        Ke->values[0] = - Am/(rho_water*g)/dt;
+
+        /*Clean up and return*/
+        return Ke;
+}
+/*}}}*/
+ElementVector* Moulin::CreatePVectorHydrologyGlaDS(void){/*{{{*/
+
+        /*If this node is not the master node (belongs to another partition of the
+         * mesh), don't add the moulin input a second time*/
+        if(node->IsClone()) return NULL;
+
+        /*Initialize Element vector*/
         ElementVector* pe=new ElementVector(&node,1,this->parameters);
 
-        this->element->GetInputValue(&moulin_load,node,HydrologyMoulinInputEnum);
-        pe->values[0]=moulin_load;
+        /*Get all inputs and parameters*/
+        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
+        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+        IssmDouble g         = element->FindParam(ConstantsGEnum);
+        IssmDouble Am        = 0.; //For now...
+
+        /*Get hydraulic potential*/
+        IssmDouble phi_old,moulin_load;
+        element->GetInputValue(&phi_old,node,HydraulicPotentialOldEnum);
+        element->GetInputValue(&moulin_load,node,HydrologyMoulinInputEnum);
+
+        pe->values[0] = moulin_load -Am/(rho_water*g) * phi_old/dt;
 
         /*Clean up and return*/

issm/trunk-jpl/src/c/classes/Loads/Moulin.h

@@ -76 +76 @@
                 /*}}}*/
 
+                ElementMatrix* CreateKMatrixHydrologyGlaDS(void);
                 ElementVector* CreatePVectorHydrologyShakti(void);
                 ElementVector* CreatePVectorHydrologyGlaDS(void);