Changeset 22199


Ignore:
Timestamp:
10/27/17 17:06:52 (7 years ago)
Author:
jbondzio
Message:

NEW: New CalvingDev for calving of tidewaterglaciers

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

Legend:

Unmodified
Added
Removed

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

    r22192 r22199  
    7575                        iomodel->FetchDataToInput(elements,"md.calving.water_height",WaterheightEnum);
    7676                        break;
     77                case CalvingDev2Enum:
     78                        iomodel->FetchDataToInput(elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum);
     79                        iomodel->FetchDataToInput(elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum);
     80                        break;
    7781                default:
    7882                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
    … …  
    100104                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.critical_fraction",CalvingCrevasseDepthEnum));
    101105                        break;
     106                case CalvingDev2Enum:
     107                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.height_above_floatation",CalvingHeightAboveFloatationEnum));
     108                        break;
    102109                default:
    103110                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
    … …  
    154161        IssmDouble vel;
    155162        IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice;
    156         IssmDouble calvingmax;
     163        IssmDouble calvingmax, calvinghaf, heaviside, haf_eps;
    157164        IssmDouble* xyz_list = NULL;
    158165
    … …  
    259266                        meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    260267                        break;
     268                case CalvingDev2Enum:
     269                        basalelement->FindParam(&calvinghaf,CalvingHeightAboveFloatationEnum);
     270                        lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     271                        if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     272                        calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
     273                        gr_input=basalelement->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input);
     274                        break;
    261275                default:
    262276                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
    … …  
    315329                                 }
    316330                                break;
    317 
     331                               
    318332                        case CalvingLevermannEnum:
    319333                                calvingratex_input->GetInputValue(&c[0],gauss);
    … …  
    394408                                 }
    395409                                break;
     410
     411                        case CalvingDev2Enum:
     412                                {
     413                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     414                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     415                                calvingrate_input->GetInputValue(&calvingrate,gauss);
     416                                gr_input->GetInputValue(&groundedice,gauss);
     417
     418                                //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function.
     419                                vel=sqrt(v[0]*v[0] + v[1]*v[1]);
     420                                haf_eps=10.;
     421                                if(groundedice-calvinghaf<=-haf_eps){
     422                                        // ice floats freely below calvinghaf: calve freely
     423                                }
     424                                else if(groundedice-calvinghaf>=haf_eps){
     425                                        // ice is well above calvinghaf -> no calving back, i.e. limit calving rate to ice velocity
     426                                        calvingrate=min(calvingrate,vel);
     427                                }
     428                                else{ // ice is close to calvinghaf: smooth transition between limitation and free calving.
     429                                        heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI);
     430                                        calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate;
     431                                }
     432
     433                                norm_dlsf=0.;
     434                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     435                                norm_dlsf=sqrt(norm_dlsf);
     436
     437                                if(norm_dlsf>1.e-10)
     438                                 for(i=0;i<dim;i++){
     439                                         c[i]=calvingrate*dlsf[i]/norm_dlsf;
     440                                         m[i]=0.;
     441                                 }
     442                                else
     443                                 for(i=0;i<dim;i++){
     444                                         c[i]=0.;
     445                                         m[i]=0.;
     446                                 }
     447                                break;
     448                                }
    396449
    397450                        default:

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

    r22197 r22199  
    23532353                                        break;
    23542354                                case CalvingVonmisesEnum:
     2355                                case CalvingDev2Enum:
    23552356                                        this->CalvingRateVonmises();
    23562357                                        break;

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

    r22192 r22199  
    3232                        break;
    3333                case CalvingVonmisesEnum:
     34                case CalvingDev2Enum:
    3435                        femmodel->CalvingRateVonmisesx();
    3536                        femmodel->ElementOperationx(&Element::CalvingRateVonmises);

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

    r22193 r22199  
    256256        CalvingCalvingrateEnum,
    257257        CalvingMeltingrateEnum,
     258        CalvingHeightAboveFloatationEnum,
    258259        CalvingLevermannEnum,
    259260        CalvingVonmisesEnum,
    … …  
    261262        CalvingHabEnum,
    262263        CalvingCrevasseDepthEnum,
     264        CalvingDev2Enum,
    263265        SurfaceCrevasseEnum,
    264266        BasalCrevasseEnum,

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

    r22193 r22199  
    262262                case CalvingCalvingrateEnum : return "CalvingCalvingrate";
    263263                case CalvingMeltingrateEnum : return "CalvingMeltingrate";
     264                case CalvingHeightAboveFloatationEnum : return "CalvingHeightAboveFloatation";
    264265                case CalvingLevermannEnum : return "CalvingLevermann";
    265266                case CalvingVonmisesEnum : return "CalvingVonmises";
    … …  
    267268                case CalvingHabEnum : return "CalvingHab";
    268269                case CalvingCrevasseDepthEnum : return "CalvingCrevasseDepth";
     270                case CalvingDev2Enum : return "CalvingDev2";
    269271                case SurfaceCrevasseEnum : return "SurfaceCrevasse";
    270272                case BasalCrevasseEnum : return "BasalCrevasse";

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

    r22193 r22199  
    268268              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
    269269              else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
     270              else if (strcmp(name,"CalvingHeightAboveFloatation")==0) return CalvingHeightAboveFloatationEnum;
    270271              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    271272              else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
    … …  
    273274              else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
    274275              else if (strcmp(name,"CalvingCrevasseDepth")==0) return CalvingCrevasseDepthEnum;
     276              else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
    275277              else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum;
    276278              else if (strcmp(name,"BasalCrevasse")==0) return BasalCrevasseEnum;
    … …  
    381383              else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum;
    382384              else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
    383               else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum;
    384               else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
     388              if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum;
     389              else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum;
     390              else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
    389391              else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum;
    390392              else if (strcmp(name,"TransientIsslr")==0) return TransientIsslrEnum;
    … …  
    504506              else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
    505507              else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
    506               else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
    507               else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
     511              if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
     512              else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
     513              else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
    512514              else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
    513515              else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
    … …  
    627629              else if (strcmp(name,"Outputdefinition10")==0) return Outputdefinition10Enum;
    628630              else if (strcmp(name,"Outputdefinition11")==0) return Outputdefinition11Enum;
    629               else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum;
    630               else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
     634              if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum;
     635              else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
     636              else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
    635637              else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
    636638              else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
    … …  
    750752              else if (strcmp(name,"Gradient")==0) return GradientEnum;
    751753              else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
    752               else if (strcmp(name,"Gset")==0) return GsetEnum;
    753               else if (strcmp(name,"Index")==0) return IndexEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"Indexed")==0) return IndexedEnum;
     757              if (strcmp(name,"Gset")==0) return GsetEnum;
     758              else if (strcmp(name,"Index")==0) return IndexEnum;
     759              else if (strcmp(name,"Indexed")==0) return IndexedEnum;
    758760              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
    759761              else if (strcmp(name,"Nodal")==0) return NodalEnum;
    … …  
    873875              else if (strcmp(name,"AmrIceFrontResolution")==0) return AmrIceFrontResolutionEnum;
    874876              else if (strcmp(name,"AmrIceFrontDistance")==0) return AmrIceFrontDistanceEnum;
    875               else if (strcmp(name,"AmrThicknessErrorResolution")==0) return AmrThicknessErrorResolutionEnum;
    876               else if (strcmp(name,"AmrThicknessErrorThreshold")==0) return AmrThicknessErrorThresholdEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"AmrDeviatoricErrorResolution")==0) return AmrDeviatoricErrorResolutionEnum;
     880              if (strcmp(name,"AmrThicknessErrorResolution")==0) return AmrThicknessErrorResolutionEnum;
     881              else if (strcmp(name,"AmrThicknessErrorThreshold")==0) return AmrThicknessErrorThresholdEnum;
     882              else if (strcmp(name,"AmrDeviatoricErrorResolution")==0) return AmrDeviatoricErrorResolutionEnum;
    881883              else if (strcmp(name,"AmrDeviatoricErrorThreshold")==0) return AmrDeviatoricErrorThresholdEnum;
    882884              else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
    … …  
    996998              else if (strcmp(name,"GiaIvinsAnalysis")==0) return GiaIvinsAnalysisEnum;
    997999              else if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum;
    998               else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
    999               else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum;
     1003              if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
     1004              else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum;
     1005              else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum;
    10041006              else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
    10051007              else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;

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

    r22196 r22199  
    153153                case 5: return CalvingHabEnum;
    154154                case 6: return CalvingCrevasseDepthEnum;
     155                case 7: return CalvingDev2Enum;
    155156                default: _error_("Marshalled Calving law code \""<<enum_in<<"\" not supported yet");
    156157        }
Note: See TracChangeset for help on using the changeset viewer.