Changeset 19033


Ignore:
Timestamp:
01/22/15 14:51:25 (10 years ago)
Author:
Mathieu Morlighem
Message:

NEW: testing new calving law

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

Legend:

Unmodified
Added
Removed

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

    r19025 r19033  
    6868                                iomodel->FetchDataToInput(elements,CalvingpiMeltingrateEnum);
    6969                                break;
     70                        case CalvingDevEnum:
     71                                iomodel->FetchDataToInput(elements,CalvingpiCoeffEnum);
     72                                iomodel->FetchDataToInput(elements,CalvingMeltingrateEnum);
     73                                break;
    7074                        default:
    7175                                _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
    … …  
    174178
    175179        /*Load velocities*/
    176         if(domaintype==Domain2DhorizontalEnum){
    177                 vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
    178                 vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
    179         }
    180         else{
    181                 if(dim==1){
     180        switch(domaintype){
     181                case Domain2DverticalEnum:
    182182                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
    183                 }
    184                 if(dim==2){
     183                        break;
     184                case Domain2DhorizontalEnum:
     185                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     186                        vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
     187                        break;
     188                case Domain3DEnum:
    185189                        vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
    186190                        vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
    187                 }
     191                        break;
     192                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    188193        }
    189194
    … …  
    198203                                break;
    199204                        case CalvingLevermannEnum:
    200                                 if(domaintype==Domain2DhorizontalEnum){
    201                                         calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    202                                         calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    203                                 }
    204                                 else{
    205                                         if(dim==1){
     205                                switch(domaintype){
     206                                        case Domain2DverticalEnum:
    206207                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    207                                         }
    208                                         if(dim==2){
     208                                                break;
     209                                        case Domain2DhorizontalEnum:
     210                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     211                                                calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     212                                                break;
     213                                        case Domain3DEnum:
    209214                                                calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
    210215                                                calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
    211                                         }
     216                                                break;
     217                                        default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    212218                                }
    213219                                meltingrate_input = basalelement->GetInput(CalvinglevermannMeltingrateEnum);     _assert_(meltingrate_input);
    214220                                break;
    215221                        case CalvingPiEnum:
    216                                 if(domaintype==Domain2DhorizontalEnum){
    217                                         calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    218                                         calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    219                                 }
    220                                 else{
    221                                         if(dim==1){
     222                                switch(domaintype){
     223                                        case Domain2DverticalEnum:
    222224                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    223                                         }
    224                                         if(dim==2){
     225                                                break;
     226                                        case Domain2DhorizontalEnum:
     227                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     228                                                calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     229                                                break;
     230                                        case Domain3DEnum:
    225231                                                calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
    226232                                                calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
    227                                         }
     233                                                break;
     234                                        default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    228235                                }
    229236                                meltingrate_input = basalelement->GetInput(CalvingpiMeltingrateEnum);     _assert_(meltingrate_input);
     237                                break;
     238                        case CalvingDevEnum:
     239                                switch(domaintype){
     240                                        case Domain2DverticalEnum:
     241                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     242                                                break;
     243                                        case Domain2DhorizontalEnum:
     244                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     245                                                calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     246                                                break;
     247                                        case Domain3DEnum:
     248                                                calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     249                                                calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     250                                                break;
     251                                        default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     252                                }
     253                                meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    230254                                break;
    231255                        default:
    … …  
    255279                GetB(B,basalelement,xyz_list,gauss);
    256280                GetBprime(Bprime,basalelement,xyz_list,gauss);
    257                 vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity
     281                vx_input->GetInputValue(&v[0],gauss);
    258282                vy_input->GetInputValue(&v[1],gauss);
    259283
    … …  
    273297                                        if(norm_dlsf>1.e-10)
    274298                                         for(i=0;i<dim;i++){
    275                                           c[i]=calvingrate*dlsf[i]/norm_dlsf;
    276                                           m[i]=meltingrate*dlsf[i]/norm_dlsf;
     299                                          c[i]=calvingrate*dlsf[i]/norm_dlsf; m[i]=meltingrate*dlsf[i]/norm_dlsf;
    277300                                         }
    278301                                        else
    279302                                         for(i=0;i<dim;i++){
    280                                                  c[i]=0.;
    281                                                  m[i]=0.;
     303                                                 c[i]=0.; m[i]=0.;
    282304                                         }
    283 
    284305                                        break;
    285306
    286307                                case CalvingLevermannEnum:
    287                                         calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity
     308                                case CalvingPiEnum:
     309                                case CalvingDevEnum:
     310                                        calvingratex_input->GetInputValue(&c[0],gauss);
    288311                                        if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
    289312                                        meltingrate_input->GetInputValue(&meltingrate,gauss);
    290 
    291313                                        norm_calving=0.;
    292314                                        for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
    293315                                        norm_calving=sqrt(norm_calving)+1.e-14;
    294 
    295316                                        for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
    296                                        
    297                                         break;
    298 
    299                                 case CalvingPiEnum:
    300                                         calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity
    301                                         if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
    302                                         meltingrate_input->GetInputValue(&meltingrate,gauss);
    303 
    304                                         norm_calving=0.;
    305                                         for(i=0;i<dim;i++) norm_calving+=c[i]*c[i];
    306                                         norm_calving=sqrt(norm_calving)+1.e-14;
    307 
    308                                         for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
    309 
    310317                                        break;
    311318

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

    r18928 r19033  
    15681568                element->NodalFunctions(basis, gauss);
    15691569
    1570                 thickness_input->GetInputValue(&thickness,gauss);
     1570                thickness_input->GetInputValue(&thickness,gauss); _assert_(thickness>0);
    15711571                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
    15721572

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

    r18968 r19033  
    171171                virtual void        CalvingRateLevermann(void)=0;
    172172                virtual void       CalvingRatePi(void)=0;
     173                virtual void       CalvingRateDev(void){_error_("not implemented yet");};
    173174                virtual IssmDouble CharacteristicLength(void)=0;
    174175                virtual void       ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;

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

    r19029 r19033  
    326326        delete gauss;
    327327
     328}
     329/*}}}*/
     330void       Tria::CalvingRateDev(){/*{{{*/
     331
     332        IssmDouble  xyz_list[NUMVERTICES][3];
     333        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     334        IssmDouble  calvingratex[NUMVERTICES];
     335        IssmDouble  calvingratey[NUMVERTICES];
     336        IssmDouble  calvingrate[NUMVERTICES];
     337        IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
     338
     339        /* Get node coordinates and dof list: */
     340        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     341
     342        /*Retrieve all inputs and parameters we will need*/
     343        Input* vx_input=inputs->GetInput(VxEnum);        _assert_(vx_input);
     344        Input* vy_input=inputs->GetInput(VyEnum);        _assert_(vy_input);
     345
     346        /* Start looping on the number of vertices: */
     347        GaussTria* gauss=new GaussTria();
     348        for(int iv=0;iv<NUMVERTICES;iv++){
     349                gauss->GaussVertex(iv);
     350
     351                /*Get velocity components and thickness*/
     352                vx_input->GetInputValue(&vx,gauss);
     353                vy_input->GetInputValue(&vy,gauss);
     354                vel=sqrt(vx*vx+vy*vy)+1.e-14;
     355
     356                /*Compute strain rate and viscosity: */
     357                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     358
     359                /*Get Eigen values*/
     360                Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
     361
     362                /*Process Eigen values*/
     363                lambda1>0? lambda1 = pow(lambda1,.3) : lambda1=0.;
     364                lambda2>0? lambda2 = pow(lambda2,.3) : lambda2=0.;
     365                lambda1 = lambda1*5.e-2;
     366                lambda2 = lambda2*5.e-2;
     367                _assert_(!xIsNan<IssmDouble>(lambda1));
     368                _assert_(!xIsNan<IssmDouble>(lambda2));
     369
     370                /*Assign values*/
     371                //calvingratex[iv]=ex*lambda1 - ey*lambda2;
     372                //calvingratey[iv]=ey*lambda1 + ex*lambda2;
     373                calvingratex[iv]=vx/vel*(lambda1 + lambda2);
     374                calvingratey[iv]=vy/vel*(lambda1 + lambda2);
     375                calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
     376        }
     377
     378        /*Add input*/
     379        this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
     380        this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
     381        this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
     382
     383        /*Clean up and return*/
     384        delete gauss;
    328385}
    329386/*}}}*/

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

    r18968 r19033  
    5252                void        AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
    5353                void                    CalvingRateLevermann();
     54                void                    CalvingRatePi();
     55                void                    CalvingRateDev();
    5456                IssmDouble  CharacteristicLength(void);
    5557                void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
    … …  
    5860                void        ComputeStressTensor();
    5961                void        ComputeSurfaceNormalVelocity();
    60                 void                    CalvingRatePi();
    6162                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
    6263                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r19019 r19033  
    669669        /*Now, update degrees of freedoms: */
    670670        NodesDofx(nodes,parameters,analysis_type);
     671
     672}
     673/*}}}*/
     674void FemModel::ElementOperationx(void (Element::*function)(void)){ /*{{{*/
     675
     676        for(int i=0;i<elements->Size();i++){
     677                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
     678                (element->*function)();
     679        }
    671680
    672681}
    … …  
    17431752}
    17441753/*}}}*/
     1754void FemModel::CalvingRateDevx(){/*{{{*/
     1755
     1756        for(int i=0;i<elements->Size();i++){
     1757                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
     1758                element->CalvingRateDev();
     1759        }
     1760}
     1761/*}}}*/
    17451762void FemModel::StrainRateparallelx(){/*{{{*/
    17461763

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r19019 r19033  
    6262
    6363                /*Modules*/
     64                void ElementOperationx(void (Element::*function)(void));
    6465                void GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,IssmDouble* ug);
    6566                void MassFluxx(IssmDouble* presponse);
    … …  
    8889                void CalvingRateLevermannx();
    8990                void CalvingRatePix();
     91                void CalvingRateDevx();
    9092                #ifdef  _HAVE_DAKOTA_
    9193                void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);

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

    r19025 r19033  
    103103
    104104                if(islevelset){
    105                         if(iscalving && calvinglaw==CalvingLevermannEnum){
    106                                 if(VerboseSolution()) _printf0_("   computing calving rate\n");
    107                                 femmodel->StrainRateparallelx();
    108                                 femmodel->StrainRateperpendicularx();
    109                                 femmodel->CalvingRateLevermannx();
     105                        if(iscalving){
     106                                switch(calvinglaw){
     107                                        case DefaultCalvingEnum:
     108                                                break;
     109                                        case CalvingLevermannEnum:
     110                                                if(VerboseSolution()) _printf0_("   computing Levermann's calving rate\n");
     111                                                femmodel->StrainRateparallelx();
     112                                                femmodel->StrainRateperpendicularx();
     113                                                femmodel->CalvingRateLevermannx();
     114                                                break;
     115                                        case CalvingPiEnum:
     116                                                if(VerboseSolution()) _printf0_("   computing Pi's calving rate\n");
     117                                                femmodel->StrainRateparallelx();
     118                                                femmodel->DeviatoricStressx();
     119                                                femmodel->CalvingRatePix();
     120                                                break;
     121                                        case CalvingDevEnum:
     122                                                femmodel->CalvingRateDevx();
     123                                                femmodel->ElementOperationx(&Element::CalvingRateDev);
     124                                                break;
     125                                        default:
     126                                                _error_("Caving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     127                                }
    110128                        }
    111                         if(iscalving && calvinglaw==CalvingPiEnum){
    112                                 if(VerboseSolution()) _printf0_("   computing calving rate\n");
    113                                 femmodel->StrainRateparallelx();
    114                                 femmodel->DeviatoricStressx();
    115                                 femmodel->CalvingRatePix();
    116                         }
    117                         if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
     129                        if(VerboseSolution()) _printf0_("   computing levelset transport\n");
    118130                        /* smoothen slope of lsf for computation of normal on ice domain*/
    119131                        levelsetfunctionslope_core(femmodel);
Note: See TracChangeset for help on using the changeset viewer.