Changeset 24080


Ignore:
Timestamp:
07/11/19 08:10:15 (6 years ago)
Author:
Mathieu Morlighem
Message:

NEW: improved GlaDS solution sequence to be consistent with Elmer

Location:
issm/trunk-jpl/src/c
Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r24060 r24080  
    307307                                        ./solutionsequences/solutionsequence_hydro_nonlinear.cpp\
    308308                                        ./solutionsequences/solutionsequence_shakti_nonlinear.cpp\
     309                                        ./solutionsequences/solutionsequence_glads_nonlinear.cpp\
    309310                                        ./cores/stressbalance_core.cpp\
    310311                                        ./solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp\

  • issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp

    r24069 r24080  
    55#include "../modules/modules.h"
    66
    7 #define AEPS 2.22e-15
     7#define AEPS 2.2204460492503131E-015
    88
    99/*Model processing*/
    … …  
    423423
    424424/*GlaDS specifics*/
     425void HydrologyGlaDSAnalysis::SetChannelCrossSectionOld(FemModel* femmodel){/*{{{*/
     426
     427        bool ischannels;
     428        femmodel->parameters->FindParam(&ischannels,HydrologyIschannelsEnum);
     429        if(!ischannels) return;
     430
     431        for(int i=0;i<femmodel->loads->Size();i++){
     432                if(femmodel->loads->GetEnum(i)==ChannelEnum){
     433                        Channel* channel=(Channel*)femmodel->loads->GetObjectByOffset(i);
     434                        channel->SetChannelCrossSectionOld();
     435                }
     436        }
     437
     438}/*}}}*/
    425439void HydrologyGlaDSAnalysis::UpdateSheetThickness(FemModel* femmodel){/*{{{*/
    426440
    … …  
    458472        Input* H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
    459473        Input* b_input  = element->GetInput(BedEnum); _assert_(b_input);
    460         Input* h_input  = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
     474        Input* hold_input  = element->GetInput(HydrologySheetThicknessOldEnum);_assert_(hold_input);
    461475        Input* B_input  = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    462476        Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
    … …  
    472486                vx_input->GetInputValue(&vx,gauss);
    473487                vy_input->GetInputValue(&vy,gauss);
    474                 h_input->GetInputValue(&h_old,gauss);
     488                hold_input->GetInputValue(&h_old,gauss);
    475489                B_input->GetInputValue(&B,gauss);
    476490                n_input->GetInputValue(&n,gauss);

  • issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h

    r23965 r24080  
    3232
    3333                /*Specific to GlaDS*/
     34                void SetChannelCrossSectionOld(FemModel* femmodel);
    3435                void UpdateSheetThickness(FemModel* femmodel);
    3536                void UpdateSheetThickness(Element*  element);

  • issm/trunk-jpl/src/c/classes/Loads/Channel.cpp

    r24070 r24080  
    2525#define ALPHA_S     5./4.
    2626#define BETA_S      3./2.
    27 #define AEPS 2.22e-15
     27#define AEPS        2.2204460492503131E-015
    2828
    2929/*Channel constructors and destructor*/
    … …  
    670670        IssmDouble qc = - ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids;
    671671
    672         /*Compute f factor*/
    673         IssmDouble fFactor = 0.;
    674         if(this->S>0. || qc*dPw>0.){
    675                 fFactor = lc * qc;
    676         }
    677 
     672        /*Ice rate factor*/
    678673        A=pow(B,-n);
    679674
    … …  
    682677        IssmDouble N = phi_0 - phi;
    683678
    684         IssmDouble alpha = 1./(rho_ice*L)*(
    685                                 fabs(Qprime*pow(this->S,ALPHA_C-1.)*dphids)
    686                                 + C*Qprime*pow(this->S,ALPHA_C-1.)*dPw
    687                                 ) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
    688 
    689         IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw );
    690 
    691         /*Solve ODE*/
    692         this->S = ODE1(alpha,beta,this->Sold,dt,1);
    693 
    694         /*Make sure Area > 0*/
    695         if(this->S<0.) this->S = 0.;
     679        bool converged  = false;
     680        int  count      = 0;
     681
     682        while(!converged){
     683
     684                IssmDouble Snew = this->S;
     685
     686                /*Compute f factor*/
     687                IssmDouble fFactor = 0.;
     688                if(this->S>0. || qc*dPw>0.){
     689                        fFactor = lc * qc;
     690                }
     691
     692                IssmDouble alpha = 1./(rho_ice*L)*(
     693                                        fabs(Qprime*pow(Snew,ALPHA_C-1.)*dphids)
     694                                        + C*Qprime*pow(Snew,ALPHA_C-1.)*dPw
     695                                        ) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
     696
     697                IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw );
     698
     699                /*Solve ODE*/
     700                this->S = ODE1(alpha,beta,this->Sold,dt,2);
     701                if(this->S<0.) this->S = 0.;
     702                _assert_(!xIsNan<IssmDouble>(this->S));
     703
     704                count++;
     705
     706                if(fabs((this->S - Snew)/(Snew+AEPS))<1e-8  || count>=10) converged = true;
     707        }
    696708
    697709        /*Clean up and return*/

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

    r23965 r24080  
    195195                HydrologyGlaDSAnalysis* analysis = new HydrologyGlaDSAnalysis();
    196196                femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum);
     197
     198                /*Set fields as old*/
    197199                InputDuplicatex(femmodel,HydraulicPotentialEnum,HydraulicPotentialOldEnum);
    198                 analysis->UpdateEffectivePressure(femmodel);
    199                 solutionsequence_shakti_nonlinear(femmodel);
     200                InputDuplicatex(femmodel,HydrologySheetThicknessEnum,HydrologySheetThicknessOldEnum);
     201                analysis->SetChannelCrossSectionOld(femmodel);
     202
     203                /*Solve for new potential*/
     204                solutionsequence_glads_nonlinear(femmodel);
    200205
    201206                if(VerboseSolution()) _printf0_("   updating effective pressure\n");
    202207                analysis->UpdateEffectivePressure(femmodel);
    203                 if(VerboseSolution()) _printf0_("   updating sheet thickness\n"); /*Uses N, so needs to come after*/
    204                 analysis->UpdateSheetThickness(femmodel);
    205                 if(VerboseSolution()) _printf0_("   updating channels cross section\n");
    206                 analysis->UpdateChannelCrossSection(femmodel);
    207208                delete analysis;
    208209        }

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r24061 r24080  
    570570syn keyword cConstant HydrologySheetConductivityEnum
    571571syn keyword cConstant HydrologySheetThicknessEnum
     572syn keyword cConstant HydrologySheetThicknessOldEnum
    572573syn keyword cConstant HydrologyWaterVxEnum
    573574syn keyword cConstant HydrologyWaterVyEnum
    … …  
    941942syn keyword cConstant ChannelEnum
    942943syn keyword cConstant ChannelAreaEnum
     944syn keyword cConstant ChannelAreaOldEnum
    943945syn keyword cConstant ClosedEnum
    944946syn keyword cConstant ColinearEnum
    … …  
    12771279syn keyword cType Cfsurfacesquare
    12781280syn keyword cType Channel
     1281syn keyword cType classes
    12791282syn keyword cType Constraint
    12801283syn keyword cType Constraints
    … …  
    12831286syn keyword cType ControlInput
    12841287syn keyword cType Covertree
     1288syn keyword cType DatasetInput
    12851289syn keyword cType DataSetParam
    1286 syn keyword cType DatasetInput
    12871290syn keyword cType Definition
    12881291syn keyword cType DependentObject
    … …  
    12971300syn keyword cType ElementHook
    12981301syn keyword cType ElementMatrix
     1302syn keyword cType Elements
    12991303syn keyword cType ElementVector
    1300 syn keyword cType Elements
    13011304syn keyword cType ExponentialVariogram
    13021305syn keyword cType ExternalResult
    … …  
    13051308syn keyword cType Friction
    13061309syn keyword cType Gauss
     1310syn keyword cType GaussianVariogram
     1311syn keyword cType gaussobjects
    13071312syn keyword cType GaussPenta
    13081313syn keyword cType GaussSeg
    13091314syn keyword cType GaussTetra
    13101315syn keyword cType GaussTria
    1311 syn keyword cType GaussianVariogram
    13121316syn keyword cType GenericExternalResult
    13131317syn keyword cType GenericOption
    … …  
    13171321syn keyword cType Input
    13181322syn keyword cType Inputs
     1323syn keyword cType IntArrayInput
    13191324syn keyword cType IntInput
    13201325syn keyword cType IntMatParam
    … …  
    13241329syn keyword cType IssmDirectApplicInterface
    13251330syn keyword cType IssmParallelDirectApplicInterface
     1331syn keyword cType krigingobjects
    13261332syn keyword cType Load
    13271333syn keyword cType Loads
    … …  
    13341340syn keyword cType Matice
    13351341syn keyword cType Matlitho
     1342syn keyword cType matrixobjects
    13361343syn keyword cType MatrixParam
    13371344syn keyword cType Misfit
    … …  
    13461353syn keyword cType Observations
    13471354syn keyword cType Option
     1355syn keyword cType Options
    13481356syn keyword cType OptionUtilities
    1349 syn keyword cType Options
    13501357syn keyword cType Param
    13511358syn keyword cType Parameters
    … …  
    13601367syn keyword cType Regionaloutput
    13611368syn keyword cType Results
     1369syn keyword cType Riftfront
    13621370syn keyword cType RiftStruct
    1363 syn keyword cType Riftfront
    13641371syn keyword cType Seg
    13651372syn keyword cType SegInput
     1373syn keyword cType Segment
    13661374syn keyword cType SegRef
    1367 syn keyword cType Segment
    13681375syn keyword cType SpcDynamic
    13691376syn keyword cType SpcStatic
    … …  
    13851392syn keyword cType Vertex
    13861393syn keyword cType Vertices
    1387 syn keyword cType classes
    1388 syn keyword cType gaussobjects
    1389 syn keyword cType krigingobjects
    1390 syn keyword cType matrixobjects
    13911394syn keyword cType AdjointBalancethickness2Analysis
    13921395syn keyword cType AdjointBalancethicknessAnalysis
    … …  
    14071410syn keyword cType FreeSurfaceBaseAnalysis
    14081411syn keyword cType FreeSurfaceTopAnalysis
     1412syn keyword cType GiaIvinsAnalysis
    14091413syn keyword cType GLheightadvectionAnalysis
    1410 syn keyword cType GiaIvinsAnalysis
    14111414syn keyword cType HydrologyDCEfficientAnalysis
    14121415syn keyword cType HydrologyDCInefficientAnalysis

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r24061 r24080  
    566566        HydrologySheetConductivityEnum,
    567567        HydrologySheetThicknessEnum,
     568        HydrologySheetThicknessOldEnum,
    568569        HydrologyWaterVxEnum,
    569570        HydrologyWaterVyEnum,
    … …  
    939940        ChannelEnum,
    940941        ChannelAreaEnum,
     942        ChannelAreaOldEnum,
    941943        ClosedEnum,
    942944        ColinearEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r24061 r24080  
    572572                case HydrologySheetConductivityEnum : return "HydrologySheetConductivity";
    573573                case HydrologySheetThicknessEnum : return "HydrologySheetThickness";
     574                case HydrologySheetThicknessOldEnum : return "HydrologySheetThicknessOld";
    574575                case HydrologyWaterVxEnum : return "HydrologyWaterVx";
    575576                case HydrologyWaterVyEnum : return "HydrologyWaterVy";
    … …  
    943944                case ChannelEnum : return "Channel";
    944945                case ChannelAreaEnum : return "ChannelArea";
     946                case ChannelAreaOldEnum : return "ChannelAreaOld";
    945947                case ClosedEnum : return "Closed";
    946948                case ColinearEnum : return "Colinear";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r24061 r24080  
    584584              else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
    585585              else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
     586              else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum;
    586587              else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
    587588              else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
    … …  
    628629              else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
    629630              else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
    630               else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Misfit")==0) return MisfitEnum;
     634              if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
     635              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
    635636              else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
    636637              else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    … …  
    751752              else if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
    752753              else if (strcmp(name,"SmbV")==0) return SmbVEnum;
    753               else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
     757              if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
     758              else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
    758759              else if (strcmp(name,"SmbW")==0) return SmbWEnum;
    759760              else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
    … …  
    874875              else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
    875876              else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum;
    876               else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum;
     880              if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum;
     881              else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum;
    881882              else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
    882883              else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
    … …  
    964965              else if (strcmp(name,"Channel")==0) return ChannelEnum;
    965966              else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
     967              else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum;
    966968              else if (strcmp(name,"Closed")==0) return ClosedEnum;
    967969              else if (strcmp(name,"Colinear")==0) return ColinearEnum;
    … …  
    996998              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
    997999              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    998               else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
    999               else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
     1003              if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
     1004              else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
     1005              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
    10041006              else if (strcmp(name,"Element")==0) return ElementEnum;
    10051007              else if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
    … …  
    11191121              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    11201122              else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
    1121               else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    1122               else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
     1126              if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
     1127              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
     1128              else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
    11271129              else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
    11281130              else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
    … …  
    12421244              else if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum;
    12431245              else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum;
    1244               else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum;
    1245               else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
     1249              if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum;
     1250              else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
     1251              else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
    12501252              else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
    12511253              else if (strcmp(name,"Tetra")==0) return TetraEnum;

  • issm/trunk-jpl/src/c/solutionsequences/convergence.cpp

    r23205 r24080  
    77#include "../shared/shared.h"
    88
    9 void convergence(bool* pconverged, Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Vector<IssmDouble>* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs){/*{{{*/
     9void convergence(bool* pconverged, Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Vector<IssmDouble>* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs){
    1010
    1111        /*output*/
    … …  
    139139        /*assign output*/
    140140        *pconverged=converged;
    141 }/*}}}*/
    142 
    143 
     141}

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_shakti_nonlinear.cpp

    r23587 r24080  
    3939        Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
    4040
    41         //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
    42         InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
    43 
    4441        for(;;){
    4542                /*save pointer to old solution*/
    … …  
    5653
    5754                convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df;
    58                 InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
    5955                InputUpdateFromSolutionx(femmodel,ug);
    6056

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h

    r23161 r24080  
    1515void solutionsequence_hydro_nonlinear(FemModel* femmodel);
    1616void solutionsequence_shakti_nonlinear(FemModel* femmodel);
     17void solutionsequence_glads_nonlinear(FemModel* femmodel);
    1718void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);
    1819void solutionsequence_newton(FemModel* femmodel);
Note: See TracChangeset for help on using the changeset viewer.