Changeset 15020

Timestamp:

05/13/13 16:45:14 (12 years ago)

Author:

bdef

Message:

NEW: Introduction of the transfer term in the Pvectors

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Tria.cpp (modified) (9 diffs)
• classes/Loads/Pengrid.cpp (modified) (1 diff)
• modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp (modified) (3 diffs)
• solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (1 diff)

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

@@ -2078 +2078 @@
                                 name==SedimentHeadOldEnum ||
                                 name==SedimentHeadEnum ||
+                                name==WaterTransferEnum ||
                                 name==BasisIntegralEnum ||
                                 name==QmuVxEnum ||
@@ -6150 +6151 @@
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble dt,scalar,water_head;
-        IssmDouble water_load;
+        IssmDouble water_load,transfer;
         IssmDouble sediment_storing;
         IssmDouble basis[numdof];
@@ -6163 +6164 @@
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
         Input* water_input=inputs->GetInput(BasalforcingsMeltingRateEnum);  _assert_(water_input);
+        Input* transfer_input=inputs->GetInput(WaterTransferEnum);  _assert_(transfer_input);
         Input* old_wh_input=NULL; 
 
@@ -6181 +6183 @@
                 /*Loading term*/
                 water_input->GetInputValue(&water_load,gauss);
-                scalar = Jdet*gauss->weight*water_load;
+                transfer_input->GetInputValue(&transfer,gauss);
+                scalar = Jdet*gauss->weight*(water_load-transfer);
                 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
                 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
@@ -6207 +6210 @@
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble dt,scalar,water_head;
-        IssmDouble residual_load;
+        IssmDouble residual_load,transfer;
         IssmDouble epl_storing;
         IssmDouble basis[numdof];
@@ -6220 +6223 @@
         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
         Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum);  _assert_(residual_input);
+        Input* transfer_input=inputs->GetInput(WaterTransferEnum);  _assert_(transfer_input);
         Input* old_wh_input=NULL; 
 
@@ -6237 +6241 @@
                 /*Loading term*/
                 residual_input->GetInputValue(&residual_load,gauss);
-                scalar = Jdet*gauss->weight*residual_load;
+                transfer_input->GetInputValue(&transfer,gauss);
+                scalar = Jdet*gauss->weight*(residual_load+transfer);
                 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
                 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
@@ -6350 +6355 @@
                 for(int i=0;i<NUMVERTICES;i++){
                         this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
-                        if(values[i]>h_max)
-                                residual[i]=kappa*(values[i]-h_max)/(dt*intbasis[i]);
+                        if(values[i]>h_max){
+                                residual[i]=kappa*(values[i]-h_max)/intbasis[i];
+                                if(reCast<bool,IssmDouble>(dt))residual[i]=residual[i]/dt;
+                        }
                         else
                                 residual[i]=0.0;
                 }
-        }       
-
+        }
         /*Add input to the element: */
 
@@ -6516 +6522 @@
                         }
                 }
+                break;
         default:
                 _error_("no case higher than 1 for the Transfer method");

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

@@ -703 +703 @@
         element->GetInputValue(&h,node,SedimentHeadEnum);
         element->GetHydrologyDCInefficientHmax(&h_max,node);
-
         if (h>h_max)
          new_active=1;

issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp

@@ -13 +13 @@
         int         hydrology_model;
         int         sedimentlimit_flag;
+        int         transfer_flag;
         bool        isefficientlayer;
         IssmDouble  sedimentlimit;
         IssmDouble  penalty_factor;
+        IssmDouble  leakagefactor;
 
         /*Get parameters: */
@@ -31 +33 @@
         iomodel->FetchData(&isefficientlayer,HydrologydcIsefficientlayerEnum);
         iomodel->FetchData(&sedimentlimit_flag,HydrologydcSedimentlimitFlagEnum);
+        iomodel->FetchData(&transfer_flag,HydrologydcTransferFlagEnum);
         iomodel->FetchData(&penalty_factor,HydrologydcPenaltyFactorEnum);
 
@@ -38 +41 @@
         }
 
+        if(transfer_flag==1){
+                iomodel->FetchData(&leakagefactor,HydrologydcLeakageFactorEnum);
+                parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor));
+        }
+
         parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor));
         parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
         parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
         parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag));
+        parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag));
 
         /*Assign output pointer: */

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp

@@ -45 +45 @@
                 femmodel->UpdateConstraintsx();
                 femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum);
+                femmodel->HydrologyTransferx();
                 for(;;){