Changeset 26292


Ignore:
Timestamp:
05/26/21 17:29:06 (4 years ago)
Author:
tsantos
Message:

CHG: added viscosity computation for Mono Layer Higher Order model, and added penalty methods for L1L2 and Mono Layer model

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

Legend:

Unmodified
Added
Removed

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

    r26290 r26292  
    2121//#define LATERALFRICTION 1
    2222//#define DISCSLOPE 1 //testing for SSA
     23#define NOSPCSHEARVEL 1 //MLHO
    2324
    2425/*Model processing*/
    … …  
    194195                                }
    195196                        }
    196                         else{//MLHO itapopo FIXME base and shear velocities are constrained
     197                        else{//MLHO
     198                                #ifdef NOSPCSHEARVEL   
     199                                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0);
     200                                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2);
     201                                #else
     202                                /*Default: apply spcs to shear vx and shear vy*/
    197203                                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0);
    198204                                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1);
    199205                                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2);
    200206                                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3);
     207                                #endif
    201208                        }
    202209                }
    … …  
    781788        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum,0.);
    782789        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum,0.);
    783         if(isMLHO){//itapopo FIXME shear and base velocities should be initialized correctly
     790        if(isMLHO){
     791                /*itapopo FIXME applying the same initialization for shear vx and shear vy for now*/
    784792                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxShearEnum,0.);
    785793                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyShearEnum,0.);
    … …  
    29072915        Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input);
    29082916        Input* surface_input  = element->GetInput(SurfaceEnum);   _assert_(surface_input);
    2909         Input* vx_input       = element->GetInput(VxEnum);        _assert_(vx_input);
    2910         Input* vy_input       = element->GetInput(VyEnum);        _assert_(vy_input);
     2917        Input* vx_input       = element->GetInput(VxEnum);        _assert_(vx_input); //vertically integrated vx
     2918        Input* vy_input       = element->GetInput(VyEnum);        _assert_(vy_input); //vertically integrated vy
     2919        Input* vxshear_input  = element->GetInput(VxShearEnum);   _assert_(vxshear_input);
     2920        Input* vyshear_input  = element->GetInput(VyShearEnum);   _assert_(vyshear_input);
    29112921        Input* n_input                   = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
    29122922
    … …  
    29232933                thickness_input->GetInputValue(&thickness, gauss);
    29242934                n_input->GetInputValue(&n,gauss);
    2925                 //FIXME testing with L1L2-type viscosity
    2926                 element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
     2935                IssmDouble dim=2;//itapopo
     2936                element->material->ViscosityMLHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vxshear_input,vyshear_input,thickness_input);
    29272937                //viscosity=10e13;//itapopo
    29282938

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

    r26271 r26292  
    42124212
    42134213}/*}}}*/
     4214void       Element::StrainRateMLHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input){/*{{{*/
     4215        /*Compute the 2d Blatter/MLHO Strain Rate (5 components):
     4216         *
     4217         * epsilon=[exx eyy exy exz eyz]
     4218         *
     4219         * with exz=1/2 du/dz
     4220         *      eyz=1/2 dv/dz
     4221         *
     4222         * the contribution of vz is neglected
     4223         */
     4224
     4225        /*Intermediaries*/
     4226        IssmDouble dvx[2]; // [dvx/dx, dvx/dy]
     4227        IssmDouble dvy[2]; // [dvy/dx, dvy/dy]
     4228        IssmDouble vxshear,vyshear,thickness;
     4229
     4230        /*Check that both inputs have been found*/
     4231        if (!vx_input || !vy_input || !vxshear_input || !vyshear_input || !thickness_input){
     4232                _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vxshear: " << vxshear_input << ", vyshear: " << vyshear_input << ", thickness: " << thickness_input << "\n");
     4233        }
     4234
     4235        /*Get strain rate assuming that epsilon has been allocated*/
     4236        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     4237        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     4238        thickness_input->GetInputValue(&thickness,gauss); _assert_(thickness>0);
     4239   vxshear_input->GetInputValue(&vxshear,gauss);
     4240   vyshear_input->GetInputValue(&vyshear,gauss);
     4241       
     4242        epsilon[0] = dvx[0];
     4243        epsilon[1] = dvy[1];
     4244        epsilon[2] = 0.5*(dvx[1] + dvy[0]);
     4245        epsilon[3] = 0.5*vxshear/thickness;
     4246        epsilon[4] = 0.5*vyshear/thickness;
     4247
     4248}/*}}}*/
    42144249void       Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
    42154250

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r26272 r26292  
    181181                void               StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    182182                void               StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     183                void               StrainRateMLHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input);
    183184                void               StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    184185                void               StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input);

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

    r26144 r26292  
    355355                                default: _error_("Approximation "<<EnumToStringx(approximation1)<<" not supported yet");
    356356                        }
     357                case MLHOApproximationEnum:
     358                        switch(approximation1){
     359                                case MLHOApproximationEnum:   return PenaltyCreateKMatrixStressbalanceSSAHO(kmax);
     360                                default: _error_("Approximation "<<EnumToStringx(approximation1)<<" not supported yet");
     361                        }
     362                case L1L2ApproximationEnum:
     363                        switch(approximation1){
     364                                case L1L2ApproximationEnum:   return PenaltyCreateKMatrixStressbalanceSSAHO(kmax);
     365                                default: _error_("Approximation "<<EnumToStringx(approximation1)<<" not supported yet");
     366                        }
    357367                case FSvelocityEnum:
    358368                        switch(approximation1){

  • issm/trunk-jpl/src/c/classes/Materials/Material.h

    r25379 r26292  
    5151                virtual void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
    5252                virtual void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
     53                virtual void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input)=0;
    5354                virtual void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0;
    5455                virtual void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;

  • issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp

    r25508 r26292  
    567567        *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged,gauss);
    568568}/*}}}*/
     569void  Matestar::ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input){/*{{{*/
     570        _error_("not implemented yet");
     571}/*}}}*/
    569572void  Matestar::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
    570573        _error_("not implemented yet");

  • issm/trunk-jpl/src/c/classes/Materials/Matestar.h

    r25911 r26292  
    7676                void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    7777                void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     78                void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input);
    7879                void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
    7980                void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r25911 r26292  
    736736        *pviscosity=viscosity;
    737737}/*}}}*/
     738void  Matice::ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input){/*{{{*/
     739
     740        /*Intermediaries*/
     741        IssmDouble viscosity;
     742        IssmDouble thickness;
     743        IssmDouble epsilon[5];  /* epsilon=[exx,eyy,exy,exz,eyz]; */
     744        IssmDouble eps_eff;
     745
     746        /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
     747        element->StrainRateMLHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vxshear_input,vyshear_input,thickness_input);
     748        eps_eff = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] +  epsilon[3]*epsilon[3] + epsilon[4]*epsilon[4] + epsilon[0]*epsilon[1]);
     749
     750        /*Get viscosity*/
     751        this->GetViscosityBar(&viscosity,eps_eff,gauss);
     752
     753        /*Assign output pointer*/
     754   *pviscosity=viscosity;
     755}/*}}}*/
    738756void  Matice::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
    739757        /*Compute the L1L2 viscosity

  • issm/trunk-jpl/src/c/classes/Materials/Matice.h

    r25911 r26292  
    7878                void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    7979                void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     80                void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input);
    8081                void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
    8182                void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);

  • issm/trunk-jpl/src/c/classes/Materials/Matlitho.h

    r26234 r26292  
    7272                void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
    7373                void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
     74                void       ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input){_error_("not supported");};
    7475                void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not supported");};
    7576                void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
Note: See TracChangeset for help on using the changeset viewer.