Changeset 26850


Ignore:
Timestamp:
02/02/22 08:24:54 (3 years ago)
Author:
Cheng Gong
Message:

NEW: added a fake calving law using linear and nonlinear parameterization

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

Legend:

Unmodified
Added
Removed

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

    r26784 r26850  
    118118                case CalvingTestEnum:
    119119                        break;
     120                case CalvingParameterizationEnum:
     121                        iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum);
     122                        iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum);
     123                        iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
     124                        break;
    120125                default:
    121126                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
    … …  
    172177                case CalvingTestEnum:
    173178                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.speedfactor",CalvingTestSpeedfactorEnum));
     179                        break;
     180                case CalvingParameterizationEnum:
     181                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum));
     182                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.use_param",CalvingUseParamEnum));
     183                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.scale_theta",CalvingScaleThetaEnum));
     184                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.amp_alpha",CalvingAmpAlphaEnum));
     185                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.midp",CalvingMidpointEnum));
     186                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.nonlinearlaw",CalvingNonlinearLawEnum));
    174187                        break;
    175188                default:
    … …  
    572585                }
    573586        }
     587        else if(calvinglaw==CalvingParameterizationEnum){
     588
     589      /*Intermediaries*/
     590      IssmDouble thickness,bed,sealevel;
     591      IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum);
     592
     593                /*Loop over all elements of this partition*/
     594                for(Object* & object : femmodel->elements->objects){
     595                        Element* element  = xDynamicCast<Element*>(object);
     596
     597                        int      numnodes = element->GetNumberOfNodes();
     598                        Gauss*   gauss    = element->NewGauss();
     599                        Input*   H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
     600                        Input*   b_input = element->GetInput(BedEnum); _assert_(b_input);
     601                        Input*   sl_input = element->GetInput(SealevelEnum); _assert_(sl_input);
     602
     603                        /*Potentially constrain nodes of this element*/
     604                        for(int in=0;in<numnodes;in++){
     605                                gauss->GaussNode(element->GetElementType(),in);
     606                                Node* node=element->GetNode(in);
     607                                if(!node->IsActive()) continue;
     608
     609                                H_input->GetInputValue(&thickness,gauss);
     610                                b_input->GetInputValue(&bed,gauss);
     611                                sl_input->GetInputValue(&sealevel,gauss);
     612                                if(thickness<min_thickness && bed<sealevel){
     613                                        node->ApplyConstraint(0,+1.);
     614                                }
     615                                else {
     616                                        /* no ice, set no spc */
     617                                        node->DofInFSet(0);
     618                                }
     619                        }
     620                        delete gauss;
     621                }
     622        }
    574623   else if(calvinglaw==CalvingHabEnum){
    575624

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

    r26836 r26850  
    34413441                                                                                                          case CalvingCrevasseDepthEnum:
    34423442                                                                                                                  this->CalvingCrevasseDepth();
     3443                                                                                                                  break;
     3444                                                                                                          case CalvingParameterizationEnum:
     3445                                                                                                                  this->CalvingRateParameterization();
    34433446                                                                                                                  break;
    34443447                                                                                                          default:

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

    r26836 r26850  
    229229                virtual void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
    230230                virtual void             BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");};
     231                virtual void       CalvingRateParameterization(void){_error_("not implemented yet");};
    231232                virtual void       CalvingRateVonmises(void){_error_("not implemented yet");};
    232233                virtual void       CalvingRateTest(void){_error_("not implemented yet");};

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

    r26830 r26850  
    837837                delete gauss;
    838838        }
     839}
     840/*}}}*/
     841void       Tria::CalvingRateParameterization(){/*{{{*/
     842
     843        IssmDouble  xyz_list[NUMVERTICES][3];
     844        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     845        IssmDouble  calvingratex[NUMVERTICES];
     846        IssmDouble  calvingratey[NUMVERTICES];
     847        IssmDouble  calvingrate[NUMVERTICES];
     848        IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
     849        IssmDouble  sigma_vm[NUMVERTICES];
     850        IssmDouble  B,sigma_max,sigma_max_floating,sigma_max_grounded,n;
     851        IssmDouble  epse_2,groundedice,bed,sealevel;
     852        IssmDouble  mrate, rho_ice, rho_water, thickness, paramX, Hab;
     853        int                     use_parameter=0;
     854        int                     nonlinear_law=0;
     855        IssmDouble  theta, alpha, midp, gamma;
     856
     857        /* Get node coordinates and dof list: */
     858        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     859
     860        /*Retrieve all inputs and parameters we will need*/
     861        Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input);
     862        Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input);
     863        Input* B_input  = this->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
     864        Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
     865        Input* bs_input = this->GetInput(BedEnum);                    _assert_(bs_input);
     866        Input* H_input  = this->GetInput(ThicknessEnum); _assert_(H_input);
     867        Input* smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
     868        Input* smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
     869        Input* n_input  = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
     870        Input* sl_input  = this->GetInput(SealevelEnum); _assert_(sl_input);
     871        Input* mrate_input  = this->GetInput(CalvingMeltingrateEnum); _assert_(mrate_input);
     872
     873        /* Ice and sea water density */
     874        this->FindParam(&rho_ice,MaterialsRhoIceEnum);
     875        this->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
     876
     877        /* Use which parameter  */
     878        this->FindParam(&use_parameter, CalvingUseParamEnum);
     879        this->FindParam(&theta, CalvingScaleThetaEnum);
     880        this->FindParam(&alpha, CalvingAmpAlphaEnum);
     881        this->FindParam(&midp, CalvingMidpointEnum);
     882        this->FindParam(&nonlinear_law, CalvingNonlinearLawEnum);
     883
     884        /* Start looping on the number of vertices: */
     885        GaussTria* gauss=new GaussTria();
     886        for(int iv=0;iv<NUMVERTICES;iv++){
     887                gauss->GaussVertex(iv);
     888
     889                /*Get velocity components and thickness*/
     890                B_input->GetInputValue(&B,gauss);
     891                n_input->GetInputValue(&n,gauss);
     892                vx_input->GetInputValue(&vx,gauss);
     893                vy_input->GetInputValue(&vy,gauss);
     894                gr_input->GetInputValue(&groundedice,gauss);
     895                bs_input->GetInputValue(&bed,gauss);
     896                H_input->GetInputValue(&thickness,gauss);
     897                smax_fl_input->GetInputValue(&sigma_max_floating,gauss);
     898                smax_gr_input->GetInputValue(&sigma_max_grounded,gauss);
     899                vel=sqrt(vx*vx+vy*vy)+1.e-14;
     900                sl_input->GetInputValue(&sealevel,gauss);
     901                mrate_input->GetInputValue(&mrate,gauss);
     902
     903                /*Compute strain rate and viscosity: */
     904                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     905
     906                /*Get Eigen values*/
     907                Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
     908                _assert_(!xIsNan<IssmDouble>(lambda1));
     909                _assert_(!xIsNan<IssmDouble>(lambda2));
     910
     911                /*Process Eigen values (only account for extension)*/
     912                lambda1 = max(lambda1,0.);
     913                lambda2 = max(lambda2,0.);
     914
     915                /*Calculate sigma_vm*/
     916                epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
     917                sigma_vm[iv]  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
     918
     919                /*Tensile stress threshold*/
     920                if(groundedice<0)
     921                 sigma_max = sigma_max_floating;
     922                else
     923                 sigma_max = sigma_max_grounded;
     924
     925                switch (use_parameter) {
     926                        case 0:
     927                                /* bed elevation */
     928                                paramX = bed;
     929                                break;
     930                        case 1:
     931                                /* Height above floatation */
     932                                if (bed>sealevel)       paramX = 0.0;
     933                                else paramX = thickness - (rho_water/rho_ice) * (sealevel-bed);
     934                                break;
     935                        case 2:
     936                                /* Thicknese */
     937                                paramX = thickness;
     938                                break;
     939                        case 4:
     940                                /* bed elevation+linear curve fitting */
     941                                paramX = bed;
     942                                break;
     943                        case -1:
     944                                /* use nothing, just the mrate*/
     945                                break;
     946                        default:
     947                                _error_("The parameter is not supported yet!");
     948                }
     949
     950                /* Compute the hyperbolic tangent function */
     951                if ((use_parameter>-0.5) & (use_parameter<3)) {
     952                        gamma = 0.5*theta*(1.0-tanh((paramX+midp)/alpha))+(1.0-theta);
     953                        if (gamma<0.0) gamma =0.0;
     954                }
     955                else if (use_parameter>=3) {
     956                        gamma = alpha*paramX + theta;
     957                        if (gamma > 1.0) gamma = 1.0;
     958                        if (gamma < 0.0) gamma = 0.0;
     959                }
     960                else gamma = 1;
     961
     962                /*-------------------------------------------*/
     963                if (nonlinear_law) {
     964                        /*This von Mises type has too strong positive feedback with vel included
     965                         * calvingrate[iv] = (mrate+sigma_vm[iv]*vel/sigma_max)*gamma;
     966                         */
     967                        Hab = thickness - (rho_water/rho_ice) * (sealevel-bed);
     968                        if (Hab < 0.) Hab = 0.;
     969                        if (bed > sealevel) Hab = 0.;
     970
     971                        calvingrate[iv] = (mrate+Hab/sigma_max)*gamma;
     972                }
     973                else {
     974                        calvingrate[iv] = mrate*gamma;
     975                }
     976                calvingratex[iv] = calvingrate[iv]*vx/vel;
     977                calvingratey[iv] = calvingrate[iv]*vy/vel;
     978        }
     979
     980        /*Add input*/
     981        this->AddInput(CalvingratexEnum,&calvingratex[0],P1DGEnum);
     982        this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
     983        this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
     984        this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
     985
     986        /*Clean up and return*/
     987        delete gauss;
    839988}
    840989/*}}}*/
    … …  
    42564405                        calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    42574406                        break;
     4407                case CalvingParameterizationEnum:
     4408                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     4409                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     4410                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
     4411                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     4412                        break;
    42584413                default:
    42594414                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
    … …  
    44184573                                for(i=0;i<dim;i++) m[i]=0.;
    44194574                                break;
     4575                        case CalvingParameterizationEnum:
     4576            lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
     4577            if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss);
     4578            calvingrate_input->GetInputValue(&calvingrate,&gauss);
     4579            meltingrate_input->GetInputValue(&meltingrate,&gauss);
     4580                                meltingrate = 0.;
     4581
     4582            norm_dlsf=0.;
     4583            for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     4584            norm_dlsf=sqrt(norm_dlsf);
     4585
     4586            if(norm_dlsf>1.e-10)
     4587             for(i=0;i<dim;i++){
     4588                c[i]=calvingrate*dlsf[i]/norm_dlsf;
     4589                m[i]=meltingrate*dlsf[i]/norm_dlsf;
     4590             }
     4591            else
     4592             for(i=0;i<dim;i++){
     4593                c[i]=0.; m[i]=0.;
     4594             }
     4595            break;
    44204596
    44214597                        default:

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

    r26830 r26850  
    6060                void                    CalvingFluxLevelset();
    6161                void                    CalvingMeltingFluxLevelset();
     62                void                    CalvingRateParameterization();
    6263                IssmDouble  CharacteristicLength(void);
    6364                void        ComputeBasalStress(void);

  • issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp

    r26784 r26850  
    3838                        femmodel->ElementOperationx(&Element::CalvingRateTest);
    3939                        break;
     40                case CalvingParameterizationEnum:
     41                        femmodel->ElementOperationx(&Element::CalvingRateParameterization);
     42                        break;
    4043                default:
    4144                        _error_("Caving law "<<EnumToStringx(calvinglaw)<<" not supported yet");

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

    r26837 r26850  
    107107        CalvingMinthicknessEnum,
    108108        CalvingTestSpeedfactorEnum,
     109        CalvingUseParamEnum,
     110        CalvingScaleThetaEnum,
     111        CalvingAmpAlphaEnum,
     112        CalvingMidpointEnum,
     113        CalvingNonlinearLawEnum,
    109114        ConfigurationTypeEnum,
    110115        ConstantsGEnum,
    … …  
    12471252        CalvingLevermannEnum,
    12481253        CalvingTestEnum,
     1254        CalvingParameterizationEnum,
    12491255        CalvingVonmisesEnum,
    12501256        CfdragcoeffabsgradEnum,

  • issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

    r26836 r26850  
    268268                case 6: return CalvingCrevasseDepthEnum;
    269269                case 7: return CalvingDev2Enum;
     270                case 9: return CalvingParameterizationEnum;
    270271                case 8: return CalvingTestEnum;
    271272                default: _error_("Marshalled Calving law code \""<<enum_in<<"\" not supported yet");
Note: See TracChangeset for help on using the changeset viewer.