Changeset 27138


Ignore:
Timestamp:
07/01/22 17:35:05 (3 years ago)
Author:
Eric.Larour
Message:

CHG: ice and ocean levelsets were not being processed well, as L2projection was hijacking, and
mme mass transport analysis was not updating it well.

Location:
issm/branches/trunk-larour-SLPS2022/src/c
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • issm/branches/trunk-larour-SLPS2022/src/c/analyses/MmemasstransportAnalysis.cpp

    r27126 r27138  
    4040        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    4141        iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0);
    42         iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    43         iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
     42        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MmemasstransportMaskIceLevelsetEnum);
     43        iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MmemasstransportMaskOceanLevelsetEnum);
    4444        iomodel->FetchDataToInput(inputs,elements,"md.mmemasstransport.thickness", MmemasstransportThicknessEnum);
    4545
    … …  
    9999void           MmemasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    100100
    101         /*retrieve thickness from MmemasstransportAnalysis spcs in our element:*/
     101        /*retrieve thickness, ice and ocean levelsets from MmemasstransportAnalysis inputs in our element:*/
    102102
    103         IssmDouble h;
     103        IssmDouble height,ice,ocean;
    104104        int       *doflist = NULL;
    105105
    106106        /*Fetch number of nodes and initialize values*/
    107107        int         numnodes = element->GetNumberOfNodes();
    108         int         numdof   = numnodes;
     108        int         numdof   = numnodes*3;
    109109        IssmDouble* values   = xNew<IssmDouble>(numdof);
    110110
    … …  
    112112        element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    113113        Input* h_input=element->GetInput(MmemasstransportThicknessEnum); _assert_(h_input);
     114        Input* i_input=element->GetInput(MmemasstransportMaskIceLevelsetEnum); _assert_(i_input);
     115        Input* o_input=element->GetInput(MmemasstransportMaskOceanLevelsetEnum); _assert_(o_input);
    114116
    115117        /*Ok, we have the thickness in inputs, fill in solution */
    … …  
    117119        for(int i=0;i<numnodes;i++){
    118120                gauss->GaussVertex(i);
    119                 h_input->GetInputValue(&h,gauss);
    120                 values[i+0]=h;
     121                h_input->GetInputValue(&height,gauss);
     122                i_input->GetInputValue(&ice,gauss);
     123                o_input->GetInputValue(&ocean,gauss);
     124                values[3*i+0]=height;
     125                values[3*i+1]=ice;
     126                values[3*i+2]=ocean;
    121127        }
    122128
    … …  
    139145        /*Fetch number of nodes and dof for this finite element*/
    140146        int numnodes = element->GetNumberOfNodes();
    141         int numdof   = numnodes;
     147        int numdof   = numnodes*3;
    142148
    143149        /*Fetch dof list and allocate solution vectors*/
    144150        element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
    145151        IssmDouble* values    = xNew<IssmDouble>(numdof);
    146         IssmDouble* h        = xNew<IssmDouble>(numnodes);
     152        IssmDouble* height        = xNew<IssmDouble>(numnodes);
     153        IssmDouble* ice        = xNew<IssmDouble>(numnodes);
     154        IssmDouble* ocean        = xNew<IssmDouble>(numnodes);
    147155
    148156        /*Use the dof list to index into the solution vector: */
    … …  
    151159        /*Retrieve h:*/
    152160        for(i=0;i<numnodes;i++){
    153                 h[i]=values[i+0];
     161                height[i]=values[3*i+0];
     162                ice[i]=values[3*i+1];
     163                ocean[i]=values[3*i+2];
    154164
    155165                /*Check solution*/
    156                 if(xIsNan<IssmDouble>(h[i])) _error_("NaN found in ice thickness solution vector");
    157                 if(xIsInf<IssmDouble>(h[i])) _error_("Inf found in ice thickness solution vector");
     166                if(xIsNan<IssmDouble>(height[i])) _error_("NaN found in ice thickness solution vector");
     167                if(xIsInf<IssmDouble>(height[i])) _error_("Inf found in ice thickness solution vector");
     168                if(xIsNan<IssmDouble>(ice[i])) _error_("NaN found in ice mask solution vector");
     169                if(xIsInf<IssmDouble>(ice[i])) _error_("Inf found in ice mask solution vector");
     170                if(xIsNan<IssmDouble>(ocean[i])) _error_("NaN found in ocean mask solution vector");
     171                if(xIsInf<IssmDouble>(ocean[i])) _error_("Inf found in ocean mask solution vector");
     172
    158173        }
    159174
    160         /*Add thickness as inputs to the tria element: */
    161         element->AddInput(ThicknessEnum,h,P1Enum);
     175        /*Add thickness ice and ocean levelsets as inputs to the tria element: */
     176        element->AddInput(ThicknessEnum,height,P1Enum);
     177        element->AddInput(MaskIceLevelsetEnum,ice,P1Enum);
     178        element->AddInput(MaskOceanLevelsetEnum,ocean,P1Enum);
    162179
    163180        /*Free ressources:*/
    164         xDelete<IssmDouble>(h);
     181        xDelete<IssmDouble>(height);
     182        xDelete<IssmDouble>(ice);
     183        xDelete<IssmDouble>(ocean);
    165184        xDelete<IssmDouble>(values);
    166185        xDelete<int>(doflist);

  • issm/branches/trunk-larour-SLPS2022/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp

    r27095 r27138  
    173173                        MmeToInputx(femmodel,distributed_values,variable_partition,npart,MmemasstransportThicknessEnum, P0Enum);
    174174
    175                 if (femmodel->inputs->GetInputObjectEnum(MaskIceLevelsetEnum)==DatasetInputEnum)
    176                         MmeToInputx(femmodel,distributed_values,variable_partition,npart,MaskIceLevelsetEnum, P1Enum);
    177 
    178                 if (femmodel->inputs->GetInputObjectEnum(MaskOceanLevelsetEnum)==DatasetInputEnum)
    179                         MmeToInputx(femmodel,distributed_values,variable_partition,npart,MaskOceanLevelsetEnum, P1Enum);
     175                if (femmodel->inputs->GetInputObjectEnum(MmemasstransportMaskIceLevelsetEnum)==DatasetInputEnum){
     176                        MmeToInputx(femmodel,distributed_values,variable_partition,npart,MmemasstransportMaskIceLevelsetEnum, P1Enum);
     177                        /*delete MaskIceLevelsetEnum which will be replaced by MmemasstransportMaskIceLevelsetEnum for each time step: */
     178                        femmodel->inputs->DeleteInput(MaskIceLevelsetEnum);
     179                }
     180
     181                if (femmodel->inputs->GetInputObjectEnum(MmemasstransportMaskOceanLevelsetEnum)==DatasetInputEnum){
     182                        MmeToInputx(femmodel,distributed_values,variable_partition,npart,MmemasstransportMaskOceanLevelsetEnum, P1Enum);
     183                       
     184                        /*delete MaskOceanLevelsetEnum which will be replaced by MmemasstransportMaskOceanLevelsetEnum for each time step: */
     185                        femmodel->inputs->DeleteInput(MaskOceanLevelsetEnum);
     186                }
    180187
    181188        } /*}}}*/

  • issm/branches/trunk-larour-SLPS2022/src/c/modules/UpdateMmesx/UpdateMmesx.cpp

    r27095 r27138  
    2525
    2626        } /*}}}*/
    27         /*Deal with ocean: {{{*/
     27        /*Deal with dsl: {{{*/
    2828        if (femmodel->inputs->Exist(OceantransportSpcbottompressureEnum) && femmodel->inputs->GetInputObjectEnum(OceantransportSpcbottompressureEnum)==DatasetInputEnum){
    2929               
    … …  
    5151                MmeToInputx(femmodel,modelids,partition,nids,MmemasstransportThicknessEnum, P0Enum);
    5252
    53                 if (femmodel->inputs->GetInputObjectEnum(MaskIceLevelsetEnum)==DatasetInputEnum)
    54                         MmeToInputx(femmodel,modelids,partition,nids,MaskIceLevelsetEnum, P1Enum);
     53                if (femmodel->inputs->GetInputObjectEnum(MmemasstransportMaskIceLevelsetEnum)==DatasetInputEnum){
     54                        MmeToInputx(femmodel,modelids,partition,nids,MmemasstransportMaskIceLevelsetEnum, P1Enum);
     55                        /*delete MaskIceLevelsetEnum which will be replaced by MmemasstransportMaskIceLevelsetEnum for each time step: */
     56                        femmodel->inputs->DeleteInput(MaskIceLevelsetEnum);
     57                }
    5558
    56                 if (femmodel->inputs->GetInputObjectEnum(MaskOceanLevelsetEnum)==DatasetInputEnum)
    57                         MmeToInputx(femmodel,modelids,partition,nids,MaskOceanLevelsetEnum, P1Enum);
     59                if (femmodel->inputs->GetInputObjectEnum(MmemasstransportMaskOceanLevelsetEnum)==DatasetInputEnum){
     60                        MmeToInputx(femmodel,modelids,partition,nids,MmemasstransportMaskOceanLevelsetEnum, P1Enum);
     61                        /*delete MaskOceanLevelsetEnum which will be replaced by MmemasstransportMaskOceanLevelsetEnum for each time step: */
     62                        femmodel->inputs->DeleteInput(MaskOceanLevelsetEnum);
     63                }
    5864
    5965                /*free ressources:*/

  • issm/branches/trunk-larour-SLPS2022/src/c/shared/Enum/Enum.vim

    r27094 r27138  
    811811syn keyword cConstant MasstransportSpcthicknessEnum
    812812syn keyword cConstant MmemasstransportThicknessEnum
     813syn keyword cConstant MmemasstransportMaskIceLevelsetEnum
     814syn keyword cConstant MmemasstransportMaskOceanLevelsetEnum
    813815syn keyword cConstant MaterialsRheologyBEnum
    814816syn keyword cConstant MaterialsRheologyBbarEnum

  • issm/branches/trunk-larour-SLPS2022/src/c/shared/Enum/EnumDefinitions.h

    r27094 r27138  
    807807        MasstransportSpcthicknessEnum,
    808808        MmemasstransportThicknessEnum,
     809        MmemasstransportMaskIceLevelsetEnum,
     810        MmemasstransportMaskOceanLevelsetEnum,
    809811        MaterialsRheologyBEnum,
    810812        MaterialsRheologyBbarEnum,

  • issm/branches/trunk-larour-SLPS2022/src/c/shared/Enum/EnumToStringx.cpp

    r27094 r27138  
    813813                case MasstransportSpcthicknessEnum : return "MasstransportSpcthickness";
    814814                case MmemasstransportThicknessEnum : return "MmemasstransportThickness";
     815                case MmemasstransportMaskIceLevelsetEnum : return "MmemasstransportMaskIceLevelset";
     816                case MmemasstransportMaskOceanLevelsetEnum : return "MmemasstransportMaskOceanLevelset";
    815817                case MaterialsRheologyBEnum : return "MaterialsRheologyB";
    816818                case MaterialsRheologyBbarEnum : return "MaterialsRheologyBbar";

  • issm/branches/trunk-larour-SLPS2022/src/c/shared/Enum/Enumjl.vim

    r27094 r27138  
    804804syn keyword juliaConstC MasstransportSpcthicknessEnum
    805805syn keyword juliaConstC MmemasstransportThicknessEnum
     806syn keyword juliaConstC MmemasstransportMaskIceLevelsetEnum
     807syn keyword juliaConstC MmemasstransportMaskOceanLevelsetEnum
    806808syn keyword juliaConstC MaterialsRheologyBEnum
    807809syn keyword juliaConstC MaterialsRheologyBbarEnum

  • issm/branches/trunk-larour-SLPS2022/src/c/shared/Enum/StringToEnumx.cpp

    r27094 r27138  
    831831              else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum;
    832832              else if (strcmp(name,"MmemasstransportThickness")==0) return MmemasstransportThicknessEnum;
     833              else if (strcmp(name,"MmemasstransportMaskIceLevelset")==0) return MmemasstransportMaskIceLevelsetEnum;
     834              else if (strcmp(name,"MmemasstransportMaskOceanLevelset")==0) return MmemasstransportMaskOceanLevelsetEnum;
    833835              else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum;
    834836              else if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum;
    … …  
    873875              else if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum;
    874876              else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum;
    875               else if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum;
    876               else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Sealevel")==0) return SealevelEnum;
     880              if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum;
     881              else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum;
     882              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
    881883              else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
    882884              else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum;
    … …  
    996998              else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
    997999              else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
    998               else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
    999               else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
     1003              if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
     1004              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
     1005              else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
    10041006              else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
    10051007              else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum;
    … …  
    11191121              else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum;
    11201122              else if (strcmp(name,"Vel")==0) return VelEnum;
    1121               else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    1122               else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Vx")==0) return VxEnum;
     1126              if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
     1127              else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
     1128              else if (strcmp(name,"Vx")==0) return VxEnum;
    11271129              else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
    11281130              else if (strcmp(name,"VxObs")==0) return VxObsEnum;
    … …  
    12421244              else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum;
    12431245              else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum;
    1244               else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum;
    1245               else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
     1249              if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum;
     1250              else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
     1251              else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
    12501252              else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
    12511253              else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
    … …  
    13651367              else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum;
    13661368              else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum;
    1367               else if (strcmp(name,"FrontalForcingsRignotAutoregression")==0) return FrontalForcingsRignotAutoregressionEnum;
    1368               else if (strcmp(name,"Fset")==0) return FsetEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
     1372              if (strcmp(name,"FrontalForcingsRignotAutoregression")==0) return FrontalForcingsRignotAutoregressionEnum;
     1373              else if (strcmp(name,"Fset")==0) return FsetEnum;
     1374              else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
    13731375              else if (strcmp(name,"GLheightadvectionAnalysis")==0) return GLheightadvectionAnalysisEnum;
    13741376              else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
    … …  
    14881490              else if (strcmp(name,"MinVy")==0) return MinVyEnum;
    14891491              else if (strcmp(name,"MinVz")==0) return MinVzEnum;
    1490               else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
    1491               else if (strcmp(name,"Moulin")==0) return MoulinEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
     1495              if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
     1496              else if (strcmp(name,"Moulin")==0) return MoulinEnum;
     1497              else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
    14961498              else if (strcmp(name,"Mpi")==0) return MpiEnum;
    14971499              else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
    … …  
    16111613              else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
    16121614              else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
    1613               else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
    1614               else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;
    16151615         else stage=14;
    16161616   }
    16171617   if(stage==14){
    1618               if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
     1618              if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
     1619              else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;
     1620              else if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
    16191621              else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum;
    16201622              else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum;

  • issm/branches/trunk-larour-SLPS2022/src/c/shared/Enum/issmenums.jl

    r27094 r27138  
    800800        MasstransportSpcthicknessEnum
    801801        MmemasstransportThicknessEnum
     802        MmemasstransportMaskIceLevelsetEnum
     803        MmemasstransportMaskOceanLevelsetEnum
    802804        MaterialsRheologyBEnum
    803805        MaterialsRheologyBbarEnum
    … …  
    24062408        if(enum==MasstransportSpcthicknessEnum) return "MasstransportSpcthickness" end
    24072409        if(enum==MmemasstransportThicknessEnum) return "MmemasstransportThickness" end
     2410        if(enum==MmemasstransportMaskIceLevelsetEnum) return "MmemasstransportMaskIceLevelset" end
     2411        if(enum==MmemasstransportMaskOceanLevelsetEnum) return "MmemasstransportMaskOceanLevelset" end
    24082412        if(enum==MaterialsRheologyBEnum) return "MaterialsRheologyB" end
    24092413        if(enum==MaterialsRheologyBbarEnum) return "MaterialsRheologyBbar" end
Note: See TracChangeset for help on using the changeset viewer.