 |
Ice Sheet System Model
4.18
Code documentation
|
#include <EnthalpyAnalysis.h>
|
void | CreateConstraints (Constraints *constraints, IoModel *iomodel) |
|
void | CreateLoads (Loads *loads, IoModel *iomodel) |
|
void | CreateNodes (Nodes *nodes, IoModel *iomodel, bool isamr=false) |
|
int | DofsPerNode (int **doflist, int domaintype, int approximation) |
|
void | UpdateElements (Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type) |
|
void | UpdateParameters (Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum) |
|
void | Core (FemModel *femmodel) |
|
ElementVector * | CreateDVector (Element *element) |
|
ElementMatrix * | CreateJacobianMatrix (Element *element) |
|
ElementMatrix * | CreateKMatrix (Element *element) |
|
ElementMatrix * | CreateKMatrixVolume (Element *element) |
|
ElementMatrix * | CreateKMatrixShelf (Element *element) |
|
ElementVector * | CreatePVector (Element *element) |
|
ElementVector * | CreatePVectorVolume (Element *element) |
|
ElementVector * | CreatePVectorSheet (Element *element) |
|
ElementVector * | CreatePVectorShelf (Element *element) |
|
void | GetBConduct (IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetSolutionFromInputs (Vector< IssmDouble > *solution, Element *element) |
|
void | GradientJ (Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index) |
|
void | InputUpdateFromSolution (IssmDouble *solution, Element *element) |
|
void | UpdateConstraints (FemModel *femmodel) |
|
virtual | ~Analysis () |
|
|
static void | ApplyBasalConstraints (IssmDouble *serial_spc, Element *element) |
|
static void | ComputeBasalMeltingrate (FemModel *femmodel) |
|
static void | ComputeBasalMeltingrate (Element *element) |
|
static void | DrainWaterfraction (FemModel *femmodel) |
|
static void | ComputeWaterfractionDrainage (FemModel *femmodel) |
|
static void | DrainageUpdateWatercolumn (FemModel *femmodel) |
|
static void | DrainageUpdateEnthalpy (FemModel *femmodel) |
|
static IssmDouble | EnthalpyDiffusionParameter (Element *element, IssmDouble enthalpy, IssmDouble pressure) |
|
static IssmDouble | EnthalpyDiffusionParameterVolume (Element *element, int enthalpy_enum) |
|
static void | GetBasalConstraints (Vector< IssmDouble > *vec_spc, Element *element) |
|
static void | GetBasalConstraintsSteadystate (Vector< IssmDouble > *vec_spc, Element *element) |
|
static void | GetBasalConstraintsTransient (Vector< IssmDouble > *vec_spc, Element *element) |
|
static int | GetThermalBasalCondition (Element *element, IssmDouble enthalpy, IssmDouble enthalpy_up, IssmDouble pressure, IssmDouble pressure_up, IssmDouble watercolumn, IssmDouble meltingrate) |
|
static IssmDouble | GetWetIceConductivity (Element *element, IssmDouble enthalpy, IssmDouble pressure) |
|
static void | PostProcessing (FemModel *femmodel) |
|
static IssmDouble | PureIceEnthalpy (Element *element, IssmDouble pressure) |
|
static IssmDouble | TMeltingPoint (Element *element, IssmDouble pressure) |
|
static void | UpdateBasalConstraints (FemModel *femmodel) |
|
Definition at line 12 of file EnthalpyAnalysis.h.
◆ CreateConstraints()
Implements Analysis.
Definition at line 10 of file EnthalpyAnalysis.cpp.
15 bool spcpresent =
false;
28 iomodel->
FindConstant(&heatcapacity,
"md.materials.heatcapacity");
29 iomodel->
FindConstant(&referencetemperature,
"md.constants.referencetemperature");
37 iomodel->
FetchData(&spcvector,&M,&N,
"md.thermal.spctemperature");
38 iomodel->
FetchData(&spcvectorstatic,&M,&N,
"md.thermal.spctemperature");
41 bool isdynamic =
false;
56 spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
57 spcvectorstatic[i*N+j] = NAN;
60 spcvector[i*N+j] = NAN;
61 spcvectorstatic[i*N+j] = heatcapacity*(spcvectorstatic[i*N+j]-referencetemperature);
65 spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
79 iomodel->
DeleteData(spcvector,
"md.thermal.spctemperature");
80 iomodel->
DeleteData(spcvectorstatic,
"md.thermal.spctemperature");
81 iomodel->
DeleteData(issurface,
"md.mesh.vertexonsurface");
82 xDelete<IssmDouble>(times);
83 xDelete<IssmDouble>(values);
◆ CreateLoads()
void EnthalpyAnalysis::CreateLoads |
( |
Loads * |
loads, |
|
|
IoModel * |
iomodel |
|
) |
| |
|
virtual |
◆ CreateNodes()
void EnthalpyAnalysis::CreateNodes |
( |
Nodes * |
nodes, |
|
|
IoModel * |
iomodel, |
|
|
bool |
isamr = false |
|
) |
| |
|
virtual |
◆ DofsPerNode()
int EnthalpyAnalysis::DofsPerNode |
( |
int ** |
doflist, |
|
|
int |
domaintype, |
|
|
int |
approximation |
|
) |
| |
|
virtual |
◆ UpdateElements()
void EnthalpyAnalysis::UpdateElements |
( |
Elements * |
elements, |
|
|
Inputs2 * |
inputs2, |
|
|
IoModel * |
iomodel, |
|
|
int |
analysis_counter, |
|
|
int |
analysis_type |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 101 of file EnthalpyAnalysis.cpp.
103 bool dakota_analysis,ismovingfront,isenthalpy;
104 int frictionlaw,basalforcing_model,materialstype;
105 int FrictionCoupling;
111 iomodel->
FindConstant(&isenthalpy,
"md.thermal.isenthalpy");
112 if(!isenthalpy)
return;
115 iomodel->
FetchData(3,
"md.initialization.temperature",
"md.initialization.waterfraction",
"md.initialization.pressure");
126 element->
Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
131 iomodel->
FindConstant(&dakota_analysis,
"md.qmu.isdakota");
132 iomodel->
FindConstant(&ismovingfront,
"md.transient.ismovingfront");
134 iomodel->
FindConstant(&materialstype,
"md.materials.type");
163 iomodel->
FindConstant(&basalforcing_model,
"md.basalforcings.model");
164 switch(basalforcing_model){
174 switch(materialstype){
196 iomodel->
FindConstant(&FrictionCoupling,
"md.friction.coupling");
200 if (FrictionCoupling==3){
202 else if(FrictionCoupling==4){
211 iomodel->
FindConstant(&FrictionCoupling,
"md.friction.coupling");
215 if (FrictionCoupling==3){
217 else if(FrictionCoupling==4){
227 iomodel->
FindConstant(&FrictionCoupling,
"md.friction.coupling");
242 iomodel->
FindConstant(&FrictionCoupling,
"md.friction.coupling");
247 if (FrictionCoupling==3){
249 else if(FrictionCoupling==4){
260 _error_(
"friction law not supported");
264 iomodel->
DeleteData(3,
"md.initialization.temperature",
"md.initialization.waterfraction",
"md.initialization.pressure");
◆ UpdateParameters()
void EnthalpyAnalysis::UpdateParameters |
( |
Parameters * |
parameters, |
|
|
IoModel * |
iomodel, |
|
|
int |
solution_enum, |
|
|
int |
analysis_enum |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 267 of file EnthalpyAnalysis.cpp.
270 char** requestedoutputs = NULL;
281 iomodel->
FindConstant(&requestedoutputs,&numoutputs,
"md.thermal.requested_outputs");
284 iomodel->
DeleteData(&requestedoutputs,numoutputs,
"md.thermal.requested_outputs");
297 if(frictionlaw==1 || frictionlaw==3 || frictionlaw==7){
◆ ApplyBasalConstraints()
void EnthalpyAnalysis::ApplyBasalConstraints |
( |
IssmDouble * |
serial_spc, |
|
|
Element * |
element |
|
) |
| |
|
static |
Definition at line 309 of file EnthalpyAnalysis.cpp.
316 bool isdynamicbasalspc;
324 if(!isdynamicbasalspc)
return;
334 for(
int i=0;i<numindices;i++){
340 node=element->
GetNode(indices[i]);
342 if(serial_spc[node->
Sid()]==1.){
352 xDelete<int>(indices);
◆ ComputeBasalMeltingrate() [1/2]
void EnthalpyAnalysis::ComputeBasalMeltingrate |
( |
FemModel * |
femmodel | ) |
|
|
static |
◆ ComputeBasalMeltingrate() [2/2]
void EnthalpyAnalysis::ComputeBasalMeltingrate |
( |
Element * |
element | ) |
|
|
static |
Definition at line 366 of file EnthalpyAnalysis.cpp.
380 int nodedown,nodeup,numnodes,numsegments;
382 IssmDouble vec_heatflux[dim],normal_base[dim],d1enthalpy[dim],d1pressure[dim];
383 IssmDouble basalfriction,alpha2,geothermalflux,heatflux;
389 int *pairindices = NULL;
420 element->
NormalBase(&normal_base[0],xyz_list_base);
422 IssmDouble* meltingrate_enthalpy = xNew<IssmDouble>(numsegments);
423 IssmDouble* heating = xNew<IssmDouble>(numsegments);
426 IssmDouble* enthalpies = xNew<IssmDouble>(numnodes);
427 IssmDouble* pressures = xNew<IssmDouble>(numnodes);
428 IssmDouble* watercolumns = xNew<IssmDouble>(numnodes);
429 IssmDouble* basalmeltingrates = xNew<IssmDouble>(numnodes);
438 for(is=0;is<numsegments;is++){
439 nodedown = pairindices[is*2+0];
440 nodeup = pairindices[is*2+1];
443 state=
GetThermalBasalCondition(element, enthalpies[nodedown], enthalpies[nodeup], pressures[nodedown], pressures[nodeup], watercolumns[nodedown], basalmeltingrates[nodedown]);
447 for(i=0;i<3;i++) vec_heatflux[i]=0.;
449 case 1:
case 2:
case 3:
454 for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
460 for(i=0;i<3;i++) vec_heatflux[i]=kappa_mix*beta*d1pressure[i];
463 _printf0_(
" unknown thermal basal state found!");
465 if(state==0) meltingrate_enthalpy[is]=0.;
469 for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
474 basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
477 heating[is]=(heatflux+basalfriction+geothermalflux);
478 meltingrate_enthalpy[is]=heating[is]/(latentheat*rho_ice);
483 for(is=0;is<numsegments;is++){
484 nodedown = pairindices[is*2+0];
485 nodeup = pairindices[is*2+1];
487 if(watercolumns[nodedown]+meltingrate_enthalpy[is]*dt<0.){
488 lambda = -watercolumns[nodedown]/(dt*meltingrate_enthalpy[is]);
_assert_(lambda>=0.);
_assert_(lambda<1.);
489 watercolumns[nodedown]=0.;
490 basalmeltingrates[nodedown]=lambda*meltingrate_enthalpy[is];
491 enthalpies[nodedown]+=(1.-lambda)*dt/yts*meltingrate_enthalpy[is]*latentheat*rho_ice;
494 basalmeltingrates[nodedown]=meltingrate_enthalpy[is];
495 watercolumns[nodedown]+=dt*meltingrate_enthalpy[is];
497 if(watercolumns[nodedown]>watercolumnupperlimit) watercolumns[nodedown]=watercolumnupperlimit;
500 basalmeltingrates[nodedown]=meltingrate_enthalpy[is];
501 if(watercolumns[nodedown]+meltingrate_enthalpy[is]<0.)
502 watercolumns[nodedown]=0.;
504 watercolumns[nodedown]+=meltingrate_enthalpy[is];
506 basalmeltingrates[nodedown]*=rho_water/rho_ice;
507 _assert_(watercolumns[nodedown]>=0.);
513 element->
AddInput2(enthalpy_enum,enthalpies,finite_element);
521 xDelete<int>(pairindices);
522 xDelete<IssmDouble>(enthalpies);
523 xDelete<IssmDouble>(pressures);
524 xDelete<IssmDouble>(watercolumns);
525 xDelete<IssmDouble>(basalmeltingrates);
526 xDelete<IssmDouble>(meltingrate_enthalpy);
527 xDelete<IssmDouble>(heating);
528 xDelete<IssmDouble>(xyz_list);
529 xDelete<IssmDouble>(xyz_list_base);
◆ Core()
void EnthalpyAnalysis::Core |
( |
FemModel * |
femmodel | ) |
|
|
virtual |
◆ CreateDVector()
◆ CreateJacobianMatrix()
◆ CreateKMatrix()
◆ CreateKMatrixVolume()
Definition at line 572 of file EnthalpyAnalysis.cpp.
582 IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
591 IssmDouble* basis = xNew<IssmDouble>(numnodes);
592 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
616 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
623 D_scalar=gauss->
weight*Jdet;
624 if(dt!=0.) D_scalar=D_scalar*dt;
627 for(
int i=0;i<numnodes;i++){
628 for(
int j=0;j<numnodes;j++){
629 Ke->
values[i*numnodes+j] += D_scalar*kappa/rho_ice*(
630 dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
639 for(
int i=0;i<numnodes;i++){
640 for(
int j=0;j<numnodes;j++){
641 Ke->
values[i*numnodes+j] += D_scalar*(
642 vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i] +vz*dbasis[2*numnodes+j]*basis[i]
649 D_scalar=gauss->
weight*Jdet;
650 for(
int i=0;i<numnodes;i++){
651 for(
int j=0;j<numnodes;j++){
652 Ke->
values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
655 D_scalar=D_scalar*dt;
659 if(stabilization==1){
661 vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
662 h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
663 K[0][0]=h/(2.*vel)*fabs(vx*vx); K[0][1]=h/(2.*vel)*fabs(vx*vy); K[0][2]=h/(2.*vel)*fabs(vx*vz);
664 K[1][0]=h/(2.*vel)*fabs(vy*vx); K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz);
665 K[2][0]=h/(2.*vel)*fabs(vz*vx); K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
666 for(
int i=0;i<3;i++)
for(
int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
668 for(
int i=0;i<numnodes;i++){
669 for(
int j=0;j<numnodes;j++){
670 Ke->
values[i*numnodes+j] += (
671 dbasis[0*numnodes+i] *(K[0][0]*dbasis[0*numnodes+j] + K[0][1]*dbasis[1*numnodes+j]+ K[0][2]*dbasis[2*numnodes+j]) +
672 dbasis[1*numnodes+i] *(K[1][0]*dbasis[0*numnodes+j] + K[1][1]*dbasis[1*numnodes+j]+ K[1][2]*dbasis[2*numnodes+j]) +
673 dbasis[2*numnodes+i] *(K[2][0]*dbasis[0*numnodes+j] + K[2][1]*dbasis[1*numnodes+j]+ K[2][2]*dbasis[2*numnodes+j])
679 else if(stabilization==2){
683 for(
int i=0;i<numnodes;i++){
684 for(
int j=0;j<numnodes;j++){
685 Ke->
values[i*numnodes+j]+=tau_parameter*D_scalar*
686 ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*
687 ((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
691 D_scalar=gauss->
weight*Jdet;
692 for(
int i=0;i<numnodes;i++){
693 for(
int j=0;j<numnodes;j++){
694 Ke->
values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
700 else if(stabilization==3){
704 tau_parameter_hor=tau_parameter_anisotropic[0];
705 tau_parameter_ver=tau_parameter_anisotropic[1];
706 for(
int i=0;i<numnodes;i++){
707 for(
int j=0;j<numnodes;j++){
708 Ke->
values[i*numnodes+j]+=D_scalar*
709 (sqrt(tau_parameter_hor)*(u-um)*dbasis[0*numnodes+i]+sqrt(tau_parameter_hor)*(v-vm)*dbasis[1*numnodes+i]+sqrt(tau_parameter_ver)*(w-wm)*dbasis[2*numnodes+i])*
710 (sqrt(tau_parameter_hor)*(u-um)*dbasis[0*numnodes+j]+sqrt(tau_parameter_hor)*(v-vm)*dbasis[1*numnodes+j]+sqrt(tau_parameter_ver)*(w-wm)*dbasis[2*numnodes+j]);
717 xDelete<IssmDouble>(xyz_list);
718 xDelete<IssmDouble>(basis);
719 xDelete<IssmDouble>(dbasis);
◆ CreateKMatrixShelf()
Definition at line 723 of file EnthalpyAnalysis.cpp.
740 IssmDouble* basis = xNew<IssmDouble>(numnodes);
754 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
760 D=gauss->
weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
761 if(reCast<bool,IssmDouble>(dt)) D=dt*D;
762 for(
int i=0;i<numnodes;i++)
for(
int j=0;j<numnodes;j++) Ke->
values[i*numnodes+j] += D*basis[i]*basis[j];
767 xDelete<IssmDouble>(basis);
768 xDelete<IssmDouble>(xyz_list_base);
◆ CreatePVector()
◆ CreatePVectorVolume()
Definition at line 788 of file EnthalpyAnalysis.cpp.
794 int i, stabilization;
797 IssmDouble enthalpypicard, d1enthalpypicard[3];
798 IssmDouble pressure, d1pressure[3], d2pressure;
800 IssmDouble kappa,tau_parameter,diameter,hx,hy,hz,kappa_w;
801 IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
803 IssmDouble scalar_def, scalar_sens ,scalar_transient;
813 IssmDouble* basis = xNew<IssmDouble>(numnodes);
814 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
831 Input2* enthalpy_input=NULL;
838 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
845 element->
ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
847 scalar_def=phi/rho_ice*Jdet*gauss->
weight;
848 if(dt!=0.) scalar_def=scalar_def*dt;
850 for(i=0;i<numnodes;i++) pe->
values[i]+=scalar_def*basis[i];
857 if(enthalpypicard>=Hpmp){
863 for(i=0;i<3;i++) d1H_d1P+=d1enthalpypicard[i]*d1pressure[i];
865 for(i=0;i<3;i++) d1P2+=pow(d1pressure[i],2.);
867 scalar_sens=-beta*((temperateiceconductivity - thermalconductivity)/latentheat*(d1H_d1P + beta*heatcapacity*d1P2))/rho_ice;
868 if(dt!=0.) scalar_sens=scalar_sens*dt;
869 for(i=0;i<numnodes;i++) pe->
values[i]+=scalar_sens*basis[i];
875 scalar_transient=enthalpy*Jdet*gauss->
weight;
876 for(i=0;i<numnodes;i++) pe->
values[i]+=scalar_transient*basis[i];
880 if(stabilization==2){
889 for(i=0;i<numnodes;i++) pe->
values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
892 for(i=0;i<numnodes;i++) pe->
values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
896 else if(stabilization==3){
904 tau_parameter_hor=tau_parameter_anisotropic[0];
905 tau_parameter_ver=tau_parameter_anisotropic[1];
907 for(i=0;i<numnodes;i++) pe->
values[i]+=scalar_def*(tau_parameter_hor*u*dbasis[0*numnodes+i]+tau_parameter_hor*v*dbasis[1*numnodes+i]+tau_parameter_ver*w*dbasis[2*numnodes+i]);
912 xDelete<IssmDouble>(basis);
913 xDelete<IssmDouble>(dbasis);
914 xDelete<IssmDouble>(xyz_list);
◆ CreatePVectorSheet()
Definition at line 919 of file EnthalpyAnalysis.cpp.
927 bool converged, isdynamicbasalspc;
931 IssmDouble enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
933 IssmDouble alpha2,basalfriction,geothermalflux,heatflux;
941 IssmDouble* basis = xNew<IssmDouble>(numnodes);
966 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
973 if(isdynamicbasalspc){
986 case 0:
case 1:
case 2:
case 3:
996 basalfriction=alpha2*(vx*vx+vy*vy+vz*vz);
997 heatflux=(basalfriction+geothermalflux)/(rho_ice);
998 scalar=gauss->
weight*Jdet*heatflux;
999 if(dt!=0.) scalar=dt*scalar;
1000 for(i=0;i<numnodes;i++)
1001 pe->
values[i]+=scalar*basis[i];
1005 for(i=0;i<numnodes;i++)
1009 _printf0_(
" unknown thermal basal state found!");
1017 xDelete<IssmDouble>(basis);
1018 xDelete<IssmDouble>(xyz_list_base);
◆ CreatePVectorShelf()
Definition at line 1022 of file EnthalpyAnalysis.cpp.
1030 IssmDouble Hpmp,dt,Jdet,scalar_ocean,pressure;
1038 IssmDouble* basis = xNew<IssmDouble>(numnodes);
1053 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
1062 scalar_ocean=gauss->
weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*Hpmp/(heatcapacity*rho_ice);
1063 if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
1065 for(
int i=0;i<numnodes;i++) pe->
values[i]+=scalar_ocean*basis[i];
1070 xDelete<IssmDouble>(basis);
1071 xDelete<IssmDouble>(xyz_list_base);
◆ DrainWaterfraction()
void EnthalpyAnalysis::DrainWaterfraction |
( |
FemModel * |
femmodel | ) |
|
|
static |
◆ ComputeWaterfractionDrainage()
void EnthalpyAnalysis::ComputeWaterfractionDrainage |
( |
FemModel * |
femmodel | ) |
|
|
static |
Definition at line 1080 of file EnthalpyAnalysis.cpp.
1091 IssmDouble* waterfractions= xNew<IssmDouble>(numnodes);
1092 IssmDouble* drainage= xNew<IssmDouble>(numnodes);
1095 for(k=0; k<numnodes;k++){
1101 xDelete<IssmDouble>(waterfractions);
1102 xDelete<IssmDouble>(drainage);
◆ DrainageUpdateWatercolumn()
void EnthalpyAnalysis::DrainageUpdateWatercolumn |
( |
FemModel * |
femmodel | ) |
|
|
static |
Definition at line 1105 of file EnthalpyAnalysis.cpp.
1107 int i,k,numnodes, numbasalnodes;
1109 int* basalnodeindices=NULL;
1124 IssmDouble* drainage_int= xNew<IssmDouble>(numnodes);
1125 IssmDouble* thicknesses= xNew<IssmDouble>(numnodes);
1129 for(k=0;k<numnodes;k++){
1130 drainage_int[k]*=thicknesses[k];
1135 xDelete<IssmDouble>(drainage_int);
1136 xDelete<IssmDouble>(thicknesses);
1147 IssmDouble* watercolumn= xNew<IssmDouble>(numnodes);
1148 IssmDouble* drainage_int= xNew<IssmDouble>(numnodes);
1153 for(k=0;k<numbasalnodes;k++){
1154 watercolumn[basalnodeindices[k]]+=dt*drainage_int[basalnodeindices[k]];
1159 xDelete<IssmDouble>(watercolumn);
1160 xDelete<IssmDouble>(drainage_int);
1161 xDelete<int>(basalnodeindices);
◆ DrainageUpdateEnthalpy()
void EnthalpyAnalysis::DrainageUpdateEnthalpy |
( |
FemModel * |
femmodel | ) |
|
|
static |
Definition at line 1164 of file EnthalpyAnalysis.cpp.
1174 IssmDouble* enthalpies= xNew<IssmDouble>(numnodes);
1175 IssmDouble* pressures= xNew<IssmDouble>(numnodes);
1176 IssmDouble* temperatures= xNew<IssmDouble>(numnodes);
1177 IssmDouble* waterfractions= xNew<IssmDouble>(numnodes);
1178 IssmDouble* drainage= xNew<IssmDouble>(numnodes);
1185 for(k=0;k<numnodes;k++){
1187 waterfractions[k]-=drainage[k];
1189 waterfractions[k]-=dt*drainage[k];
1191 element->
ThermalToEnthalpy(&enthalpies[k], temperatures[k], waterfractions[k], pressures[k]);
1197 xDelete<IssmDouble>(enthalpies);
1198 xDelete<IssmDouble>(pressures);
1199 xDelete<IssmDouble>(temperatures);
1200 xDelete<IssmDouble>(waterfractions);
1201 xDelete<IssmDouble>(drainage);
◆ EnthalpyDiffusionParameter()
Definition at line 1204 of file EnthalpyAnalysis.cpp.
1211 return thermalconductivity/heatcapacity;
1213 return temperateiceconductivity/heatcapacity;
◆ EnthalpyDiffusionParameterVolume()
IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume |
( |
Element * |
element, |
|
|
int |
enthalpy_enum |
|
) |
| |
|
static |
Definition at line 1215 of file EnthalpyAnalysis.cpp.
1224 int effectiveconductivity_averaging;
1225 IssmDouble* pressures = xNew<IssmDouble>(numvertices);
1226 IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
1227 IssmDouble* PIE = xNew<IssmDouble>(numvertices);
1228 IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
1233 for(iv=0;iv<numvertices;iv++){
1235 dHpmp[iv] = enthalpies[iv]-PIE[iv];
1238 bool allequalsign =
true;
1240 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0.));
1243 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0.));
1254 for(iv=0; iv<numvertices;iv++){
1255 if(enthalpies[iv]<PIE[iv])
1256 Hc+=(PIE[iv]-enthalpies[iv]);
1258 Ht+=(enthalpies[iv]-PIE[iv]);
1261 lambda = Hc/(Hc+Ht);
1265 if(effectiveconductivity_averaging==0){
1267 kappa=kappa_c*lambda+(1.-lambda)*kappa_t;
1269 else if(effectiveconductivity_averaging==1){
1271 kappa=kappa_c*kappa_t/(lambda*kappa_t+(1.-lambda)*kappa_c);
1273 else if(effectiveconductivity_averaging==2){
1275 kappa=pow(kappa_c,lambda)*pow(kappa_t,1.-lambda);
1278 _error_(
"effectiveconductivity_averaging not supported yet");
1283 xDelete<IssmDouble>(PIE);
1284 xDelete<IssmDouble>(dHpmp);
1285 xDelete<IssmDouble>(pressures);
1286 xDelete<IssmDouble>(enthalpies);
◆ GetBasalConstraints()
◆ GetBasalConstraintsSteadystate()
Definition at line 1307 of file EnthalpyAnalysis.cpp.
1317 int numindices, numindicesup, state;
1318 int *indices = NULL, *indicesup = NULL;
1319 IssmDouble enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
1334 for(
int i=0;i<numindices;i++){
1368 _printf0_(
" unknown thermal basal state found!");
1373 xDelete<int>(indices);
1374 xDelete<int>(indicesup);
◆ GetBasalConstraintsTransient()
Definition at line 1378 of file EnthalpyAnalysis.cpp.
1388 int numindices, numindicesup, state;
1389 int *indices = NULL, *indicesup = NULL;
1390 IssmDouble enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
1406 for(
int i=0;i<numindices;i++){
1441 _printf0_(
" unknown thermal basal state found!");
1447 xDelete<int>(indices);
1448 xDelete<int>(indicesup);
◆ GetBConduct()
Definition at line 1452 of file EnthalpyAnalysis.cpp.
1468 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
1472 for(
int i=0;i<numnodes;i++){
1473 B[numnodes*0+i] = dbasis[0*numnodes+i];
1474 B[numnodes*1+i] = dbasis[1*numnodes+i];
1475 B[numnodes*2+i] = dbasis[2*numnodes+i];
1479 xDelete<IssmDouble>(dbasis);
◆ GetSolutionFromInputs()
◆ GetThermalBasalCondition()
Definition at line 1484 of file EnthalpyAnalysis.cpp.
1490 if(!(element->
IsOnBase()))
return -1;
1500 if(watercolumn<=0.) state=0;
1505 if((dt==0.) && (meltingrate<0.)) state=2;
◆ GetWetIceConductivity()
◆ GradientJ()
void EnthalpyAnalysis::GradientJ |
( |
Vector< IssmDouble > * |
gradient, |
|
|
Element * |
element, |
|
|
int |
control_type, |
|
|
int |
control_index |
|
) |
| |
|
virtual |
◆ InputUpdateFromSolution()
void EnthalpyAnalysis::InputUpdateFromSolution |
( |
IssmDouble * |
solution, |
|
|
Element * |
element |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 1526 of file EnthalpyAnalysis.cpp.
1530 IssmDouble B_average,s_average,T_average=0.,P_average=0.;
1531 int *doflist = NULL;
1539 IssmDouble* values = xNew<IssmDouble>(numnodes);
1540 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
1541 IssmDouble* surface = xNew<IssmDouble>(numnodes);
1543 IssmDouble* temperature = xNew<IssmDouble>(numnodes);
1544 IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
1547 for(i=0;i<numnodes;i++){
1548 values[i]=solution[doflist[i]];
1551 if(xIsNan<IssmDouble>(values[i]))
_error_(
"NaN found in solution vector");
1552 if(xIsInf<IssmDouble>(values[i]))
_error_(
"Inf found in solution vector");
1560 for(i=0;i<numnodes;i++){
1561 element->
EnthalpyToThermal(&temperature[i],&waterfraction[i],values[i],pressure[i]);
1562 if(waterfraction[i]<0.)
_error_(
"Negative water fraction found in solution vector");
1571 for(i=0;i<numnodes;i++) n[i]=3.;
1581 switch(rheology_law){
1586 for(i=0;i<numnodes;i++) B[i]=
BuddJacka(temperature[i]);
1590 for(i=0;i<numnodes;i++) B[i]=
Cuffey(temperature[i]);
1594 for(i=0;i<numnodes;i++) B[i]=
CuffeyTemperate(temperature[i], waterfraction[i],n[i]);
1598 for(i=0;i<numnodes;i++) B[i]=
Paterson(temperature[i]);
1602 for(i=0;i<numnodes;i++) B[i]=
NyeH2O(values[i]);
1606 for(i=0;i<numnodes;i++) B[i]=
NyeCO2(values[i]);
1611 for(i=0;i<numnodes;i++) B[i]=
Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],n[i]);
1622 xDelete<IssmDouble>(n);
1629 xDelete<IssmDouble>(values);
1630 xDelete<IssmDouble>(pressure);
1631 xDelete<IssmDouble>(surface);
1632 xDelete<IssmDouble>(B);
1633 xDelete<IssmDouble>(temperature);
1634 xDelete<IssmDouble>(waterfraction);
1635 xDelete<IssmDouble>(xyz_list);
1636 xDelete<int>(doflist);
◆ PostProcessing()
void EnthalpyAnalysis::PostProcessing |
( |
FemModel * |
femmodel | ) |
|
|
static |
Definition at line 1638 of file EnthalpyAnalysis.cpp.
1641 bool computebasalmeltingrates=
true;
1642 bool isdrainicecolumn;
1648 if(isdrainicecolumn){
1651 if(computebasalmeltingrates){
◆ PureIceEnthalpy()
◆ TMeltingPoint()
◆ UpdateBasalConstraints()
void EnthalpyAnalysis::UpdateBasalConstraints |
( |
FemModel * |
femmodel | ) |
|
|
static |
◆ UpdateConstraints()
void EnthalpyAnalysis::UpdateConstraints |
( |
FemModel * |
femmodel | ) |
|
|
virtual |
The documentation for this class was generated from the following files:
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
@ MaterialsThermalconductivityEnum
static void DrainageUpdateWatercolumn(FemModel *femmodel)
@ ThermalWatercolumnUpperlimitEnum
@ FrictionCoefficientEnum
void depthaverage_core(FemModel *femmodel)
virtual int GetElementType(void)=0
IssmDouble Cuffey(IssmDouble temperature)
IssmDouble Paterson(IssmDouble temperature)
void GetInputListOnNodes(IssmDouble *pvalue, int enumtype)
@ ThermalStabilizationEnum
void extrudefrombase_core(FemModel *femmodel)
@ InputToDepthaverageOutEnum
virtual int GetNumberOfNodes(void)=0
@ ThermalIsdrainicecolumnEnum
@ ThermalRequestedOutputsEnum
static void UpdateBasalConstraints(FemModel *femmodel)
#define _printf0_(StreamArgs)
@ FrictionEffectivePressureEnum
void FindParam(bool *pvalue, int paramenum)
virtual Gauss * NewGaussBase(int order)=0
@ MaterialsRheologyEcEnum
@ MaterialsRhoFreshwaterEnum
@ MaterialsMeltingpointEnum
@ TimesteppingTimeStepEnum
@ ThermalNumRequestedOutputsEnum
@ ThermalIsdynamicbasalspcEnum
static void ComputeBasalMeltingrate(FemModel *femmodel)
@ MaterialsLatentheatEnum
virtual void GaussNode(int finitelement, int iv)=0
IssmDouble PureIceEnthalpy(IssmDouble pressure)
virtual void JacobianDeterminantBase(IssmDouble *Jdet, IssmDouble *xyz_list_base, Gauss *gauss)=0
virtual void NormalBase(IssmDouble *normal, IssmDouble *xyz_list)=0
void BasalNodeIndices(int *pnumindices, int **pindices, int finiteelement)
ElementMatrix * CreateKMatrixVolume(Element *element)
@ ConstantsReferencetemperatureEnum
static int GetThermalBasalCondition(Element *element, IssmDouble enthalpy, IssmDouble enthalpy_up, IssmDouble pressure, IssmDouble pressure_up, IssmDouble watercolumn, IssmDouble meltingrate)
void AddObject(Param *newparam)
static void PostProcessing(FemModel *femmodel)
ElementVector * CreatePVectorShelf(Element *element)
void ThermalToEnthalpy(IssmDouble *penthalpy, IssmDouble temperature, IssmDouble waterfraction, IssmDouble pressure)
ElementVector * CreatePVectorVolume(Element *element)
void IoModelToDynamicConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
ElementVector * CreatePVectorSheet(Element *element)
virtual Input2 * GetInput2(int inputenum)=0
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
static IssmDouble GetWetIceConductivity(Element *element, IssmDouble enthalpy, IssmDouble pressure)
IssmDouble NyeCO2(IssmDouble temperature)
virtual void ElementSizes(IssmDouble *phx, IssmDouble *phy, IssmDouble *phz)=0
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
void SetParam(bool boolean, int enum_type)
void DeleteData(int num,...)
static void DrainWaterfraction(FemModel *femmodel)
Param * CopyConstantObject(const char *constant_name, int param_enum)
@ BasalforcingsGroundediceMeltingRateEnum
IssmDouble BuddJacka(IssmDouble temperature)
void SurfaceNodeIndices(int *pnumindices, int **pindices, int finiteelement)
void UpdateConstraintsx(void)
static void GetBasalConstraints(Vector< IssmDouble > *vec_spc, Element *element)
IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.)
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
virtual Gauss * NewGauss(void)=0
@ FrictionPressureAdjustedTemperatureEnum
@ MaterialsEffectiveconductivityAveragingEnum
virtual void VerticalSegmentIndicesBase(int **pindices, int *pnumseg)=0
@ MaterialsThermalExchangeVelocityEnum
const char * EnumToStringx(int enum_in)
static IssmDouble TMeltingPoint(Element *element, IssmDouble pressure)
@ MantlePlumeGeothermalFluxEnum
void FindConstant(bool *pvalue, const char *constant_name)
void solutionsequence_thermal_nonlinear(FemModel *femmodel)
IssmDouble NyeH2O(IssmDouble temperature)
virtual void StabilizationParameterAnisotropic(IssmDouble *tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa)=0
void GetVerticesCoordinates(IssmDouble **xyz_list)
@ MaterialsHeatcapacityEnum
Node * GetNode(int nodeindex)
ElementMatrix * CreateKMatrixShelf(Element *element)
static void ComputeWaterfractionDrainage(FemModel *femmodel)
@ MaterialsRheologyLawEnum
void FetchData(bool *pboolean, const char *data_name)
@ MaterialsRheologyEsEnum
@ MaterialsRhoSeawaterEnum
static void GetBasalConstraintsSteadystate(Vector< IssmDouble > *vec_spc, Element *element)
static void DrainageUpdateEnthalpy(FemModel *femmodel)
IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat)
static IssmDouble EnthalpyDiffusionParameter(Element *element, IssmDouble enthalpy, IssmDouble pressure)
virtual int ObjectEnum()=0
void GetDofListLocal(int **pdoflist, int approximation_enum, int setenum)
virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0
void GetSolutionFromInputsOneDof(Vector< IssmDouble > *solution, int solutionenum)
@ InputToDepthaverageInEnum
virtual void BasalNodeIndices(int *pnumindices, int **pindices, int finiteelement)
virtual void GetVerticesCoordinatesBase(IssmDouble **xyz_list)=0
#define _error_(StreamArgs)
void SetActiveNodesLSMx(FemModel *femmodel)
virtual int begin(void)=0
bool VerboseSolution(void)
Object * GetObjectByOffset(int offset)
@ MeshVertexonsurfaceEnum
IssmDouble Arrhenius(IssmDouble temperature, IssmDouble depth, IssmDouble n)
void GaussNode(int finitelement, int iv)
@ WaterfractionDrainageEnum
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
@ FrictionEffectivePressureLimitEnum
virtual void GaussPoint(int ig)=0
virtual void Update(Inputs2 *inputs2, int index, IoModel *iomodel, int analysis_counter, int analysis_type, int finite_element)=0
void SetCurrentConfiguration(int configuration_type)
void GetInputValue(bool *pvalue, int enum_type)
void FindParam(bool *pinteger, int enum_type)
void GetAlpha2(IssmDouble *palpha2, Gauss *gauss)
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
@ FrictionCoefficientcoulombEnum
static IssmDouble PureIceEnthalpy(Element *element, IssmDouble pressure)
static void GetBasalConstraintsTransient(Vector< IssmDouble > *vec_spc, Element *element)
@ WaterfractionDrainageIntegratedEnum
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
void ApplyConstraint(int dof, IssmDouble value)
@ BasalforcingsGeothermalfluxEnum
virtual void ViscousHeating(IssmDouble *pphi, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input, Input2 *vz_input)
static IssmDouble EnthalpyDiffusionParameterVolume(Element *element, int enthalpy_enum)
virtual Gauss * NewGaussTop(int order)=0
doubletype * ToMPISerial(void)
@ MaterialsTemperateiceconductivityEnum
virtual int GetNumberOfVertices(void)=0
void GetInputListOnVertices(IssmDouble *pvalue, int enumtype)
static void ApplyBasalConstraints(IssmDouble *serial_spc, Element *element)
IssmDouble CuffeyTemperate(IssmDouble temperature, IssmDouble waterfraction, IssmDouble stressexp)
void SetValue(int dof, doubletype value, InsMode mode)
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
@ MaterialsMixedLayerCapacityEnum
void EnthalpyToThermal(IssmDouble *ptemperature, IssmDouble *pwaterfraction, IssmDouble enthalpy, IssmDouble pressure)
virtual IssmDouble MinEdgeLength(IssmDouble *xyz_list)=0
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)