Changeset 17444

Timestamp:

03/17/14 08:58:00 (11 years ago)

Author:

bdef

Message:

NEW:Adding moulin Input capability to the Hydro Model

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (4 diffs)
• classes/Loads/Pengrid.cpp (modified) (2 diffs)
• classes/Loads/Pengrid.h (modified) (1 diff)
• modules/SystemMatricesx/SystemMatricesx.cpp (modified) (3 diffs)

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

@@ -124 +124 @@
         iomodel->FetchData(1,MeshVertexonbedEnum);
 
-        //create penalties for nodes: no node can have a temperature over the melting point
+        //create penalties for nodes: no node can have water above the max
         CreateSingleNodeToElementConnectivity(iomodel);
         for(int i=0;i<iomodel->numberofvertices;i++){
-                if (iomodel->meshtype==Mesh3DEnum){
+                if (iomodel->meshtype!=Mesh3DEnum){
                         /*keep only this partition's nodes:*/
                         if(iomodel->my_vertices[i]){
@@ -280 +280 @@
         bool       active_element,isefficientlayer;
         IssmDouble dt,scalar;
-        IssmDouble moulin_load,water_head;
+        IssmDouble water_head;
         IssmDouble water_load,transfer;
         IssmDouble Jdet;
@@ -307 +307 @@
         Input* bed_input         = basalelement->GetInput(BedEnum);
         Input* water_input       = basalelement->GetInput(BasalforcingsMeltingRateEnum);    _assert_(water_input);
-        Input* moulin_input      = basalelement->GetInput(HydrologydcBasalMoulinInputEnum); _assert_(moulin_input);
         if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);             _assert_(old_wh_input);}
 
@@ -327 +326 @@
                 /*Loading term*/
                 water_input->GetInputValue(&water_load,gauss);
-                moulin_input->GetInputValue(&moulin_load,gauss);
-        
+
                 scalar = Jdet*gauss->weight*(water_load);
-                scalar = scalar + Jdet* moulin_load;
                 if(dt!=0.) scalar = scalar*dt;
                 for(int i=0;i<numnodes;i++){

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

@@ -181 +181 @@
 void  Pengrid::CreatePVector(Vector<IssmDouble>* pf){
 
-        /*No loads applied, do nothing: */
-        return;
+        ElementVector* pe=NULL;
+        int analysis_type;
+        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+        switch(analysis_type){
+                
+        case HydrologyDCInefficientAnalysisEnum:
+                pe = CreatePVectorHydrologyDCInefficient();
+                break;
+        default:
+                /*No loads applied, do nothing: */
+                return;
+        }
+        if(pe){
+                pe->AddToGlobal(pf);
+                delete pe;
+        }
 
 }
@@ -793 +808 @@
 }
 /*}}}*/
+/*FUNCTION Pengrid::CreatePVectorHydrologyDCInefficient {{{*/
+ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){
+
+        IssmDouble moulin_load,dt;
+
+        /*Initialize Element matrix*/
+        ElementVector* pe=new ElementVector(&node,1,this->parameters);
+
+        this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
+        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
+        if(dt!=0.0) pe->values[0]=moulin_load*dt;
+
+        /*Clean up and return*/
+        return pe;
+ }
+/*}}}*/
 /*FUNCTION Pengrid::ResetConstraint {{{*/
 void  Pengrid::ResetConstraint(void){

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

@@ -94 +94 @@
                 void           ConstraintActivateHydrologyDCInefficient(int* punstable);
                 void  ConstraintActivate(int* punstable);
+                ElementVector* CreatePVectorHydrologyDCInefficient(void);
                 void  ResetConstraint(void);
                 /*}}}*/

issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp

@@ -76 +76 @@
                 element->ReduceMatrices(Ke,pe);
                 if(Ke) Ke->AddToGlobal(Kff,Kfs);
-                if(pe) pe->AddToGlobal(pf);
+                if(pe){
+                        pe->AddToGlobal(pf);
+                        /* printf("-------------------USUAL \n"); */
+                        /* pf->Echo(); */
+                }
                 delete Ke;
                 delete pe;
@@ -87 +91 @@
                         load->CreateKMatrix(Kff,Kfs);
                         load->CreatePVector(pf);
+                        /* printf("-------------------LOADING \n"); */
+                        /* pf->Echo(); */
                 }
         }
@@ -97 +103 @@
                                 load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
                                 load->PenaltyCreatePVector(pf,kmax);
+                                /* printf("-------------------PENALTY \n"); */
+                                /* pf->Echo(); */
                         }
                 }