Changeset 15214

Timestamp:

06/07/13 14:48:34 (12 years ago)

Author:

bdef

Message:

BUG:fix on the EPL mask definition

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Tria.cpp (modified) (10 diffs)
• solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (2 diffs)

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

@@ -4703 +4703 @@
 
                 /*Compute SurfaceLogVelMisfit:
-                 *                 [        vel + eps     ] 2
+                 *        4         [        vel + eps     ] 2
                  * J = 4 \bar{v}^2 | log ( -----------  ) |  
                  *                 [       vel   + eps    ]
@@ -6178 +6178 @@
                 water_input->GetInputValue(&water_load,gauss);
                 transfer_input->GetInputValue(&transfer,gauss);
-                scalar = Jdet*gauss->weight*(water_load-transfer);
+                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];
@@ -6236 +6236 @@
                 /*Loading term*/
                 transfer_input->GetInputValue(&transfer,gauss);
-                scalar = Jdet*gauss->weight*transfer;
+                scalar = Jdet*gauss->weight*(-transfer);
                 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
                 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
@@ -6418 +6418 @@
         this->inputs->GetInputValue(&converged,ConvergedEnum);
 
+        /*Get inputs*/
+        if(isefficientlayer){
+                GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
+        }
         if(converged){
                 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
                 this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
                 this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
-
-                if(isefficientlayer)GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
+                
                 
                 kappa=kmax*pow(10.,penalty_factor);
@@ -6441 +6444 @@
                 }
         }
+        else{
+                for(int i=0;i<NUMVERTICES;i++){ 
+                        if(isefficientlayer){
+                                if(activeEpl[i]==1.0)activate[i]=1.0;
+                        }                               
+                }
+        }
+
         /*Add input to the element: */
-        if(isefficientlayer) this->inputs->AddInput(new TriaP1Input(HydrologydcMaskEplactiveEnum,activate));
+        if(isefficientlayer) this->inputs->AddInput(new TriaP1Input(HydrologydcMaskEplactiveEnum,&activate[0]));
         this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values));
         this->inputs->AddInput(new TriaP1Input(SedimentHeadResidualEnum,residual));
 
+        /*              if(id==1717){
+                        printf("6 \n");
+                        Input* EplMask=inputs->GetInput(HydrologydcMaskEplactiveEnum);
+                        EplMask->Echo();
+                }
+        */
         /*Free ressources:*/
         xDelete<int>(doflist);
@@ -6506 +6523 @@
                                 _error_("no case higher than 3 for SedimentlimitFlag");
                 }
-
                 /*Assign output pointer*/
                 *ph_max=h_max;
@@ -6520 +6536 @@
         int        transfermethod;
         IssmDouble sed_trans,sed_thick;
-        IssmDouble leakage;
+        IssmDouble leakage,h_max;
         IssmDouble activeEpl[numdof];
         IssmDouble wh_trans[numdof]={0.0};
@@ -6534 +6550 @@
                 /*Also get the flag to the transfer method*/
                 this->parameters->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
-                
+
                 /*Switch between the different transfer methods cases*/
                 switch(transfermethod){
@@ -6552 +6568 @@
                         
                         for(int i=0;i<numdof;i++){
-                        
-                                if(activeEpl[i]==0.0)break;
-                                
+                                if(activeEpl[i]==0.0)continue;
+                
                                 if(sed_head[i]>epl_head[i]){
                                         storing[i]=matpar->GetSedimentStoring();
-                                        wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick);
+                                        wh_trans[i]=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
                                 }
                                 else{
+                                        this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
                                         storing[i]=matpar->GetEplStoring();
-                                        wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick);
+                                        wh_trans[i]=sed_trans*storing[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
+                                        if(sed_head[i]>h_max){
+                                                wh_trans[i]=0.0;
+                                        }
+                                        if((sed_head[i]+wh_trans[i])>h_max){
+                                                wh_trans[i]=h_max-sed_head[i];
+                                        }
                                 }
+                                transfer->SetValue(doflist[i],wh_trans[i],INS_VAL);
                         }
                         break;
@@ -6570 +6593 @@
         }
         /*Assign output pointer*/
-        transfer->SetValues(numdof,doflist,&wh_trans[0],ADD_VAL);
+        /*transfer->SetValues(numdof,doflist,&wh_trans[0],INS_VAL);*/
 }
 

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

@@ -73 +73 @@
                 femmodel->UpdateConstraintsx();
                 femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum);
-                femmodel->HydrologyTransferx();
-
+                /*              femmodel->HydrologyTransferx();
+                 */
                 /*Iteration on the sediment layer*/
                 for(;;){
+                        femmodel->HydrologyTransferx();
                         femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &sediment_kmax);
                         CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
@@ -117 +118 @@
                 /*Iteration on the EPL layer*/
                 for(;;){
+                        femmodel->HydrologyTransferx();
                         femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df,NULL);
                         CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);