Changeset 23546

Timestamp:

12/13/18 07:05:22 (6 years ago)

Author:

bdef

Message:

CHG:clean up of old hydro substeping and change in reset of inactive EPL

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCEfficientAnalysis.cpp (modified) (4 diffs)
• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (1 diff)
• classes/FemModel.cpp (modified) (1 diff)
• classes/FemModel.h (modified) (1 diff)
• classes/Inputs/TransientInput.cpp (modified) (4 diffs)
• classes/Inputs/TransientInput.h (modified) (2 diffs)
• cores/hydrology_core.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp ¶

@@ -198 +198 @@
         /*Check that all nodes are active, else return empty matrix*/
         if(!active_element) {
-        if(domaintype!=Domain2DhorizontalEnum){
+                if(domaintype!=Domain2DhorizontalEnum){
                         basalelement->DeleteMaterials();
                         delete basalelement;
@@ -435 +435 @@
         IssmDouble* eplHeads    = xNew<IssmDouble>(numnodes);
         IssmDouble* basevalue    = xNew<IssmDouble>(numnodes);
-        IssmDouble* initvalue    = xNew<IssmDouble>(numnodes);
         basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 
@@ -452 +451 @@
         }
         else{
-                //tradeof between keeping initial condition and not getting to far from head at deactivation
+                //Fixing epl head at bedrock elevation when deactivating
                 basalelement->GetInputListOnVertices(&basevalue[0],BaseEnum);
-                basalelement->GetInputListOnVertices(&initvalue[0],EplHeadHydrostepEnum);
                 for(int i=0;i<numnodes;i++){
-                        eplHeads[i]=max(basevalue[i],initvalue[i]);
+                        eplHeads[i]=basevalue[i];
                         if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
                         if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector");
@@ -466 +464 @@
         xDelete<IssmDouble>(eplHeads);
         xDelete<IssmDouble>(basevalue);
-        xDelete<IssmDouble>(initvalue);
         xDelete<int>(doflist);
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

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

@@ -205 +205 @@
         /*Check that all nodes are active, else return empty matrix*/
         if(!thawed_element) {
-        if(domaintype!=Domain2DhorizontalEnum){
+                if(domaintype!=Domain2DhorizontalEnum){
                         basalelement->DeleteMaterials();
                         delete basalelement;

issm/trunk-jpl/src/c/classes/FemModel.cpp ¶

@@ -4821 +4821 @@
 }
 /*}}}*/
-void FemModel::InitTransientOutputx(int* stackedinput_enum,int numoutputs){ /*{{{*/
-
-  for(int i=0;i<numoutputs;i++){
-                if(stackedinput_enum[i]<0){
-                        _error_("Can't deal with non enum fields for result Stack");
-                }
-                else{
-                        for(int j=0;j<elements->Size();j++){
-                                TransientInput* transient_input = new TransientInput(stackedinput_enum[i]);
-                                /*Intermediaries*/
-                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
-                                //transient_input->Configure(element->parameters);
-                                element->inputs->AddInput(transient_input);
-                        }
-                }
-        }
-}
-/*}}}*/
-void FemModel::StackTransientOutputx(int* input_enum,int* stackedinput_enum,IssmDouble hydrotime,int numoutputs){ /*{{{*/
-
-  for(int i=0;i<numoutputs;i++){
-                if(input_enum[i]<0){
-                        _error_("Can't deal with non enum fields for result Stack");
-                }
-                else{
-                        for(int j=0;j<elements->Size();j++){
-                                /*Intermediaries*/
-                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
-                                TransientInput* stacking_input=NULL;
-                                Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input); //this is the enum stack
-                                stacking_input=dynamic_cast<TransientInput*>(input);
-
-                                int  numvertices = element->GetNumberOfVertices();
-                                IssmDouble* N=xNew<IssmDouble>(numvertices);
-                                element->GetInputListOnVertices(&N[0],input_enum[i]);   //this is the enum to stack
-                                switch(element->ObjectEnum()){
-                                case TriaEnum:
-                                        stacking_input->AddTimeInput(new TriaInput(stackedinput_enum[i],&N[0],P1Enum),hydrotime);
-                                        break;
-                                case PentaEnum:
-                                        stacking_input->AddTimeInput(new PentaInput(stackedinput_enum[i],&N[0],P1Enum),hydrotime);
-                                        break;
-                                case TetraEnum:
-                                        stacking_input->AddTimeInput(new TetraInput(stackedinput_enum[i],&N[0],P1Enum),hydrotime);
-                                        break;
-                                default: _error_("Not implemented yet");
-                                }
-                                xDelete<IssmDouble>(N);
-                        }
-                }
-        }
-}
-/*}}}*/
-void FemModel::AverageTransientOutputx(int* stackedinput_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs){ /*{{{*/
-
-  for(int i=0;i<numoutputs;i++){
-                if(stackedinput_enum[i]<0){
-                        _error_("Can't deal with non enum fields for result Stack");
-                }
-                else{
-                        for(int j=0;j<elements->Size();j++){
-                                /*Intermediaries*/
-                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
-                                int  numvertices = element->GetNumberOfVertices();
-                                IssmDouble* time_averaged=NULL;
-
-                                Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input);
-                                TransientInput* stacking_input=NULL;
-                                stacking_input=dynamic_cast<TransientInput*>(input);
-                                stacking_input->GetInputAverageOnTimes(&time_averaged,init_time);
-
-                                element->AddInput(averagedinput_enum[i],&time_averaged[0],P1Enum);
-                                xDelete<IssmDouble>(time_averaged);
-                        }
-                }
-        }
-}
-/*}}}*/
 void FemModel::InitMeanOutputx(int* stackedinput_enum,int numoutputs){ /*{{{*/
 

issm/trunk-jpl/src/c/classes/FemModel.h ¶

@@ -157 +157 @@
                 void UpdateConstraintsExtrudeFromTopx();
                 void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count);
-                void InitTransientOutputx(int* input_enum, int numoutputs);
-                void StackTransientOutputx(int* input_enum,int* stackedinput_enum, IssmDouble hydrotime, int numoutputs);
-                void AverageTransientOutputx(int* stackedinput_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs);
                 void InitMeanOutputx(int* stackedinput_enum,int numoutputs);
                 void SumOutputx(int* input_enum,int* stackedinput_enum,int numoutputs);

issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp ¶

@@ -87 +87 @@
         _printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
         _printf_("   numtimesteps: " << this->numtimesteps << "\n");
-        _printf_("---inputs: \n"); 
+        _printf_("---inputs: \n");
         for(i=0;i<this->numtimesteps;i++){
                 _printf_("   time: " << this->timesteps[i]<<"  ");
@@ -313 +313 @@
 }
 /*}}}*/
-void TransientInput::GetInputAverageOnTimes(IssmDouble** pvalues, IssmDouble init_time){/*{{{*/
-
-        int                                     numnodes;
-        IssmDouble  subtimestep;
-        IssmDouble* values= NULL;
-        Input* input=NULL;
-
-        /*Initialize numnode from transientinput out of time loop*/
-        for(int i=0;i<this->numtimesteps;i++){
-                /*First compute the lengt of the sub-timestep*/
-                if(i==0){
-                        subtimestep = this->timesteps[i]-init_time;
-                }
-                else{
-                        subtimestep = this->timesteps[i]-this->timesteps[i-1];
-                }       
-                /*Figure out type of input and get it*/
-                input = xDynamicCast<Input*>(this->inputs->GetObjectByOffset(i)); _assert_(input);
-                switch(input->ObjectEnum()){
-                case TriaInputEnum:{
-                        TriaInput* triainput = (TriaInput*)this->inputs->GetObjectByOffset(i); _assert_(triainput);
-                        if(i==0){
-                                numnodes= triainput->NumberofNodes(triainput->interpolation_type);
-                                values=xNewZeroInit<IssmDouble>(numnodes);
-                        }
-                        /*Sum the values weighted by their respective sub-timestep*/
-                        for(int j=0;j<numnodes;j++){
-                                values[j]+=(triainput->values[j]*subtimestep);
-                        }
-                        break;
-                }
-                case PentaInputEnum:{
-                        PentaInput*     pentainput = (PentaInput*)this->inputs->GetObjectByOffset(i); _assert_(pentainput);
-                        if(i==0){
-                                numnodes= pentainput->NumberofNodes(pentainput->interpolation_type);
-                                values=xNewZeroInit<IssmDouble>(numnodes);
-                        }
-                        /*Sum the values weighted by their respective sub-timestep*/
-                        for(int j=0;j<numnodes;j++){
-                                values[j]+=(pentainput->values[j]*subtimestep);
-                        }
-                        break;
-                }
-                case TetraInputEnum:{
-                        TetraInput*     tetrainput = (TetraInput*)this->inputs->GetObjectByOffset(i); _assert_(tetrainput);
-                        if(i==0){
-                                numnodes= tetrainput->NumberofNodes(tetrainput->interpolation_type);
-                                values=xNewZeroInit<IssmDouble>(numnodes);
-                        }
-                        /*Sum the values weighted by their respective sub-timestep*/
-                        for(int j=0;j<numnodes;j++){
-                                values[j]+=(tetrainput->values[j]*subtimestep);
-                        }
-                        break;
-                }
-                default: _error_("not implemented yet");
-                }
-        }
-        /*Compute the average based on the length of the full timestep*/
-        for(int j=0;j<numnodes;j++){
-                values[j]=values[j]/(this->timesteps[this->numtimesteps-1]-init_time);
-        }
-        *pvalues=values;
-}
-/*}}}*/
 
 /*Intermediary*/
@@ -476 +411 @@
         Input *input2 = NULL;
 
-        /*go through the timesteps, and figure out which interval we 
+        /*go through the timesteps, and figure out which interval we
          *fall within. Then interpolate the values on this interval: */
         found=binary_search(&offset,intime,this->timesteps,this->numtimesteps);
@@ -498 +433 @@
                 alpha1=(1.0-alpha2);
 
-                input1=(Input*)this->inputs->GetObjectByOffset(offset); 
+                input1=(Input*)this->inputs->GetObjectByOffset(offset);
                 input2=(Input*)this->inputs->GetObjectByOffset(offset+1);
 

issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h ¶

@@ -1 @@
-/*! \file TransientInput.h 
+/*! \file TransientInput.h
  *  \brief: header file for transientinput object
  */
@@ -64 +64 @@
                 void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");};
                 void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime);
-                void GetInputAverageOnTimes(IssmDouble** pvalues, IssmDouble init_time);
                 int  GetInputInterpolationType(){_error_("not implemented yet!");};
                 IssmDouble  GetTimeByOffset(int offset);

issm/trunk-jpl/src/c/cores/hydrology_core.cpp ¶

@@ -75 +75 @@
                                         int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum};
                                         int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum};
-                                        //femmodel->InitTransientOutputx(&stackedinput[0],4);
                                         femmodel->InitMeanOutputx(&stackedinput[0],4);
                                         while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
@@ -89 +88 @@
                                                 solutionsequence_hydro_nonlinear(femmodel);
                                                 /*If we have a sub-timestep we stack the variables here*/
-                                                //femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4);
                                                 femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],4);
                                         }
-                                        //femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4);
                                         femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],4);
                                 }
@@ -99 +96 @@
                                         int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum};
                                         int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum};
-                                        //femmodel->InitTransientOutputx(&stackedinput[0],2);
                                         femmodel->InitMeanOutputx(&stackedinput[0],2);
                                         while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
@@ -108 +104 @@
                                                 solutionsequence_hydro_nonlinear(femmodel);
                                                 /*If we have a sub-timestep we stack the variables here*/
-                                                //femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2);
                                                 femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],2);
                                         }
-                                        //femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2);
                                         femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],2);
                                 }