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

Index: ../trunk-jpl/src/c/analyses/hydrology_core.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/hydrology_core.cpp      (revision 15005)
+++ ../trunk-jpl/src/c/analyses/hydrology_core.cpp      (revision 15006)
@@ -62,6 +62,7 @@
                        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,7 +78,9 @@
                }
 
                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)){
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp  (revision 15005)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp  (revision 15006)
@@ -42,6 +42,7 @@
        iomodel->FetchDataToInput(elements,MaskElementonwaterEnum);
        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
        iomodel->FetchDataToInput(elements,SedimentHeadEnum);
+       /*      iomodel->FetchDataToInput(elements,SedimentHeadOldEnum);*/
 
        /*Free data: */
        iomodel->DeleteData(1,MeshElementsEnum);
Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp   (revision 15005)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp   (revision 15006)
@@ -14,7 +14,7 @@
        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,15 +36,8 @@
        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,6 +45,7 @@
                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,7 +68,6 @@
 
                        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);
Index: ../trunk-jpl/src/c/classes/objects/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/objects/Elements/Tria.cpp        (revision 15005)
+++ ../trunk-jpl/src/c/classes/objects/Elements/Tria.cpp        (revision 15006)
@@ -2074,7 +2074,9 @@
                                name==MaterialsRheologyZbarEnum ||
                                name==GradientEnum ||
                                name==OldGradientEnum ||
-                name==ConvergedEnum ||
+                               name==ConvergedEnum ||
+                               name==SedimentHeadOldEnum ||
+                               name==SedimentHeadEnum ||
                                name==BasisIntegralEnum ||
                                name==QmuVxEnum ||
                                name==QmuVyEnum ||
@@ -6022,7 +6024,6 @@
                                                                                 &Ke->values[0],1);
                }
        }
-
        /*Clean up and return*/
        delete gauss;
        return Ke;
@@ -6164,11 +6165,12 @@
        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,7 +6191,6 @@
                        for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
                }
        }
-
        /*Clean up and return*/
        delete gauss;
        return pe;
@@ -6339,7 +6340,7 @@
        /*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,8 +6357,10 @@
        }       
 
        /*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);