Changeset 20723


Ignore:
Timestamp:
06/10/16 11:15:57 (9 years ago)
Author:
Mathieu Morlighem
Message:

CHG: spc based min thickness calving

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

Legend:

Unmodified
Added
Removed

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

    r20719 r20723  
    113113        if(save_results){
    114114                if(VerboseSolution()) _printf0_("   saving results\n");
    115                 int outputs[2] = {MaskIceLevelsetEnum, CalvingCalvingrateEnum};
    116                 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
     115                int outputs[1] = {MaskIceLevelsetEnum};
     116                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
    117117        }
    118118}/*}}}*/
    … …  
    198198                case DefaultCalvingEnum:
    199199                case CalvingDevEnum:
    200                 case CalvingMinthicknessEnum:
    201200                        lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    202201                        if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    … …  
    221220                        meltingrate_input = basalelement->GetInput(CalvinglevermannMeltingrateEnum);     _assert_(meltingrate_input);
    222221                        break;
     222                case CalvingMinthicknessEnum:
     223                        lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     224                        if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     225                        meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     226                        break;
    223227                default:
    224228                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
    … …  
    253257                        case DefaultCalvingEnum:
    254258                        case CalvingDevEnum:
    255                         case CalvingMinthicknessEnum:
    256259                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    257260                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    … …  
    281284                                norm_calving=sqrt(norm_calving)+1.e-14;
    282285                                for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
     286                                break;
     287
     288                        case CalvingMinthicknessEnum:
     289                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     290                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     291                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     292
     293                                norm_dlsf=0.;
     294                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     295                                norm_dlsf=sqrt(norm_dlsf);
     296
     297                                if(norm_dlsf>1.e-10)
     298                                 for(i=0;i<dim;i++){
     299                                         c[i]=0.;
     300                                         m[i]=meltingrate*dlsf[i]/norm_dlsf;
     301                                 }
     302                                else
     303                                 for(i=0;i<dim;i++){
     304                                         c[i]=0.;
     305                                         m[i]=0.;
     306                                 }
    283307                                break;
    284308
    … …  
    621645}/*}}}*/
    622646void           LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
     647
     648        /*Intermediaries*/
     649        int         calvinglaw;
     650        IssmDouble  min_thickness,thickness;
     651        femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
     652
     653        if(calvinglaw==CalvingMinthicknessEnum){
     654
     655                /*Get minimum thickness threshold*/
     656                femmodel->parameters->FindParam(&min_thickness,CalvingMinthicknessEnum);
     657
     658                /*Loop over all elements of this partition*/
     659                for(int i=0;i<femmodel->elements->Size();i++){
     660                        Element* element  = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     661
     662                        int      numnodes = element->GetNumberOfNodes();       
     663                        Gauss*   gauss    = element->NewGauss();
     664                        Input*   H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
     665
     666                        /*Potentially constrain nodes of this element*/
     667                        for(int in=0;in<numnodes;in++){
     668                                gauss->GaussNode(element->GetElementType(),in);
     669                                Node* node=element->GetNode(in);
     670                                H_input->GetInputValue(&thickness,gauss);
     671                                if(thickness<min_thickness){
     672                                        node->ApplyConstraint(0,+1.);
     673                                }
     674                                else {
     675                                        /* no ice, set no spc */
     676                                        node->DofInFSet(0);
     677                                }
     678                        }
     679                        delete gauss;
     680                }
     681        }
    623682        /*Default, do nothing*/
    624683        return;

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

    r20719 r20723  
    20692069                        switch(calvinglaw){
    20702070                                case DefaultCalvingEnum:
     2071                                case CalvingMinthicknessEnum:
    20712072                                        //do nothing
    20722073                                        break;
    … …  
    20762077                                case CalvingDevEnum:
    20772078                                        this->CalvingRateDev();
    2078                                         break;
    2079                                 case CalvingMinthicknessEnum:
    2080                                         this->CalvingRateMinthickness();
    20812079                                        break;
    20822080                                default:

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

    r20719 r20723  
    176176                virtual void        CalvingRateLevermann(void)=0;
    177177                virtual void       CalvingRateDev(void){_error_("not implemented yet");};
    178                 virtual void       CalvingRateMinthickness(void){_error_("not implemented yet");};
    179178                virtual void       WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");};
    180179                virtual void       ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");};

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

    r20719 r20723  
    327327                calvingratey[iv]=vy*sigma_vm/sigma_max;
    328328                calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
    329         }
    330 
    331         /*Add input*/
    332         this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
    333         this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
    334         this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
    335 
    336         /*Clean up and return*/
    337         delete gauss;
    338 }
    339 /*}}}*/
    340 void       Tria::CalvingRateMinthickness(){/*{{{*/
    341 
    342         IssmDouble  calvingratex[NUMVERTICES];
    343         IssmDouble  calvingratey[NUMVERTICES];
    344         IssmDouble  calvingrate[NUMVERTICES];
    345         IssmDouble  vx,vy,H,minH;
    346 
    347         /*Retrieve all inputs and parameters we will need*/
    348         Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
    349         Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
    350         Input* H_input  = inputs->GetInput(ThicknessEnum); _assert_(H_input);
    351         parameters->FindParam(&minH,CalvingMinthicknessEnum);
    352 
    353         /* Start looping on the number of vertices: */
    354         GaussTria* gauss=new GaussTria();
    355         for(int iv=0;iv<NUMVERTICES;iv++){
    356                 gauss->GaussVertex(iv);
    357 
    358                 /*Get velocity components and thickness*/
    359                 vx_input->GetInputValue(&vx,gauss);
    360                 vy_input->GetInputValue(&vy,gauss);
    361                 H_input->GetInputValue(&H,gauss);
    362 
    363                 /*Assign values*/
    364                 if(H>minH){
    365                         calvingratex[iv]=0.;
    366                         calvingratey[iv]=0.;
    367                         calvingrate[iv]=0.;
    368                 }
    369                 else{
    370                         calvingratex[iv]=vx+1e-2; //counter balance advance + add some retreat
    371                         calvingratey[iv]=vy+1e-2;
    372                         calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
    373                 }
    374329        }
    375330

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

    r20719 r20723  
    5454                void                    CalvingRateLevermann();
    5555                void                    CalvingRateDev();
    56                 void                    CalvingRateMinthickness();
    5756                void                    WriteLevelsetSegment(DataSet* segments);
    5857                void        ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);

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

    r20719 r20723  
    21662166}
    21672167/*}}}*/
    2168 void FemModel::CalvingRateMinthicknessx(){/*{{{*/
    2169 
    2170         for(int i=0;i<elements->Size();i++){
    2171                 Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
    2172                 element->CalvingRateMinthickness();
    2173         }
    2174 }
    2175 /*}}}*/
    21762168void FemModel::StrainRateparallelx(){/*{{{*/
    21772169

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

    r20719 r20723  
    9595                void CalvingRateLevermannx();
    9696                void CalvingRateDevx();
    97                 void CalvingRateMinthicknessx();
    9897                void ResetLevelset();
    9998                #ifdef  _HAVE_DAKOTA_

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

    r20719 r20723  
    1616        switch(calvinglaw){
    1717                case DefaultCalvingEnum:
     18                case CalvingMinthicknessEnum:
    1819                        break;
    1920                case CalvingLevermannEnum:
    … …  
    2728                        femmodel->ElementOperationx(&Element::CalvingRateDev);
    2829                        break;
    29                 case CalvingMinthicknessEnum:
    30                         femmodel->CalvingRateMinthicknessx();
    31                         break;
    3230                default:
    3331                        _error_("Caving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
Note: See TracChangeset for help on using the changeset viewer.