Changeset 26204


Ignore:
Timestamp:
04/22/21 06:08:36 (4 years ago)
Author:
tsantos
Message:

CHG: working on MLHO, adding new Enums: vxBase, vyBase, vxSurface, vySurface

Location:
issm/trunk-jpl/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r26192 r26204  
    783783        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum,0.);
    784784        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum,0.);
    785         if(isMLHO){//itapopo FIXME the shear velocities should be initialized correctly
     785        if(isMLHO){//itapopo FIXME shear and base velocities should be initialized correctly
    786786                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxShearEnum,0.);
    787787                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyShearEnum,0.);
     788                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum,0.);
     789                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum,0.);
    788790        }
    789791        iomodel->FetchDataToInput(inputs,elements,"md.stressbalance.loadingforcex",LoadingforceXEnum);
    … …  
    32023204        /*Fetch number of nodes and dof for this finite element*/
    32033205        int numnodes = basalelement->GetNumberOfNodes();
    3204         int numdof   = numnodes*2*2; //4 DOFs per node
     3206        int numdof   = numnodes*4; //4 DOFs per node
    32053207
    32063208        /*Fetch dof list and allocate solution vectors*/
    … …  
    32083210        //basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
    32093211        IssmDouble* values    = xNew<IssmDouble>(numdof);
    3210         IssmDouble* vx        = xNew<IssmDouble>(numnodes);
    3211         IssmDouble* vy        = xNew<IssmDouble>(numnodes);
     3212        IssmDouble* vbx       = xNew<IssmDouble>(numnodes);
     3213        IssmDouble* vby       = xNew<IssmDouble>(numnodes);
    32123214   IssmDouble* vshx      = xNew<IssmDouble>(numnodes);
    32133215   IssmDouble* vshy      = xNew<IssmDouble>(numnodes);
     3216        IssmDouble* vsx       = xNew<IssmDouble>(numnodes);
     3217        IssmDouble* vsy       = xNew<IssmDouble>(numnodes);
    32143218        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
    32153219        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
    … …  
    32293233   /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    32303234   for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written
    3231       vx[i] =values[i*4+0]; //base vx
     3235      vbx[i] =values[i*4+0]; //base vx
    32323236      vshx[i]=values[i*4+1]; //shear vx
    3233                 if(xIsNan<IssmDouble>(vx[i]))           _error_("NaN found in solution vector");
    3234       if(xIsInf<IssmDouble>(vx[i]))             _error_("Inf found in solution vector");
     3237                vsx[i] =vbx[i]+vshx[i]; //surface vx
     3238                if(xIsNan<IssmDouble>(vbx[i]))  _error_("NaN found in solution vector");
     3239      if(xIsInf<IssmDouble>(vbx[i]))    _error_("Inf found in solution vector");
    32353240                if(xIsNan<IssmDouble>(vshx[i])) _error_("NaN found in solution vector");
    32363241      if(xIsInf<IssmDouble>(vshx[i]))   _error_("Inf found in solution vector");
    3237                 vy[i] =values[i*4+2]; //base vy
     3242                vby[i] =values[i*4+2]; //base vy
    32383243                vshy[i]=values[i*4+3]; //shear vy
    3239                 if(xIsNan<IssmDouble>(vy[i]))           _error_("NaN found in solution vector");
    3240                 if(xIsInf<IssmDouble>(vy[i]))           _error_("Inf found in solution vector");
     3244                vsy[i] =vby[i]+vshy[i]; //surface vy
     3245                if(xIsNan<IssmDouble>(vby[i]))  _error_("NaN found in solution vector");
     3246                if(xIsInf<IssmDouble>(vby[i]))  _error_("Inf found in solution vector");
    32413247                if(xIsNan<IssmDouble>(vshy[i])) _error_("NaN found in solution vector");
    32423248                if(xIsInf<IssmDouble>(vshy[i])) _error_("Inf found in solution vector");
    … …  
    32493255        /*Get Vz and compute vel (base)*/
    32503256        basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
    3251         for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     3257        for(i=0;i<numnodes;i++) vel[i]=sqrt(vbx[i]*vbx[i] + vby[i]*vby[i] + vz[i]*vz[i]);
    32523258
    32533259        /*Add vx and vy as inputs to the tria element (base velocities): */
    3254         element->AddBasalInput(VxEnum,vx,element->GetElementType());
    3255         element->AddBasalInput(VyEnum,vy,element->GetElementType());
    3256         element->AddBasalInput(VelEnum,vel,element->GetElementType());
     3260        element->AddBasalInput(VxBaseEnum,vbx,element->GetElementType());
     3261        element->AddBasalInput(VyBaseEnum,vby,element->GetElementType());
     3262        element->AddBasalInput(VelEnum,vel,element->GetElementType()); //itapopo check this
     3263       
     3264        /*Add vx and vy as inputs to the tria element (surface velocities): */
     3265        element->AddBasalInput(VxSurfaceEnum,vsx,element->GetElementType());
     3266        element->AddBasalInput(VySurfaceEnum,vsy,element->GetElementType());
    32573267
    32583268        /*Free ressources:*/
    32593269        xDelete<IssmDouble>(vel);
    32603270        xDelete<IssmDouble>(vz);
    3261         xDelete<IssmDouble>(vy);
    3262         xDelete<IssmDouble>(vx);
     3271        xDelete<IssmDouble>(vby);
     3272        xDelete<IssmDouble>(vbx);
     3273        xDelete<IssmDouble>(vsy);
     3274        xDelete<IssmDouble>(vsx);
    32633275        xDelete<IssmDouble>(vshy);
    32643276        xDelete<IssmDouble>(vshx);
    … …  
    32703282void           StressbalanceAnalysis::GetSolutionFromInputsMLHO(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    32713283
    3272         IssmDouble   vx,vy,vshx,vshy;
     3284        IssmDouble   vbx,vby,vshx,vshy;
    32733285        int          domaintype,dim,approximation,dofpernode;
    32743286        int*         doflist = NULL;
    … …  
    32933305
    32943306        /*Get inputs*/
    3295         Input* vx_input         =element->GetInput(VxEnum);             _assert_(vx_input);
     3307        Input* vxbase_input     =element->GetInput(VxBaseEnum); _assert_(vxbase_input);
    32963308        Input* vxshear_input    =element->GetInput(VxShearEnum);        _assert_(vxshear_input);
    3297         Input* vy_input         =element->GetInput(VyEnum);             _assert_(vy_input);
     3309        Input* vybase_input     =element->GetInput(VyBaseEnum); _assert_(vybase_input);
    32983310        Input* vyshear_input    =element->GetInput(VyShearEnum);        _assert_(vyshear_input);
    32993311
    … …  
    33043316
    33053317                /*Recover vx and vy*/
    3306                 vx_input->GetInputValue(&vx,gauss);                     //base vx
     3318                vxbase_input->GetInputValue(&vbx,gauss);        //base vx
    33073319                vxshear_input->GetInputValue(&vshx,gauss);//shear vx
    3308                 values[i*4+0]=vx;   //base vx
     3320                values[i*4+0]=vbx;  //base vx
    33093321                values[i*4+1]=vshx; //shear vx
    3310                 vy_input->GetInputValue(&vy,gauss);                     //base vy
     3322                vybase_input->GetInputValue(&vby,gauss);        //base vy
    33113323                vyshear_input->GetInputValue(&vshy,gauss);//shear vy
    3312                 values[i*4+2]=vy;               //base vy
     3324                values[i*4+2]=vby; //base vy
    33133325                values[i*4+3]=vshy;//shear vy 
    33143326        }

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

    r26175 r26204  
    852852}/*}}}*/
    853853IssmDouble Friction::VelMag(Gauss* gauss){/*{{{*/
    854         /*Get effective pressure as a function of  flag */
    855 
     854        /*Get the velocity magnitude as a function of flag */
    856855
    857856        /*diverse*/
     857        bool isMLHO=false; //employing mono layer higher order model? Default is no
    858858        IssmDouble vx,vy,vz,vmag;
     859       
     860        element->parameters->FindParam(&isMLHO,FlowequationIsMLHOEnum);
     861        if(isMLHO) _assert_(dim==2);
    859862
    860863        switch(dim){
    … …  
    864867                        break;
    865868                case 2:
    866                         element->GetInputValue(&vx,gauss,VxEnum);
    867                         element->GetInputValue(&vy,gauss,VyEnum);
     869                        if(isMLHO){
     870                                element->GetInputValue(&vx,gauss,VxBaseEnum);
     871                                element->GetInputValue(&vy,gauss,VyBaseEnum);
     872                        }
     873                        else{
     874                                element->GetInputValue(&vx,gauss,VxEnum);
     875                                element->GetInputValue(&vy,gauss,VyEnum);
     876                        }               
    868877                        vmag=sqrt(vx*vx+vy*vy);
    869878                        break;

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

    r26196 r26204  
    936936syn keyword cConstant VelEnum
    937937syn keyword cConstant VxAverageEnum
     938syn keyword cConstant VxBaseEnum
    938939syn keyword cConstant VxEnum
    939940syn keyword cConstant VxMeshEnum
    940941syn keyword cConstant VxObsEnum
    941942syn keyword cConstant VxShearEnum
     943syn keyword cConstant VxSurfaceEnum
    942944syn keyword cConstant VyAverageEnum
     945syn keyword cConstant VyBaseEnum
    943946syn keyword cConstant VyEnum
    944947syn keyword cConstant VyMeshEnum
    945948syn keyword cConstant VyObsEnum
    946949syn keyword cConstant VyShearEnum
     950syn keyword cConstant VySurfaceEnum
    947951syn keyword cConstant VzEnum
    948952syn keyword cConstant VzFSEnum

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

    r26196 r26204  
    933933        VelEnum,
    934934        VxAverageEnum,
     935        VxBaseEnum,
    935936        VxEnum,
    936937        VxMeshEnum,
    937938        VxObsEnum,
    938939        VxShearEnum,
     940        VxSurfaceEnum,
    939941        VyAverageEnum,
     942        VyBaseEnum,
    940943        VyEnum,
    941944        VyMeshEnum,
    942945        VyObsEnum,
    943946        VyShearEnum,
     947        VySurfaceEnum,
    944948        VzEnum,
    945949        VzFSEnum,

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

    r26196 r26204  
    938938                case VelEnum : return "Vel";
    939939                case VxAverageEnum : return "VxAverage";
     940                case VxBaseEnum : return "VxBase";
    940941                case VxEnum : return "Vx";
    941942                case VxMeshEnum : return "VxMesh";
    942943                case VxObsEnum : return "VxObs";
    943944                case VxShearEnum : return "VxShear";
     945                case VxSurfaceEnum : return "VxSurface";
    944946                case VyAverageEnum : return "VyAverage";
     947                case VyBaseEnum : return "VyBase";
    945948                case VyEnum : return "Vy";
    946949                case VyMeshEnum : return "VyMesh";
    947950                case VyObsEnum : return "VyObs";
    948951                case VyShearEnum : return "VyShear";
     952                case VySurfaceEnum : return "VySurface";
    949953                case VzEnum : return "Vz";
    950954                case VzFSEnum : return "VzFS";

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

    r26196 r26204  
    959959              else if (strcmp(name,"Vel")==0) return VelEnum;
    960960              else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
     961              else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
    961962              else if (strcmp(name,"Vx")==0) return VxEnum;
    962963              else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
    963964              else if (strcmp(name,"VxObs")==0) return VxObsEnum;
    964965              else if (strcmp(name,"VxShear")==0) return VxShearEnum;
     966              else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum;
    965967              else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
     968              else if (strcmp(name,"VyBase")==0) return VyBaseEnum;
    966969              else if (strcmp(name,"Vy")==0) return VyEnum;
    967970              else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
    968971              else if (strcmp(name,"VyObs")==0) return VyObsEnum;
    969972              else if (strcmp(name,"VyShear")==0) return VyShearEnum;
     973              else if (strcmp(name,"VySurface")==0) return VySurfaceEnum;
    970974              else if (strcmp(name,"Vz")==0) return VzEnum;
    971975              else if (strcmp(name,"VzFS")==0) return VzFSEnum;
    … …  
    994998              else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
    995999              else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
    996               else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
    9971004              else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
    9981005              else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
    9991006              else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
     1007              else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
    10041008              else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
    10051009              else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
    … …  
    11171121              else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    11181122              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    1119               else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    11201127              else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
    11211128              else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
    11221129              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
     1130              else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
    11271131              else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
    11281132              else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
    … …  
    12401244              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    12411245              else if (strcmp(name,"ElementInput")==0) return ElementInputEnum;
    1242               else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
    12431250              else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
    12441251              else if (strcmp(name,"IntParam")==0) return IntParamEnum;
    12451252              else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"Inputs")==0) return InputsEnum;
     1253              else if (strcmp(name,"Inputs")==0) return InputsEnum;
    12501254              else if (strcmp(name,"Internal")==0) return InternalEnum;
    12511255              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
    … …  
    13631367              else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum;
    13641368              else if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum;
    1365               else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
    13661373              else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
    13671374              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    13681375              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
     1376              else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
    13731377              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
    13741378              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;

  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m

    r26192 r26204  
    228228                field = field*yts;
    229229        elseif strcmp(fieldname,'VyShear'),
     230                field = field*yts;
     231        elseif strcmp(fieldname,'VxBase'),
     232                field = field*yts;
     233        elseif strcmp(fieldname,'VyBase'),
     234                field = field*yts;
     235        elseif strcmp(fieldname,'VxSurface'),
     236                field = field*yts;
     237        elseif strcmp(fieldname,'VySurface'),
    230238                field = field*yts;
    231239        elseif strcmp(fieldname,'BasalforcingsGroundediceMeltingRate'),

  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py

    r26193 r26204  
    179179            field = field * yts
    180180        elif fieldname == 'VyShear':
     181            field = field * yts
     182        elif fieldname == 'VxBase':
     183            field = field * yts
     184        elif fieldname == 'VyBase':
     185            field = field * yts
     186        elif fieldname == 'VxSurface':
     187            field = field * yts
     188        elif fieldname == 'VySurface':
    181189            field = field * yts
    182190        elif fieldname == 'BasalforcingsGroundediceMeltingRate':
Note: See TracChangeset for help on using the changeset viewer.