![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
source:
issm/oecreview/Archive/25834-26739/ISSM-26462-26463.diff
| Last change on this file was 26740, checked in by , 3 years ago | |
|---|---|
| File size: 15.2 KB | |
-
../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
408 408 Input* enthalpy_input = element->GetInput(enthalpy_enum); _assert_(enthalpy_input); 409 409 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input); 410 410 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); 411 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);412 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);413 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);414 411 415 412 /*Build friction element, needed later: */ 416 Friction* friction=new Friction(element, dim);413 Friction* friction=new Friction(element,3); 417 414 418 415 /******** MELTING RATES ************************************//*{{{*/ 419 416 element->NormalBase(&normal_base[0],xyz_list_base); … … 469 466 470 467 /*basal friction*/ 471 468 friction->GetAlpha2(&alpha2,gauss); 472 vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); vz_input->GetInputValue(&vz,gauss);469 friction->GetBasalSlidingSpeeds(&vx, &vy, &vz, gauss); 473 470 basalfriction=alpha2*(vx*vx + vy*vy + vz*vz); 474 471 geothermalflux_input->GetInputValue(&geothermalflux,gauss); 475 472 /* -Mb= Fb-(q-q_geo)/((1-w)*L*rho), and (1-w)*rho=rho_ice, cf Aschwanden 2012, eqs.1, 2, 66*/ … … 946 943 element->GetInputValue(&converged,ConvergedEnum); 947 944 if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum; // use enthalpy from last iteration 948 945 else enthalpy_enum=EnthalpyEnum; // use enthalpy from last time step 949 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);950 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);951 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);952 946 Input* enthalpy_input = element->GetInput(enthalpy_enum); _assert_(enthalpy_input); 953 947 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input); 954 948 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input); … … 987 981 // although one node might have become temperate. So keep heat flux switched on. 988 982 geothermalflux_input->GetInputValue(&geothermalflux,gauss); 989 983 friction->GetAlpha2(&alpha2,gauss); 990 vx_input->GetInputValue(&vx,gauss); 991 vy_input->GetInputValue(&vy,gauss); 992 vz_input->GetInputValue(&vz,gauss); 984 friction->GetBasalSlidingSpeeds(&vx, &vy, &vz, gauss); 993 985 basalfriction=alpha2*(vx*vx+vy*vy+vz*vz); 994 986 heatflux=(basalfriction+geothermalflux)/(rho_ice); 995 987 scalar=gauss->weight*Jdet*heatflux; -
../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
359 359 IssmDouble g = element->FindParam(ConstantsGEnum); 360 360 IssmDouble e_v = element->FindParam(HydrologyEnglacialVoidRatioEnum); 361 361 Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); 362 Input* vx_input = element->GetInput(VxBaseEnum);_assert_(vx_input);363 Input* vy_input = element->GetInput(VyBaseEnum);_assert_(vy_input);364 362 Input* h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input); 365 363 Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); 366 364 Input* b_input = element->GetInput(BedEnum); _assert_(b_input); … … 382 380 element->NodalFunctions(basis,gauss); 383 381 384 382 /*Get input values at gauss points*/ 385 vx_input->GetInputValue(&vx,gauss);386 vy_input->GetInputValue(&vy,gauss);387 383 h_input->GetInputValue(&h,gauss); 388 384 G_input->GetInputValue(&G,gauss); 389 385 B_input->GetInputValue(&B,gauss); … … 394 390 H_input->GetInputValue(&H,gauss); 395 391 396 392 /*Get basal velocity*/ 393 friction->GetBasalSlidingSpeeds(&vx, &vy ,gauss); 397 394 ub = sqrt(vx*vx + vy*vy); 398 395 399 396 /*Compute cavity opening w*/ -
../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp
118 118 iomodel->FetchDataToInput(inputs,elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum); 119 119 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum); 120 120 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum); 121 if(iomodel->domaintype==Domain2DhorizontalEnum){ 122 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum); 123 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum); 124 } 121 125 iomodel->FindConstant(&frictionlaw,"md.friction.law"); 122 126 123 127 /*Friction law variables*/ … … 253 257 Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); 254 258 Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); 255 259 Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); 256 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);257 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);258 260 Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); 259 261 Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); 260 262 Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input); … … 286 288 englacial_input->GetInputValue(&ieb,gauss); 287 289 lr_input->GetInputValue(&lr,gauss); 288 290 br_input->GetInputValue(&br,gauss); 289 vx_input->GetInputValue(&vx,gauss);290 vy_input->GetInputValue(&vy,gauss);291 291 headold_input->GetInputValue(&head_old,gauss); 292 292 293 293 /*Get ice A parameter*/ … … 303 303 304 304 /*Compute frictional heat flux*/ 305 305 friction->GetAlpha2(&alpha2,gauss); 306 vx_input->GetInputValue(&vx,gauss); 307 vy_input->GetInputValue(&vy,gauss); 306 friction->GetBasalSlidingSpeeds(&vx, &vy, gauss); 308 307 frictionheat=alpha2*(vx*vx+vy*vy); 309 308 310 309 /*Get water and ice pressures*/ … … 493 492 Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); 494 493 Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); 495 494 Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); 496 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);497 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);498 495 Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); 499 496 Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); 500 497 … … 522 519 head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss); 523 520 lr_input->GetInputValue(&lr,gauss); 524 521 br_input->GetInputValue(&br,gauss); 525 vx_input->GetInputValue(&vx,gauss);526 vy_input->GetInputValue(&vy,gauss);527 522 528 523 /*Get ice A parameter*/ 529 524 B_input->GetInputValue(&B,gauss); … … 538 533 539 534 /*Compute frictional heat flux*/ 540 535 friction->GetAlpha2(&alpha2,gauss); 541 vx_input->GetInputValue(&vx,gauss); 542 vy_input->GetInputValue(&vy,gauss); 536 friction->GetBasalSlidingSpeeds(&vx, &vy, gauss); 543 537 frictionheat=alpha2*(vx*vx+vy*vy); 544 538 545 539 /*Get water and ice pressures*/ -
../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
3488 3488 3489 3489 /*build friction object, used later on: */ 3490 3490 /*dim=4 is special for HO, which is actually 2.5D*/ 3491 Friction* friction=new Friction(element,dim==3? 4:1);3491 Friction* friction=new Friction(element,dim==3?2.5:1); 3492 3492 3493 3493 /*Recover portion of element that is grounded*/ 3494 3494 if(!(friction_style==SubelementFriction2Enum)) phi=element->GetGroundedPortion(xyz_list_base); -
../trunk-jpl/src/c/classes/Loads/Friction.cpp
17 17 /*Constructors/destructors*/ 18 18 Friction::Friction(){/*{{{*/ 19 19 this->element=NULL; 20 this->dim=0;21 20 this->law=0; 21 this->apply_dim = 1; 22 this->domaintype=-1; 23 this->vx_input=NULL; 24 this->vy_input=NULL; 25 this->vz_input=NULL; 22 26 } 23 27 /*}}}*/ 24 Friction::Friction(Element* element_in,int dim_in){/*{{{*/ 28 Friction::Friction(Element* element_in){/*{{{*/ 29 /* Determine the dimension according to the domain type automatically. 30 * There are exceptions, e.g. HO, which needs the user to specify the dimension used in Friciton.*/ 31 this->element=element_in; 25 32 33 /* Load necessary parameters */ 34 element_in->FindParam(&this->law,FrictionLawEnum); 35 element_in->FindParam(&this->domaintype,DomainTypeEnum); 36 37 /* Load VxBase and VyBase for this special case */ 38 switch(this->domaintype){ 39 case Domain2DhorizontalEnum: 40 this->apply_dim = 2; 41 this->vx_input = element_in->GetInput(VxBaseEnum); _assert_(this->vx_input); 42 this->vy_input = element_in->GetInput(VyBaseEnum); _assert_(this->vy_input); 43 this->vz_input = NULL; 44 break; 45 case Domain2DverticalEnum: 46 this->apply_dim = 2; 47 this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input); 48 this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input); 49 this->vz_input = NULL; 50 break; 51 case Domain3DEnum: 52 this->apply_dim = 3; 53 this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input); 54 this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input); 55 this->vz_input = element_in->GetInput(VzEnum); _assert_(this->vz_input); 56 break; 57 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 58 } 59 } 60 /*}}}*/ 61 Friction::Friction(Element* element_in, double dim){/*{{{*/ 62 26 63 this->element=element_in; 27 this->dim=dim_in; 64 this->apply_dim = dim; 65 66 /* Load necessary parameters */ 28 67 element_in->FindParam(&this->law,FrictionLawEnum); 68 element_in->FindParam(&this->domaintype,DomainTypeEnum); 69 70 /* Load VxBase and VyBase for this special case */ 71 switch(this->domaintype){ 72 case Domain2DhorizontalEnum: 73 this->vx_input = element_in->GetInput(VxBaseEnum); _assert_(this->vx_input); 74 this->vy_input = element_in->GetInput(VyBaseEnum); _assert_(this->vy_input); 75 this->vz_input = NULL; 76 break; 77 case Domain2DverticalEnum: 78 this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input); 79 this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input); 80 this->vz_input = NULL; 81 break; 82 case Domain3DEnum: 83 this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input); 84 this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input); 85 this->vz_input = element_in->GetInput(VzEnum); _assert_(this->vz_input); 86 break; 87 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 88 } 29 89 } 30 90 /*}}}*/ 31 91 Friction::~Friction(){/*{{{*/ … … 35 95 /*methods: */ 36 96 void Friction::Echo(void){/*{{{*/ 37 97 _printf_("Friction:\n"); 38 _printf_(" dim: " << this->dim<< "\n");98 _printf_(" Domain type: " << this->domaintype<< "\n"); 39 99 } 40 100 /*}}}*/ 41 101 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/ … … 854 914 855 915 /*diverse*/ 856 916 IssmDouble vx,vy,vz,vmag; 917 918 this->vx_input->GetInputValue(&vx, gauss); 919 this->vy_input->GetInputValue(&vy, gauss); 920 921 if ((this->vz_input == NULL) || (this->apply_dim<3)) vz = 0.0; 922 else this->vz_input->GetInputValue(&vz, gauss); 857 923 858 switch(dim){ 859 case 1: 860 element->GetInputValue(&vx,gauss,VxEnum); 861 vmag=sqrt(vx*vx); 862 break; 863 case 2: 864 element->GetInputValue(&vx,gauss,VxBaseEnum); 865 element->GetInputValue(&vy,gauss,VyBaseEnum); 866 vmag=sqrt(vx*vx+vy*vy); 867 break; 868 case 3: 869 element->GetInputValue(&vx,gauss,VxEnum); 870 element->GetInputValue(&vy,gauss,VyEnum); 871 element->GetInputValue(&vz,gauss,VzEnum); 872 vmag=sqrt(vx*vx+vy*vy+vz*vz); 873 break; 874 case 4: 875 /* This is for HO 3D case, since it requires the horizontal velocities */ 876 element->GetInputValue(&vx,gauss,VxEnum); 877 element->GetInputValue(&vy,gauss,VyEnum); 878 vmag=sqrt(vx*vx+vy*vy); 879 break; 880 default: 881 _error_("not supported"); 882 } 924 if (this->apply_dim<2) vy = 0.0; 883 925 926 vmag = sqrt(vx*vx+vy*vy+vz*vz); 884 927 return vmag; 928 }/*}}}*/ 929 void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss){/*{{{*/ 885 930 931 this->vx_input->GetInputValue(pvx, gauss); 932 /*Checks*/ 933 _assert_(!xIsNan<IssmDouble>(*pvx)); 934 _assert_(!xIsInf<IssmDouble>(*pvx)); 886 935 }/*}}}*/ 936 void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss){/*{{{*/ 937 938 this->vx_input->GetInputValue(pvx, gauss); 939 this->vy_input->GetInputValue(pvy, gauss); 940 /*Checks*/ 941 _assert_(!xIsNan<IssmDouble>(*pvx)); 942 _assert_(!xIsInf<IssmDouble>(*pvx)); 943 _assert_(!xIsNan<IssmDouble>(*pvy)); 944 _assert_(!xIsInf<IssmDouble>(*pvy)); 945 }/*}}}*/ 946 void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss){/*{{{*/ 947 948 this->vx_input->GetInputValue(pvx, gauss); 949 this->vy_input->GetInputValue(pvy, gauss); 950 this->vz_input->GetInputValue(pvz, gauss); 951 /*Checks*/ 952 _assert_(!xIsNan<IssmDouble>(*pvx)); 953 _assert_(!xIsInf<IssmDouble>(*pvx)); 954 _assert_(!xIsNan<IssmDouble>(*pvy)); 955 _assert_(!xIsInf<IssmDouble>(*pvy)); 956 _assert_(!xIsNan<IssmDouble>(*pvz)); 957 _assert_(!xIsInf<IssmDouble>(*pvz)); 958 }/*}}}*/ -
../trunk-jpl/src/c/classes/Loads/Friction.h
14 14 15 15 public: 16 16 Element* element; 17 int dim;18 17 int law; 18 int domaintype; 19 double apply_dim; 20 Input* vx_input; 21 Input* vy_input; 22 Input* vz_input; 19 23 20 24 /*methods: */ 21 25 Friction(); 22 Friction(Element* element_in,int dim_in); 26 Friction(Element* element_in); 27 Friction(Element* element_in, double dim); 23 28 ~Friction(); 24 29 25 30 void Echo(void); … … 46 51 47 52 IssmDouble EffectivePressure(Gauss* gauss); 48 53 IssmDouble VelMag(Gauss* gauss); 54 void GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss); 55 void GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss); 56 void GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss); 49 57 }; 50 58 51 59 #endif /* _FRICTION_H_ */
Note:
See TracBrowser
for help on using the repository browser.
