Changeset 24097


Ignore:
Timestamp:
07/17/19 01:18:19 (6 years ago)
Author:
Mathieu Morlighem
Message:

BUG: fixing new mass transport scheme

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r24091 r24097  
    808808        if(!element->IsOnBase()) return;
    809809
    810         /*Get basal element*/
    811         int domaintype; element->FindParam(&domaintype,DomainTypeEnum);
    812    Element* basalelement=element;
    813         if(domaintype!=Domain2DhorizontalEnum) basalelement = element->SpawnBasalElement();
    814 
    815         /*Fetch number of nodes and dof for this finite element*/
    816         int numnodes    = basalelement->GetNumberOfNodes();
    817    int numvertices = basalelement->GetNumberOfVertices();
    818 
    819    /*Keep old thickness for later*/
    820         IssmDouble* oldthickness = xNew<IssmDouble>(numvertices);
    821    basalelement->GetInputListOnNodes(&oldthickness[0],ThicknessEnum);
    822 
    823810        /*Fetch dof list and allocate solution vector*/
    824811        int *doflist = NULL;
    825812        element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
     813
     814        int numnodes = element->GetNumberOfNodes();
    826815        IssmDouble* newthickness = xNew<IssmDouble>(numnodes);
    827816
    828817        /*Use the dof list to index into the solution vector: */
    829    IssmDouble minthickness = basalelement->FindParam(MasstransportMinThicknessEnum);
     818   IssmDouble minthickness = element->FindParam(MasstransportMinThicknessEnum);
    830819        for(int i=0;i<numnodes;i++){
    831820                newthickness[i]=solution[doflist[i]];
    … …  
    835824      if(newthickness[i]<minthickness) newthickness[i]=minthickness;
    836825        }
    837 
    838826        element->AddBasalInput(ThicknessEnum,newthickness,element->GetElementType());
    839 
    840         /*Free ressources:*/
     827        xDelete<int>(doflist);
    841828        xDelete<IssmDouble>(newthickness);
    842         xDelete<int>(doflist);
     829
     830        /*Update bed and surface accordingly*/
     831
     832        /*Get basal element*/
     833        int domaintype; element->FindParam(&domaintype,DomainTypeEnum);
     834   Element* basalelement=element;
     835        if(domaintype!=Domain2DhorizontalEnum) basalelement = element->SpawnBasalElement();
     836
     837        /*Fetch number of nodes and dof for this finite element*/
     838   int numvertices = basalelement->GetNumberOfVertices();
    843839
    844840        /*Now, we need to do some "processing"*/
    845841   newthickness  = xNew<IssmDouble>(numvertices);
     842        IssmDouble* oldthickness      = xNew<IssmDouble>(numvertices);
    846843        IssmDouble* cumdeltathickness = xNew<IssmDouble>(numvertices);
    847844        IssmDouble* deltathickness    = xNew<IssmDouble>(numvertices);
    … …  
    856853        /*Get previous base, thickness, surfac and current sealevel and bed:*/
    857854   basalelement->GetInputListOnVertices(&newthickness[0],ThicknessEnum);
     855        basalelement->GetInputListOnVertices(&oldthickness[0],ThicknessOldEnum);
    858856        basalelement->GetInputListOnVertices(&oldbase[0],BaseEnum);
    859857        basalelement->GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
    … …  
    869867        /*What is the delta thickness forcing the sea-level rise core: cumulated over time, hence the +=:*/
    870868        for(int i=0;i<numvertices;i++){
    871                 cumdeltathickness[i]+=newthickness[i]-oldthickness[i];
    872                 deltathickness[i]=newthickness[i]-oldthickness[i];
     869                cumdeltathickness[i] += newthickness[i]-oldthickness[i];
     870                deltathickness[i]     = newthickness[i]-oldthickness[i];
    873871        }
    874872
    … …  
    880878
    881879        for(int i=0;i<numvertices;i++) {
    882                 if (phi[i]>0.){ //this is an ice sheet: just add thickness to base.
    883                         /*Update! actually, the bed has moved too!:*/
     880                if (phi[i]>0.){ //this is grounded ice: just add thickness to base.
    884881                        if(isgroundingline){
    885882                                newsurface[i] = bed[i]+newthickness[i]; //surface = bed + newthickness
    … …  
    888885                        else{
    889886                                 newsurface[i] = oldbase[i]+newthickness[i]; //surface = oldbase + newthickness
    890                                  newbase[i]     = oldbase[i];                 //same base: do nothing
     887                                 newbase[i]    = oldbase[i];                 //same base: do nothing
    891888                        }
    892889                }

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r24089 r24097  
    18081808        /*Fetch dof list and allocate solution vector*/
    18091809        GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
    1810         IssmDouble* values    = xNew<IssmDouble>(numnodes);
     1810        IssmDouble* values = xNew<IssmDouble>(numnodes);
    18111811
    18121812        /*Use the dof list to index into the solution vector: */

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

    r23221 r24097  
    6161                }
    6262                femmodel->SetCurrentConfiguration(MasstransportAnalysisEnum);
     63                InputDuplicatex(femmodel,ThicknessEnum,ThicknessOldEnum);
    6364                if(stabilization==4){
    6465                        solutionsequence_fct(femmodel);

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

    r24088 r24097  
    781781syn keyword cConstant ThicknessAlongGradientEnum
    782782syn keyword cConstant ThicknessEnum
     783syn keyword cConstant ThicknessOldEnum
    783784syn keyword cConstant ThicknessPositiveEnum
    784785syn keyword cConstant VelEnum
    … …  
    12811282syn keyword cType Cfsurfacesquare
    12821283syn keyword cType Channel
     1284syn keyword cType classes
    12831285syn keyword cType Constraint
    12841286syn keyword cType Constraints
    … …  
    12871289syn keyword cType ControlInput
    12881290syn keyword cType Covertree
     1291syn keyword cType DatasetInput
    12891292syn keyword cType DataSetParam
    1290 syn keyword cType DatasetInput
    12911293syn keyword cType Definition
    12921294syn keyword cType DependentObject
    … …  
    13011303syn keyword cType ElementHook
    13021304syn keyword cType ElementMatrix
     1305syn keyword cType Elements
    13031306syn keyword cType ElementVector
    1304 syn keyword cType Elements
    13051307syn keyword cType ExponentialVariogram
    13061308syn keyword cType ExternalResult
    … …  
    13091311syn keyword cType Friction
    13101312syn keyword cType Gauss
     1313syn keyword cType GaussianVariogram
     1314syn keyword cType gaussobjects
    13111315syn keyword cType GaussPenta
    13121316syn keyword cType GaussSeg
    13131317syn keyword cType GaussTetra
    13141318syn keyword cType GaussTria
    1315 syn keyword cType GaussianVariogram
    13161319syn keyword cType GenericExternalResult
    13171320syn keyword cType GenericOption
    … …  
    13211324syn keyword cType Input
    13221325syn keyword cType Inputs
     1326syn keyword cType IntArrayInput
    13231327syn keyword cType IntInput
    13241328syn keyword cType IntMatParam
    … …  
    13281332syn keyword cType IssmDirectApplicInterface
    13291333syn keyword cType IssmParallelDirectApplicInterface
     1334syn keyword cType krigingobjects
    13301335syn keyword cType Load
    13311336syn keyword cType Loads
    … …  
    13381343syn keyword cType Matice
    13391344syn keyword cType Matlitho
     1345syn keyword cType matrixobjects
    13401346syn keyword cType MatrixParam
    13411347syn keyword cType Misfit
    … …  
    13501356syn keyword cType Observations
    13511357syn keyword cType Option
     1358syn keyword cType Options
    13521359syn keyword cType OptionUtilities
    1353 syn keyword cType Options
    13541360syn keyword cType Param
    13551361syn keyword cType Parameters
    … …  
    13641370syn keyword cType Regionaloutput
    13651371syn keyword cType Results
     1372syn keyword cType Riftfront
    13661373syn keyword cType RiftStruct
    1367 syn keyword cType Riftfront
    13681374syn keyword cType Seg
    13691375syn keyword cType SegInput
     1376syn keyword cType Segment
    13701377syn keyword cType SegRef
    1371 syn keyword cType Segment
    13721378syn keyword cType SpcDynamic
    13731379syn keyword cType SpcStatic
    … …  
    13891395syn keyword cType Vertex
    13901396syn keyword cType Vertices
    1391 syn keyword cType classes
    1392 syn keyword cType gaussobjects
    1393 syn keyword cType krigingobjects
    1394 syn keyword cType matrixobjects
    13951397syn keyword cType AdjointBalancethickness2Analysis
    13961398syn keyword cType AdjointBalancethicknessAnalysis
    … …  
    14111413syn keyword cType FreeSurfaceBaseAnalysis
    14121414syn keyword cType FreeSurfaceTopAnalysis
     1415syn keyword cType GiaIvinsAnalysis
    14131416syn keyword cType GLheightadvectionAnalysis
    1414 syn keyword cType GiaIvinsAnalysis
    14151417syn keyword cType HydrologyDCEfficientAnalysis
    14161418syn keyword cType HydrologyDCInefficientAnalysis

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

    r24088 r24097  
    777777        ThicknessAlongGradientEnum,
    778778        ThicknessEnum,
     779        ThicknessOldEnum,
    779780        ThicknessPositiveEnum,
    780781        VelEnum,

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

    r24088 r24097  
    783783                case ThicknessAlongGradientEnum : return "ThicknessAlongGradient";
    784784                case ThicknessEnum : return "Thickness";
     785                case ThicknessOldEnum : return "ThicknessOld";
    785786                case ThicknessPositiveEnum : return "ThicknessPositive";
    786787                case VelEnum : return "Vel";

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

    r24088 r24097  
    801801              else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
    802802              else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
     803              else if (strcmp(name,"ThicknessOld")==0) return ThicknessOldEnum;
    803804              else if (strcmp(name,"ThicknessPositive")==0) return ThicknessPositiveEnum;
    804805              else if (strcmp(name,"Vel")==0) return VelEnum;
    … …  
    874875              else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum;
    875876              else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
    876               else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum;
     880              if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum;
     881              else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum;
    881882              else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum;
    882883              else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
    … …  
    997998              else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
    998999              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
    999               else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
     1003              if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
     1004              else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
    10041005              else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
    10051006              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
    … …  
    11201121              else if (strcmp(name,"Matestar")==0) return MatestarEnum;
    11211122              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    1122               else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
     1126              if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
     1127              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    11271128              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
    11281129              else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
    … …  
    12431244              else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
    12441245              else if (strcmp(name,"StringParam")==0) return StringParamEnum;
    1245               else if (strcmp(name,"SubelementFriction1")==0) return SubelementFriction1Enum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum;
     1249              if (strcmp(name,"SubelementFriction1")==0) return SubelementFriction1Enum;
     1250              else if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum;
    12501251              else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum;
    12511252              else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum;
Note: See TracChangeset for help on using the changeset viewer.