source: issm/oecreview/Archive/14312-15392/ISSM-15005-15006.diff@ 15393

../trunk-jpl/src/c/analyses/hydrology_core.cpp

@@ -62 +62 @@
                         solutionsequence_nonlinear(femmodel,modify_loads);
 
                         /*transfer water column thickness to old water column thickness: */
+                        
                         InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WatercolumnEnum,WaterColumnOldEnum);
 
                         if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
@@ -77 +78 @@
                 }
 
                 else if (hydrology_model==HydrologydcEnum){
+                                InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SedimentHeadEnum,SedimentHeadOldEnum);
                         femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
+                        
                         if(VerboseSolution()) _pprintLine_("   computing water head");
                         solutionsequence_hydro_nonlinear(femmodel);
                         if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){

../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp

@@ -42 +42 @@
         iomodel->FetchDataToInput(elements,MaskElementonwaterEnum);
         iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
         iomodel->FetchDataToInput(elements,SedimentHeadEnum);
+        /*      iomodel->FetchDataToInput(elements,SedimentHeadOldEnum);*/
 
         /*Free data: */
         iomodel->DeleteData(1,MeshElementsEnum);

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

@@ -14 +14 @@
         Vector<IssmDouble>* uf=NULL; 
         Vector<IssmDouble>* uf_old=NULL; 
         Vector<IssmDouble>* ys=NULL; 
-        IssmDouble sediment_kmax;
+        IssmDouble sediment_kmax,time;
 
         /*intermediary: */
         Matrix<IssmDouble>* Kff=NULL;
@@ -36 +36 @@
         count=1;
         converged=false;
 
+        femmodel->parameters->FindParam(&time,TimeEnum);
         for(;;){
-                
-                /*Computing the transfer term
-                        
-
-                 */
-
-                
-
                 /*First layer*/
                 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
                 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum);
@@ -52 +45 @@
                 femmodel->UpdateConstraintsx();
                 femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum);
                 for(;;){
+
                         femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &sediment_kmax);
                         CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
                         Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf;
@@ -74 +68 @@
 
                         if(converged){
                                 femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum);
-
                                 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum);
                                 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
                                 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum);

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

@@ -2074 +2074 @@
                                 name==MaterialsRheologyZbarEnum ||
                                 name==GradientEnum ||
                                 name==OldGradientEnum ||
-                 name==ConvergedEnum ||
+                                name==ConvergedEnum ||
+                                name==SedimentHeadOldEnum ||
+                                name==SedimentHeadEnum ||
                                 name==BasisIntegralEnum ||
                                 name==QmuVxEnum ||
                                 name==QmuVyEnum ||
@@ -6022 +6024 @@
                                                                                  &Ke->values[0],1);
                 }
         }
-
         /*Clean up and return*/
         delete gauss;
         return Ke;
@@ -6164 +6165 @@
         Input* old_wh_input=NULL; 
 
         if(reCast<bool,IssmDouble>(dt)){
-                old_wh_input=inputs->GetInput(SedimentHeadEnum); _assert_(old_wh_input);
+                old_wh_input=inputs->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input);
         }
 
         /* Start  looping on the number of gaussian points: */
         gauss=new GaussTria(2);
+        
         for(int ig=gauss->begin();ig<gauss->end();ig++){
 
                 gauss->GaussPoint(ig);
@@ -6189 +6191 @@
                         for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
                 }
         }
-
         /*Clean up and return*/
         delete gauss;
         return pe;
@@ -6339 +6340 @@
         /*If converged keep the residual in mind*/
         this->inputs->GetInputValue(&converged,ConvergedEnum);
         GetInputListOnVertices(&intbasis[0],BasisIntegralEnum);
-        
+
         if(converged){
                 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
                 this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
@@ -6356 +6357 @@
         }       
 
         /*Add input to the element: */
+
         this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values));
         this->inputs->AddInput(new TriaP1Input(SedimentHeadResidualEnum,residual));
+        if(converged)this->inputs->AddInput(new TriaP1Input(SedimentHeadOldEnum,values));
 
         /*Free ressources:*/
         xDelete<int>(doflist);