source: issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp@ 26047

Last change on this file since 26047 was 26047, checked in by Eric.Larour, 4 years ago

CHG: huge commit on solid earth capability rewrite. Complete cleanup of the sea level core.
New mass transport capabilities for ocean and tws. No more giacore. GiaIvins folded into
the sea level core. Debugging of Materials.

File size: 67.9 KB
RevLine 
[16534]1#include "./EnthalpyAnalysis.h"
2#include "../toolkits/toolkits.h"
3#include "../classes/classes.h"
4#include "../shared/shared.h"
5#include "../modules/modules.h"
[18591]6#include "../solutionsequences/solutionsequences.h"
[21721]7#include "../cores/cores.h"
[16534]8
9/*Model processing*/
[18930]10void EnthalpyAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
[16604]11
[18930]12 /*Intermediary*/
13 int count;
14 int M,N;
15 bool spcpresent = false;
[21546]16 int finiteelement;
[18930]17 IssmDouble heatcapacity;
18 IssmDouble referencetemperature;
[16604]19
[18930]20 /*Output*/
21 IssmDouble *spcvector = NULL;
[22484]22 IssmDouble *spcvectorstatic = NULL;
[18930]23 IssmDouble* times=NULL;
24 IssmDouble* values=NULL;
[22484]25 IssmDouble* issurface = NULL;
[16604]26
[18930]27 /*Fetch parameters: */
[20690]28 iomodel->FindConstant(&heatcapacity,"md.materials.heatcapacity");
29 iomodel->FindConstant(&referencetemperature,"md.constants.referencetemperature");
[21546]30 iomodel->FindConstant(&finiteelement,"md.thermal.fe");
[18732]31
[18930]32 /*return if 2d mesh*/
33 if(iomodel->domaintype==Domain2DhorizontalEnum) return;
34
35 /*Fetch data: */
[22484]36 iomodel->FetchData(&issurface,&M,&N,"md.mesh.vertexonsurface"); _assert_(N>0); _assert_(M==iomodel->numberofvertices);
[20690]37 iomodel->FetchData(&spcvector,&M,&N,"md.thermal.spctemperature");
[22484]38 iomodel->FetchData(&spcvectorstatic,&M,&N,"md.thermal.spctemperature");
[18930]39
[21548]40 /*Specific case for PDD, we want the constaints to be updated by the PDD scheme itself*/
41 bool isdynamic = false;
42 if (iomodel->solution_enum==TransientSolutionEnum){
43 int smb_model;
44 iomodel->FindConstant(&smb_model,"md.smb.model");
[23317]45 if(smb_model==SMBpddEnum) isdynamic=true;
46 if(smb_model==SMBd18opddEnum) isdynamic=true;
47 if(smb_model==SMBpddSicopolisEnum) isdynamic=true;
[21548]48 }
49
[22484]50 /*Convert spcs from temperatures to enthalpy*/
51 _assert_(N>0); _assert_(M>=iomodel->numberofvertices);
52 for(int i=0;i<iomodel->numberofvertices;i++){
53 for(int j=0;j<N;j++){
54 if (isdynamic){
55 if (issurface[i]==1){
56 spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
[22485]57 spcvectorstatic[i*N+j] = NAN;
[22484]58 }
59 else{
[22485]60 spcvector[i*N+j] = NAN;
[22484]61 spcvectorstatic[i*N+j] = heatcapacity*(spcvectorstatic[i*N+j]-referencetemperature);
62 }
63 }
64 else{
65 spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
66 }
67 }
68 }
69
[21548]70 if(isdynamic){
[22484]71 IoModelToDynamicConstraintsx(constraints,iomodel,spcvector,iomodel->numberofvertices,1,EnthalpyAnalysisEnum,finiteelement);
72 IoModelToConstraintsx(constraints,iomodel,spcvectorstatic,M,N,EnthalpyAnalysisEnum,finiteelement);
[21548]73 }
74 else{
75 IoModelToConstraintsx(constraints,iomodel,spcvector,M,N,EnthalpyAnalysisEnum,finiteelement);
76 }
77
[18930]78 /*Free ressources:*/
[20690]79 iomodel->DeleteData(spcvector,"md.thermal.spctemperature");
[22484]80 iomodel->DeleteData(spcvectorstatic,"md.thermal.spctemperature");
81 iomodel->DeleteData(issurface,"md.mesh.vertexonsurface");
[18930]82 xDelete<IssmDouble>(times);
83 xDelete<IssmDouble>(values);
[16539]84}/*}}}*/
[18930]85void EnthalpyAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
86
87 /*No loads */
88}/*}}}*/
[23585]89void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
[18930]90
[21542]91 int finiteelement;
92 iomodel->FindConstant(&finiteelement,"md.thermal.fe");
93
[20690]94 if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
[21542]95 ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,finiteelement);
[20690]96 iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
[18930]97}/*}}}*/
98int EnthalpyAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
99 return 1;
100}/*}}}*/
[25379]101void EnthalpyAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
[16539]102
[20459]103 bool dakota_analysis,ismovingfront,isenthalpy;
[21382]104 int frictionlaw,basalforcing_model,materialstype;
[19161]105 int FrictionCoupling;
[23066]106
[16539]107 /*Now, is the model 3d? otherwise, do nothing: */
[17700]108 if(iomodel->domaintype==Domain2DhorizontalEnum)return;
[16539]109
110 /*Is enthalpy requested?*/
[20690]111 iomodel->FindConstant(&isenthalpy,"md.thermal.isenthalpy");
[16539]112 if(!isenthalpy) return;
113
114 /*Fetch data needed: */
[20690]115 iomodel->FetchData(3,"md.initialization.temperature","md.initialization.waterfraction","md.initialization.pressure");
[16539]116
[21542]117 /*Finite element type*/
118 int finiteelement;
119 iomodel->FindConstant(&finiteelement,"md.thermal.fe");
120
[16539]121 /*Update elements: */
122 int counter=0;
123 for(int i=0;i<iomodel->numberofelements;i++){
124 if(iomodel->my_elements[i]){
125 Element* element=(Element*)elements->GetObjectByOffset(counter);
[25379]126 element->Update(inputs,i,iomodel,analysis_counter,analysis_type,finiteelement);
[16539]127 counter++;
128 }
129 }
130
[20690]131 iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota");
132 iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
133 iomodel->FindConstant(&frictionlaw,"md.friction.law");
[21382]134 iomodel->FindConstant(&materialstype,"md.materials.type");
[16539]135
[25379]136 iomodel->FetchDataToInput(inputs,elements,"md.geometry.thickness",ThicknessEnum);
137 iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
[26047]138 iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0);
[25379]139 iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
140 iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
141 iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
[17886]142 if(iomodel->domaintype!=Domain2DhorizontalEnum){
[25379]143 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
144 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
[17886]145 }
[25379]146 iomodel->FetchDataToInput(inputs,elements,"md.initialization.pressure",PressureEnum);
147 iomodel->FetchDataToInput(inputs,elements,"md.initialization.temperature",TemperatureEnum);
148 iomodel->FetchDataToInput(inputs,elements,"md.initialization.waterfraction",WaterfractionEnum);
149 iomodel->FetchDataToInput(inputs,elements,"md.initialization.enthalpy",EnthalpyEnum);
150 iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",WatercolumnEnum);
151 iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
152 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
153 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
154 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vz",VzEnum);
155 InputUpdateFromConstantx(inputs,elements,0.,VxMeshEnum);
156 InputUpdateFromConstantx(inputs,elements,0.,VyMeshEnum);
157 InputUpdateFromConstantx(inputs,elements,0.,VzMeshEnum);
[20459]158 if(ismovingfront){
[25379]159 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); // required for updating active nodes
[17434]160 }
[20020]161
162 /*Basal forcings variables*/
[20690]163 iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
[20020]164 switch(basalforcing_model){
165 case MantlePlumeGeothermalFluxEnum:
166 break;
167 default:
[25379]168 iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
[20020]169 break;
170 }
[21382]171
172 /*Rheology type*/
[25379]173 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
[21382]174 switch(materialstype){
[21389]175 case MatenhancediceEnum:
[25379]176 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
177 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_E",MaterialsRheologyEEnum);
[21389]178 break;
[21382]179 case MatdamageiceEnum:
[25379]180 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
[21382]181 break;
182 case MatestarEnum:
[25379]183 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
184 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
[21382]185 break;
186 case MaticeEnum:
[25379]187 iomodel->FetchDataToInput(inputs,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
[21382]188 break;
189 default:
190 _error_("not supported");
191 }
[23066]192
[17952]193 /*Friction law variables*/
194 switch(frictionlaw){
195 case 1:
[21740]196 iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
[25379]197 iomodel->FetchDataToInput(inputs,elements,"md.friction.coefficient",FrictionCoefficientEnum);
198 iomodel->FetchDataToInput(inputs,elements,"md.friction.p",FrictionPEnum);
199 iomodel->FetchDataToInput(inputs,elements,"md.friction.q",FrictionQEnum);
[23620]200 if (FrictionCoupling==3){
[25379]201 iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
[23620]202 else if(FrictionCoupling==4){
[25379]203 iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",EffectivePressureEnum);
[21740]204 }
[17952]205 break;
206 case 2:
[25379]207 iomodel->FetchDataToInput(inputs,elements,"md.friction.C",FrictionCEnum);
208 iomodel->FetchDataToInput(inputs,elements,"md.friction.m",FrictionMEnum);
[17952]209 break;
[18778]210 case 3:
[20690]211 iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
[25379]212 iomodel->FetchDataToInput(inputs,elements,"md.friction.C",FrictionCEnum);
213 iomodel->FetchDataToInput(inputs,elements,"md.friction.As",FrictionAsEnum);
214 iomodel->FetchDataToInput(inputs,elements,"md.friction.q",FrictionQEnum);
[23620]215 if (FrictionCoupling==3){
[25379]216 iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
[23620]217 else if(FrictionCoupling==4){
[25379]218 iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",EffectivePressureEnum);
[19161]219 }
[18778]220 break;
[18732]221 case 4:
[25379]222 iomodel->FetchDataToInput(inputs,elements,"md.friction.coefficient",FrictionCoefficientEnum);
223 iomodel->FetchDataToInput(inputs,elements,"md.friction.p",FrictionPEnum);
224 iomodel->FetchDataToInput(inputs,elements,"md.friction.q",FrictionQEnum);
225 iomodel->FetchDataToInput(inputs,elements,"md.initialization.pressure",PressureEnum);
226 iomodel->FetchDataToInput(inputs,elements,"md.initialization.temperature",TemperatureEnum);
[23475]227 iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
[18732]228 break;
[18772]229 case 5:
[25379]230 iomodel->FetchDataToInput(inputs,elements,"md.friction.coefficient",FrictionCoefficientEnum);
231 iomodel->FetchDataToInput(inputs,elements,"md.friction.p",FrictionPEnum);
232 iomodel->FetchDataToInput(inputs,elements,"md.friction.q",FrictionQEnum);
233 iomodel->FetchDataToInput(inputs,elements,"md.friction.water_layer",FrictionWaterLayerEnum);
[18772]234 break;
[18804]235 case 6:
[25379]236 iomodel->FetchDataToInput(inputs,elements,"md.friction.C",FrictionCEnum);
237 iomodel->FetchDataToInput(inputs,elements,"md.friction.m",FrictionMEnum);
238 iomodel->FetchDataToInput(inputs,elements,"md.initialization.pressure",PressureEnum);
239 iomodel->FetchDataToInput(inputs,elements,"md.initialization.temperature",TemperatureEnum);
[18804]240 break;
[22048]241 case 7:
242 iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
[25379]243 iomodel->FetchDataToInput(inputs,elements,"md.friction.coefficient",FrictionCoefficientEnum);
244 iomodel->FetchDataToInput(inputs,elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum);
245 iomodel->FetchDataToInput(inputs,elements,"md.friction.p",FrictionPEnum);
246 iomodel->FetchDataToInput(inputs,elements,"md.friction.q",FrictionQEnum);
[23620]247 if (FrictionCoupling==3){
[25379]248 iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
[23620]249 else if(FrictionCoupling==4){
[25379]250 iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",EffectivePressureEnum);
[22048]251 }
252 break;
[21556]253 case 9:
[25379]254 iomodel->FetchDataToInput(inputs,elements,"md.friction.coefficient",FrictionCoefficientEnum);
255 iomodel->FetchDataToInput(inputs,elements,"md.friction.pressure_adjusted_temperature",FrictionPressureAdjustedTemperatureEnum);
256 InputUpdateFromConstantx(inputs,elements,1.,FrictionPEnum);
257 InputUpdateFromConstantx(inputs,elements,1.,FrictionQEnum);
[21556]258 break;
[17952]259 default:
[21873]260 _error_("friction law not supported");
[17952]261 }
[20690]262
[16539]263 /*Free data: */
[20690]264 iomodel->DeleteData(3,"md.initialization.temperature","md.initialization.waterfraction","md.initialization.pressure");
265
[16539]266}/*}}}*/
[18930]267void EnthalpyAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
[16539]268
[18930]269 int numoutputs;
270 char** requestedoutputs = NULL;
271
[20690]272 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.stabilization",ThermalStabilizationEnum));
273 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.maxiter",ThermalMaxiterEnum));
274 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.reltol",ThermalReltolEnum));
275 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isenthalpy",ThermalIsenthalpyEnum));
276 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdynamicbasalspc",ThermalIsdynamicbasalspcEnum));
[24140]277 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdrainicecolumn",ThermalIsdrainicecolumnEnum));
278 parameters->AddObject(iomodel->CopyConstantObject("md.thermal.watercolumn_upperlimit",ThermalWatercolumnUpperlimitEnum));
[20690]279 parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
[18930]280
[20690]281 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.thermal.requested_outputs");
[18930]282 parameters->AddObject(new IntParam(ThermalNumRequestedOutputsEnum,numoutputs));
283 if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs));
[20690]284 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.thermal.requested_outputs");
[18930]285
286 /*Deal with friction parameters*/
287 int frictionlaw;
[20690]288 iomodel->FindConstant(&frictionlaw,"md.friction.law");
[24666]289 if(frictionlaw==6){
[21776]290 parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
291 }
[24666]292 if(frictionlaw==4){
293 parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
[21776]294 parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
[24666]295 parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
[21776]296 }
[24666]297 if(frictionlaw==1 || frictionlaw==3 || frictionlaw==7){
298 parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
299 parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
300 }
[21776]301 if(frictionlaw==9){
302 parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
[25539]303 parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
[21778]304 parameters->AddObject(new IntParam(FrictionCouplingEnum,0));
[21776]305 }
[16539]306}/*}}}*/
307
[18930]308/*Finite Element Analysis*/
309void EnthalpyAnalysis::ApplyBasalConstraints(IssmDouble* serial_spc,Element* element){/*{{{*/
[16539]310
[21481]311 /* Do not check if ice in element, this may lead to inconsistencies between cpu partitions */
[20213]312 /* Only update constraints at the base. */
313 if(!(element->IsOnBase())) return;
[16539]314
[18930]315 /*Intermediary*/
316 bool isdynamicbasalspc;
317 int numindices;
318 int *indices = NULL;
319 Node* node = NULL;
320 IssmDouble pressure;
[16539]321
[18930]322 /*Check wether dynamic basal boundary conditions are activated */
323 element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
324 if(!isdynamicbasalspc) return;
[16539]325
[18930]326 /*Get parameters and inputs: */
[25379]327 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
[16539]328
[18930]329 /*Fetch indices of basal & surface nodes for this finite element*/
330 Penta *penta = (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
331 penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
[16539]332
[18930]333 GaussPenta* gauss=new GaussPenta();
334 for(int i=0;i<numindices;i++){
335 gauss->GaussNode(element->GetElementType(),indices[i]);
[16539]336
[18930]337 pressure_input->GetInputValue(&pressure,gauss);
[16539]338
[18930]339 /*apply or release spc*/
340 node=element->GetNode(indices[i]);
[21481]341 if(!node->IsActive()) continue;
[18930]342 if(serial_spc[node->Sid()]==1.){
343 pressure_input->GetInputValue(&pressure, gauss);
344 node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
[16539]345 }
[21481]346 else {
[18930]347 node->DofInFSet(0);
[21481]348 }
[16539]349 }
350
[18930]351 /*Free ressources:*/
352 xDelete<int>(indices);
353 delete gauss;
354}/*}}}*/
355void EnthalpyAnalysis::ComputeBasalMeltingrate(FemModel* femmodel){/*{{{*/
356 /*Compute basal melting rates: */
[25539]357 for(Object* & object : femmodel->elements->objects){
358 Element* element=xDynamicCast<Element*>(object);
[18930]359 ComputeBasalMeltingrate(element);
360 }
[21895]361
362 /*extrude inputs*/
[21894]363 femmodel->parameters->SetParam(BasalforcingsGroundediceMeltingRateEnum,InputToExtrudeEnum);
364 extrudefrombase_core(femmodel);
[18930]365}/*}}}*/
366void EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/
367 /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
368 /* melting rate is positive when melting, negative when refreezing*/
[16539]369
[18930]370 /* Check if ice in element */
371 if(!element->IsIceInElement()) return;
[16539]372
[18930]373 /* Only compute melt rates at the base of grounded ice*/
374 if(!element->IsOnBase() || element->IsFloating()) return;
[16539]375
[18930]376 /* Intermediaries */
377 bool converged;
378 const int dim=3;
379 int i,is,state;
[21718]380 int nodedown,nodeup,numnodes,numsegments;
[18930]381 int enthalpy_enum;
382 IssmDouble vec_heatflux[dim],normal_base[dim],d1enthalpy[dim],d1pressure[dim];
383 IssmDouble basalfriction,alpha2,geothermalflux,heatflux;
384 IssmDouble dt,yts;
385 IssmDouble melting_overshoot,lambda;
386 IssmDouble vx,vy,vz;
387 IssmDouble *xyz_list = NULL;
388 IssmDouble *xyz_list_base = NULL;
389 int *pairindices = NULL;
[16539]390
[18930]391 /*Fetch parameters*/
392 element->GetVerticesCoordinates(&xyz_list);
393 element->GetVerticesCoordinatesBase(&xyz_list_base);
394 element->GetInputValue(&converged,ConvergedEnum);
395 element->FindParam(&dt,TimesteppingTimeStepEnum);
396 element->FindParam(&yts, ConstantsYtsEnum);
397
398 if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum;
399 else enthalpy_enum=EnthalpyEnum;
400
[23644]401 IssmDouble latentheat = element->FindParam(MaterialsLatentheatEnum);
402 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
403 IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
404 IssmDouble beta = element->FindParam(MaterialsBetaEnum);
[18930]405 IssmDouble kappa = EnthalpyDiffusionParameterVolume(element,enthalpy_enum); _assert_(kappa>=0.);
406 IssmDouble kappa_mix;
407
408 /*retrieve inputs*/
[25379]409 Input* enthalpy_input = element->GetInput(enthalpy_enum); _assert_(enthalpy_input);
410 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
411 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
412 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
413 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
414 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
[18930]415
416 /*Build friction element, needed later: */
417 Friction* friction=new Friction(element,dim);
418
419 /******** MELTING RATES ************************************//*{{{*/
420 element->NormalBase(&normal_base[0],xyz_list_base);
[21718]421 element->VerticalSegmentIndicesBase(&pairindices,&numsegments);
[18930]422 IssmDouble* meltingrate_enthalpy = xNew<IssmDouble>(numsegments);
[24240]423 IssmDouble* heating = xNew<IssmDouble>(numsegments);
[18930]424
[21718]425 numnodes=element->GetNumberOfNodes();
426 IssmDouble* enthalpies = xNew<IssmDouble>(numnodes);
427 IssmDouble* pressures = xNew<IssmDouble>(numnodes);
428 IssmDouble* watercolumns = xNew<IssmDouble>(numnodes);
429 IssmDouble* basalmeltingrates = xNew<IssmDouble>(numnodes);
430 element->GetInputListOnNodes(enthalpies,enthalpy_enum);
431 element->GetInputListOnNodes(pressures,PressureEnum);
432 element->GetInputListOnNodes(watercolumns,WatercolumnEnum);
433 element->GetInputListOnNodes(basalmeltingrates,BasalforcingsGroundediceMeltingRateEnum);
[18930]434
[24140]435 IssmDouble watercolumnupperlimit = element->FindParam(ThermalWatercolumnUpperlimitEnum);
[25539]436
[18930]437 Gauss* gauss=element->NewGauss();
438 for(is=0;is<numsegments;is++){
[21718]439 nodedown = pairindices[is*2+0];
440 nodeup = pairindices[is*2+1];
441 gauss->GaussNode(element->GetElementType(),nodedown);
[18930]442
[21718]443 state=GetThermalBasalCondition(element, enthalpies[nodedown], enthalpies[nodeup], pressures[nodedown], pressures[nodeup], watercolumns[nodedown], basalmeltingrates[nodedown]);
[18930]444 switch (state) {
445 case 0:
446 // cold, dry base: apply basal surface forcing
447 for(i=0;i<3;i++) vec_heatflux[i]=0.;
448 break;
[24240]449 case 1: case 2: case 3:
450 // case 1 : cold, wet base: keep at pressure melting point
[18930]451 // case 2: temperate, thin refreezing base: release spc
452 // case 3: temperate, thin melting base: set spc
453 enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],xyz_list,gauss);
454 for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
455 break;
456 case 4:
457 // temperate, thick melting base: set grad H*n=0
[21718]458 kappa_mix=GetWetIceConductivity(element, enthalpies[nodedown], pressures[nodedown]);
[18930]459 pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
460 for(i=0;i<3;i++) vec_heatflux[i]=kappa_mix*beta*d1pressure[i];
461 break;
462 default:
463 _printf0_(" unknown thermal basal state found!");
464 }
465 if(state==0) meltingrate_enthalpy[is]=0.;
466 else{
467 /*heat flux along normal*/
468 heatflux=0.;
469 for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
470
471 /*basal friction*/
472 friction->GetAlpha2(&alpha2,gauss);
473 vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); vz_input->GetInputValue(&vz,gauss);
474 basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
475 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
476 /* -Mb= Fb-(q-q_geo)/((1-w)*L*rho), and (1-w)*rho=rho_ice, cf Aschwanden 2012, eqs.1, 2, 66*/
477 heating[is]=(heatflux+basalfriction+geothermalflux);
478 meltingrate_enthalpy[is]=heating[is]/(latentheat*rho_ice); // m/s water equivalent
479 }
480 }/*}}}*/
481
482 /******** UPDATE MELTINGRATES AND WATERCOLUMN **************//*{{{*/
483 for(is=0;is<numsegments;is++){
[21718]484 nodedown = pairindices[is*2+0];
485 nodeup = pairindices[is*2+1];
[18930]486 if(dt!=0.){
[24240]487 if(watercolumns[nodedown]+meltingrate_enthalpy[is]*dt<0.){ // prevent too much freeze on
[21718]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]; // restrict freeze on only to size of watercolumn
491 enthalpies[nodedown]+=(1.-lambda)*dt/yts*meltingrate_enthalpy[is]*latentheat*rho_ice; // use rest of energy to cool down base: dE=L*m, m=(1-lambda)*meltingrate*rho_ice
[16539]492 }
[18930]493 else{
[21718]494 basalmeltingrates[nodedown]=meltingrate_enthalpy[is];
[24240]495 watercolumns[nodedown]+=dt*meltingrate_enthalpy[is];
[18930]496 }
[24140]497 if(watercolumns[nodedown]>watercolumnupperlimit) watercolumns[nodedown]=watercolumnupperlimit;
[16539]498 }
[18930]499 else{
[21718]500 basalmeltingrates[nodedown]=meltingrate_enthalpy[is];
501 if(watercolumns[nodedown]+meltingrate_enthalpy[is]<0.)
502 watercolumns[nodedown]=0.;
[18930]503 else
[21718]504 watercolumns[nodedown]+=meltingrate_enthalpy[is];
[24240]505 }
[21718]506 basalmeltingrates[nodedown]*=rho_water/rho_ice; // convert meltingrate from water to ice equivalent
507 _assert_(watercolumns[nodedown]>=0.);
[18930]508 }/*}}}*/
509
510 /*feed updated variables back into model*/
[24486]511 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
[18930]512 if(dt!=0.){
[25379]513 element->AddInput(enthalpy_enum,enthalpies,finite_element);
514 element->AddInput(WatercolumnEnum,watercolumns,finite_element);
[16539]515 }
[25379]516 element->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,P1DGEnum);
[16539]517
[18930]518 /*Clean up and return*/
519 delete gauss;
520 delete friction;
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);
[16539]530}/*}}}*/
[18930]531void EnthalpyAnalysis::Core(FemModel* femmodel){/*{{{*/
[21160]532
533 IssmDouble dt;
534 bool isdynamicbasalspc;
535
536 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
537 femmodel->parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
538
[18591]539 if(VerboseSolution()) _printf0_(" computing enthalpy\n");
540 femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
[21160]541 if((dt>0.) && isdynamicbasalspc) UpdateBasalConstraints(femmodel);
[18591]542 solutionsequence_thermal_nonlinear(femmodel);
543
544 /*transfer enthalpy to enthalpy picard for the next step: */
545 InputDuplicatex(femmodel,EnthalpyEnum,EnthalpyPicardEnum);
546
[21161]547 PostProcessing(femmodel);
[18591]548
[17005]549}/*}}}*/
[26047]550void EnthalpyAnalysis::PreCore(FemModel* femmodel){/*{{{*/
551 _error_("not implemented");
552}/*}}}*/
[17000]553ElementVector* EnthalpyAnalysis::CreateDVector(Element* element){/*{{{*/
554 /*Default, return NULL*/
555 return NULL;
556}/*}}}*/
[16992]557ElementMatrix* EnthalpyAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
558_error_("Not implemented");
559}/*}}}*/
[16782]560ElementMatrix* EnthalpyAnalysis::CreateKMatrix(Element* element){/*{{{*/
[16888]561
[17434]562 /* Check if ice in element */
563 if(!element->IsIceInElement()) return NULL;
564
[16888]565 /*compute all stiffness matrices for this element*/
566 ElementMatrix* Ke1=CreateKMatrixVolume(element);
567 ElementMatrix* Ke2=CreateKMatrixShelf(element);
568 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
569
570 /*clean-up and return*/
571 delete Ke1;
572 delete Ke2;
573 return Ke;
[16782]574}/*}}}*/
[16888]575ElementMatrix* EnthalpyAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
576
[17434]577 /* Check if ice in element */
578 if(!element->IsIceInElement()) return NULL;
579
[16888]580 /*Intermediaries */
581 int stabilization;
582 IssmDouble Jdet,dt,u,v,w,um,vm,wm,vel;
583 IssmDouble h,hx,hy,hz,vx,vy,vz;
584 IssmDouble tau_parameter,diameter;
[25539]585 IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
[16888]586 IssmDouble D_scalar;
587 IssmDouble* xyz_list = NULL;
588
589 /*Fetch number of nodes and dof for this finite element*/
590 int numnodes = element->GetNumberOfNodes();
591
592 /*Initialize Element vector and other vectors*/
593 ElementMatrix* Ke = element->NewElementMatrix();
594 IssmDouble* basis = xNew<IssmDouble>(numnodes);
595 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
596 IssmDouble K[3][3];
597
598 /*Retrieve all inputs and parameters*/
599 element->GetVerticesCoordinates(&xyz_list);
600 element->FindParam(&dt,TimesteppingTimeStepEnum);
601 element->FindParam(&stabilization,ThermalStabilizationEnum);
[23644]602 IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
603 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
604 IssmDouble gravity = element->FindParam(ConstantsGEnum);
605 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
606 IssmDouble thermalconductivity = element->FindParam(MaterialsThermalconductivityEnum);
[25379]607 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
608 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
609 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
610 Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input);
611 Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
612 Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input);
[16888]613
614 /*Enthalpy diffusion parameter*/
[17027]615 IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
[16888]616
617 /* Start looping on the number of gaussian points: */
[19637]618 Gauss* gauss=element->NewGauss(4);
[25437]619 while(gauss->next()){
[16888]620
621 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
[22511]622 element->NodalFunctions(basis,gauss);
623 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
624
[16888]625 D_scalar=gauss->weight*Jdet;
626 if(dt!=0.) D_scalar=D_scalar*dt;
627
628 /*Conduction: */
[22511]629 for(int i=0;i<numnodes;i++){
630 for(int j=0;j<numnodes;j++){
631 Ke->values[i*numnodes+j] += D_scalar*kappa/rho_ice*(
632 dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
633 );
634 }
635 }
[16888]636
637 /*Advection: */
638 vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
639 vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
640 vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
[22511]641 for(int i=0;i<numnodes;i++){
642 for(int j=0;j<numnodes;j++){
643 Ke->values[i*numnodes+j] += D_scalar*(
644 vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i] +vz*dbasis[2*numnodes+j]*basis[i]
645 );
646 }
647 }
[16888]648
649 /*Transient: */
650 if(dt!=0.){
651 D_scalar=gauss->weight*Jdet;
[22511]652 for(int i=0;i<numnodes;i++){
653 for(int j=0;j<numnodes;j++){
654 Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
655 }
656 }
[16888]657 D_scalar=D_scalar*dt;
658 }
659
[21382]660 /*Artificial diffusivity*/
[16888]661 if(stabilization==1){
662 element->ElementSizes(&hx,&hy,&hz);
663 vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
664 h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
[18484]665 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);
666 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);
667 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);
[16888]668 for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
669
[24918]670 for(int i=0;i<numnodes;i++){
671 for(int j=0;j<numnodes;j++){
672 Ke->values[i*numnodes+j] += (
673 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]) +
674 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]) +
[25539]675 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])
[24918]676 );
677 }
678 }
[16888]679 }
[24136]680 /*SUPG*/
[16888]681 else if(stabilization==2){
682 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
[25539]683 diameter=element->MinEdgeLength(xyz_list);
[16888]684 tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice);
685 for(int i=0;i<numnodes;i++){
686 for(int j=0;j<numnodes;j++){
687 Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
[24136]688 ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*
689 ((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
[16888]690 }
691 }
692 if(dt!=0.){
[16896]693 D_scalar=gauss->weight*Jdet;
[16888]694 for(int i=0;i<numnodes;i++){
695 for(int j=0;j<numnodes;j++){
[16895]696 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]);
[16888]697 }
698 }
699 }
700 }
[24136]701 /*anisotropic SUPG*/
702 else if(stabilization==3){
703 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
[25539]704 element->ElementSizes(&hx,&hy,&hz);
[24136]705 element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u-um,v-vm,w-wm,hx,hy,hz,kappa/rho_ice);
706 tau_parameter_hor=tau_parameter_anisotropic[0];
707 tau_parameter_ver=tau_parameter_anisotropic[1];
708 for(int i=0;i<numnodes;i++){
709 for(int j=0;j<numnodes;j++){
710 Ke->values[i*numnodes+j]+=D_scalar*
711 (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])*
712 (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]);
713 }
714 }
715 }
[16888]716 }
717
718 /*Clean up and return*/
719 xDelete<IssmDouble>(xyz_list);
720 xDelete<IssmDouble>(basis);
721 xDelete<IssmDouble>(dbasis);
722 delete gauss;
723 return Ke;
724}/*}}}*/
725ElementMatrix* EnthalpyAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/
726
[17434]727 /* Check if ice in element */
728 if(!element->IsIceInElement()) return NULL;
729
[16888]730 /*Initialize Element matrix and return if necessary*/
[17585]731 if(!element->IsOnBase() || !element->IsFloating()) return NULL;
[16888]732
[16986]733 /*Intermediaries*/
[16888]734 IssmDouble dt,Jdet,D;
735 IssmDouble *xyz_list_base = NULL;
736
737 /*Fetch number of nodes for this finite element*/
738 int numnodes = element->GetNumberOfNodes();
739
740 /*Initialize vectors*/
741 ElementMatrix* Ke = element->NewElementMatrix();
742 IssmDouble* basis = xNew<IssmDouble>(numnodes);
743
744 /*Retrieve all inputs and parameters*/
745 element->GetVerticesCoordinatesBase(&xyz_list_base);
746 element->FindParam(&dt,TimesteppingTimeStepEnum);
[23644]747 IssmDouble gravity = element->FindParam(ConstantsGEnum);
748 IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
749 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
750 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
751 IssmDouble mixed_layer_capacity= element->FindParam(MaterialsMixedLayerCapacityEnum);
752 IssmDouble thermal_exchange_vel= element->FindParam(MaterialsThermalExchangeVelocityEnum);
[16888]753
754 /* Start looping on the number of gaussian points: */
[19637]755 Gauss* gauss=element->NewGaussBase(4);
[25437]756 while(gauss->next()){
[16888]757
758 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
759 element->NodalFunctions(basis,gauss);
760
761 D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
762 if(reCast<bool,IssmDouble>(dt)) D=dt*D;
[24918]763 for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
[16888]764 }
765
766 /*Clean up and return*/
767 delete gauss;
768 xDelete<IssmDouble>(basis);
769 xDelete<IssmDouble>(xyz_list_base);
770 return Ke;
771}/*}}}*/
[16782]772ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/
[16812]773
[17434]774 /* Check if ice in element */
775 if(!element->IsIceInElement()) return NULL;
776
[16812]777 /*compute all load vectors for this element*/
778 ElementVector* pe1=CreatePVectorVolume(element);
779 ElementVector* pe2=CreatePVectorSheet(element);
780 ElementVector* pe3=CreatePVectorShelf(element);
781 ElementVector* pe =new ElementVector(pe1,pe2,pe3);
782
783 /*clean-up and return*/
784 delete pe1;
785 delete pe2;
786 delete pe3;
787 return pe;
[16782]788}/*}}}*/
[16812]789ElementVector* EnthalpyAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
790
[17434]791 /* Check if ice in element */
792 if(!element->IsIceInElement()) return NULL;
793
[16812]794 /*Intermediaries*/
[17014]795 int i, stabilization;
[16812]796 IssmDouble Jdet,phi,dt;
[17014]797 IssmDouble enthalpy, Hpmp;
798 IssmDouble enthalpypicard, d1enthalpypicard[3];
799 IssmDouble pressure, d1pressure[3], d2pressure;
800 IssmDouble waterfractionpicard;
[24136]801 IssmDouble kappa,tau_parameter,diameter,hx,hy,hz,kappa_w;
802 IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
[16812]803 IssmDouble u,v,w;
[17014]804 IssmDouble scalar_def, scalar_sens ,scalar_transient;
[16812]805 IssmDouble* xyz_list = NULL;
[17014]806 IssmDouble d1H_d1P, d1P2;
[16812]807
808 /*Fetch number of nodes and dof for this finite element*/
809 int numnodes = element->GetNumberOfNodes();
810 int numvertices = element->GetNumberOfVertices();
811
812 /*Initialize Element vector*/
813 ElementVector* pe = element->NewElementVector();
814 IssmDouble* basis = xNew<IssmDouble>(numnodes);
815 IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
816
817 /*Retrieve all inputs and parameters*/
818 element->GetVerticesCoordinates(&xyz_list);
[23644]819 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
820 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
821 IssmDouble thermalconductivity = element->FindParam(MaterialsThermalconductivityEnum);
822 IssmDouble temperateiceconductivity = element->FindParam(MaterialsTemperateiceconductivityEnum);
823 IssmDouble beta = element->FindParam(MaterialsBetaEnum);
824 IssmDouble latentheat = element->FindParam(MaterialsLatentheatEnum);
[16812]825 element->FindParam(&dt,TimesteppingTimeStepEnum);
826 element->FindParam(&stabilization,ThermalStabilizationEnum);
[25379]827 Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
828 Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
829 Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
830 Input* enthalpypicard_input=element->GetInput(EnthalpyPicardEnum); _assert_(enthalpypicard_input);
831 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
832 Input* enthalpy_input=NULL;
[24556]833 if(dt>0.){
[25379]834 enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
[24556]835 }
[16812]836
837 /* Start looping on the number of gaussian points: */
[19637]838 Gauss* gauss=element->NewGauss(4);
[25437]839 while(gauss->next()){
[16812]840
841 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
842 element->NodalFunctions(basis,gauss);
[23066]843
[17014]844 /*viscous dissipation*/
[16812]845 element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
846
847 scalar_def=phi/rho_ice*Jdet*gauss->weight;
[16895]848 if(dt!=0.) scalar_def=scalar_def*dt;
[16812]849
[17014]850 for(i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
[16812]851
[17014]852 /*sensible heat flux in temperate ice*/
853 enthalpypicard_input->GetInputValue(&enthalpypicard,gauss);
854 pressure_input->GetInputValue(&pressure,gauss);
855 Hpmp=this->PureIceEnthalpy(element, pressure);
856
857 if(enthalpypicard>=Hpmp){
858 enthalpypicard_input->GetInputDerivativeValue(&d1enthalpypicard[0],xyz_list,gauss);
859 pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
860 d2pressure=0.; // for linear elements, 2nd derivative is zero
[23066]861
[17014]862 d1H_d1P=0.;
863 for(i=0;i<3;i++) d1H_d1P+=d1enthalpypicard[i]*d1pressure[i];
864 d1P2=0.;
865 for(i=0;i<3;i++) d1P2+=pow(d1pressure[i],2.);
866
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];
[24240]870 }
[17014]871
[16812]872 /* Build transient now */
[24556]873 if(dt>0.){
[16812]874 enthalpy_input->GetInputValue(&enthalpy, gauss);
875 scalar_transient=enthalpy*Jdet*gauss->weight;
[17014]876 for(i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
[16812]877 }
[25539]878
[24136]879 /* SUPG */
[16812]880 if(stabilization==2){
881 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
[24136]882 diameter=element->MinEdgeLength(xyz_list);
883 kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
[16812]884 vx_input->GetInputValue(&u,gauss);
885 vy_input->GetInputValue(&v,gauss);
886 vz_input->GetInputValue(&w,gauss);
[16895]887 tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa/rho_ice);
[16812]888
[17014]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]);
[16895]890
891 if(dt!=0.){
[17014]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]);
[16812]893 }
894 }
[24136]895 /* anisotropic SUPG */
896 else if(stabilization==3){
897 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
898 element->ElementSizes(&hx,&hy,&hz);
899 kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.);
[24556]900 vx_input->GetInputValue(&u,gauss);
901 vy_input->GetInputValue(&v,gauss);
902 vz_input->GetInputValue(&w,gauss);
903 element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa/rho_ice);
904 tau_parameter_hor=tau_parameter_anisotropic[0];
905 tau_parameter_ver=tau_parameter_anisotropic[1];
906
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]);
[24136]908 }
[16812]909 }
910
911 /*Clean up and return*/
912 xDelete<IssmDouble>(basis);
913 xDelete<IssmDouble>(dbasis);
914 xDelete<IssmDouble>(xyz_list);
915 delete gauss;
916 return pe;
917
918}/*}}}*/
919ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
[16888]920
[17434]921 /* Check if ice in element */
922 if(!element->IsIceInElement()) return NULL;
923
[17014]924 /* implementation of the basal condition decision chart of Aschwanden 2012, Fig.5 */
[17585]925 if(!element->IsOnBase() || element->IsFloating()) return NULL;
[16888]926
[20272]927 bool converged, isdynamicbasalspc;
[18612]928 int i, state;
[20272]929 int enthalpy_enum;
[18612]930 IssmDouble dt,Jdet,scalar;
931 IssmDouble enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
932 IssmDouble vx,vy,vz;
933 IssmDouble alpha2,basalfriction,geothermalflux,heatflux;
[16888]934 IssmDouble *xyz_list_base = NULL;
935
936 /*Fetch number of nodes for this finite element*/
937 int numnodes = element->GetNumberOfNodes();
938
939 /*Initialize vectors*/
940 ElementVector* pe = element->NewElementVector();
941 IssmDouble* basis = xNew<IssmDouble>(numnodes);
942
943 /*Retrieve all inputs and parameters*/
944 element->GetVerticesCoordinatesBase(&xyz_list_base);
945 element->FindParam(&dt,TimesteppingTimeStepEnum);
[18622]946 element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
[20272]947 element->GetInputValue(&converged,ConvergedEnum);
948 if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum; // use enthalpy from last iteration
949 else enthalpy_enum=EnthalpyEnum; // use enthalpy from last time step
[25379]950 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
951 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
952 Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
953 Input* enthalpy_input = element->GetInput(enthalpy_enum); _assert_(enthalpy_input);
954 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
955 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
956 Input* meltingrate_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(meltingrate_input);
957 Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
[23644]958 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
[16888]959
960 /*Build friction element, needed later: */
961 Friction* friction=new Friction(element,3);
962
963 /* Start looping on the number of gaussian points: */
[19637]964 Gauss* gauss=element->NewGaussBase(4);
965 Gauss* gaussup=element->NewGaussTop(4);
[25437]966 while(gauss->next() && gaussup->next()){
[16888]967
968 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
969 element->NodalFunctions(basis,gauss);
970
[18622]971 if(isdynamicbasalspc){
972 enthalpy_input->GetInputValue(&enthalpy,gauss);
973 enthalpy_input->GetInputValue(&enthalpyup,gaussup);
974 pressure_input->GetInputValue(&pressure,gauss);
975 pressure_input->GetInputValue(&pressureup,gaussup);
976 watercolumn_input->GetInputValue(&watercolumn,gauss);
977 meltingrate_input->GetInputValue(&meltingrate,gauss);
978 state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
979 }
980 else
981 state=0;
[16888]982
[18612]983 switch (state) {
[20272]984 case 0: case 1: case 2: case 3:
[24240]985 // cold, dry base; cold, wet base; refreezing temperate base; thin temperate base:
[20272]986 // Apply basal surface forcing.
[24240]987 // Interpolated values of enthalpy on gauss nodes may indicate cold base,
[20272]988 // although one node might have become temperate. So keep heat flux switched on.
[18612]989 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
990 friction->GetAlpha2(&alpha2,gauss);
991 vx_input->GetInputValue(&vx,gauss);
992 vy_input->GetInputValue(&vy,gauss);
993 vz_input->GetInputValue(&vz,gauss);
994 basalfriction=alpha2*(vx*vx+vy*vy+vz*vz);
995 heatflux=(basalfriction+geothermalflux)/(rho_ice);
996 scalar=gauss->weight*Jdet*heatflux;
997 if(dt!=0.) scalar=dt*scalar;
[24240]998 for(i=0;i<numnodes;i++)
[18612]999 pe->values[i]+=scalar*basis[i];
1000 break;
1001 case 4:
1002 // temperate, thick melting base: set grad H*n=0
[24240]1003 for(i=0;i<numnodes;i++)
[18612]1004 pe->values[i]+=0.;
1005 break;
1006 default:
1007 _printf0_(" unknown thermal basal state found!");
[16888]1008 }
1009 }
1010
1011 /*Clean up and return*/
1012 delete gauss;
1013 delete gaussup;
1014 delete friction;
1015 xDelete<IssmDouble>(basis);
1016 xDelete<IssmDouble>(xyz_list_base);
1017 return pe;
1018
[16812]1019}/*}}}*/
1020ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
1021
[17434]1022 /* Check if ice in element */
1023 if(!element->IsIceInElement()) return NULL;
1024
[16888]1025 /*Get basal element*/
[17585]1026 if(!element->IsOnBase() || !element->IsFloating()) return NULL;
[16888]1027
[18612]1028 IssmDouble Hpmp,dt,Jdet,scalar_ocean,pressure;
[16812]1029 IssmDouble *xyz_list_base = NULL;
1030
1031 /*Fetch number of nodes for this finite element*/
1032 int numnodes = element->GetNumberOfNodes();
1033
1034 /*Initialize vectors*/
1035 ElementVector* pe = element->NewElementVector();
1036 IssmDouble* basis = xNew<IssmDouble>(numnodes);
1037
1038 /*Retrieve all inputs and parameters*/
1039 element->GetVerticesCoordinatesBase(&xyz_list_base);
1040 element->FindParam(&dt,TimesteppingTimeStepEnum);
[25379]1041 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
[23644]1042 IssmDouble gravity = element->FindParam(ConstantsGEnum);
1043 IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
1044 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
1045 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
1046 IssmDouble mixed_layer_capacity= element->FindParam(MaterialsMixedLayerCapacityEnum);
1047 IssmDouble thermal_exchange_vel= element->FindParam(MaterialsThermalExchangeVelocityEnum);
[16812]1048
1049 /* Start looping on the number of gaussian points: */
[19637]1050 Gauss* gauss=element->NewGaussBase(4);
[25437]1051 while(gauss->next()){
[16812]1052
1053 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
1054 element->NodalFunctions(basis,gauss);
1055
1056 pressure_input->GetInputValue(&pressure,gauss);
[18612]1057 Hpmp=element->PureIceEnthalpy(pressure);
[16812]1058
[18612]1059 scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*Hpmp/(heatcapacity*rho_ice);
[16812]1060 if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
1061
1062 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_ocean*basis[i];
1063 }
1064
1065 /*Clean up and return*/
1066 delete gauss;
1067 xDelete<IssmDouble>(basis);
1068 xDelete<IssmDouble>(xyz_list_base);
1069 return pe;
1070}/*}}}*/
[18930]1071void EnthalpyAnalysis::DrainWaterfraction(FemModel* femmodel){/*{{{*/
1072 /*Drain excess water fraction in ice column: */
[21721]1073 ComputeWaterfractionDrainage(femmodel);
1074 DrainageUpdateWatercolumn(femmodel);
1075 DrainageUpdateEnthalpy(femmodel);
[17002]1076}/*}}}*/
[21721]1077void EnthalpyAnalysis::ComputeWaterfractionDrainage(FemModel* femmodel){/*{{{*/
[17002]1078
[25539]1079 int k,numnodes;
[21721]1080 IssmDouble dt;
1081 Element* element= NULL;
[17434]1082
[21721]1083 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
[17434]1084
[25539]1085 for(Object* & object : femmodel->elements->objects){
1086 element=xDynamicCast<Element*>(object);
[21721]1087 numnodes=element->GetNumberOfNodes();
1088 IssmDouble* waterfractions= xNew<IssmDouble>(numnodes);
1089 IssmDouble* drainage= xNew<IssmDouble>(numnodes);
[17014]1090
[21721]1091 element->GetInputListOnNodes(waterfractions,WaterfractionEnum);
1092 for(k=0; k<numnodes;k++){
1093 drainage[k]=DrainageFunctionWaterfraction(waterfractions[k], dt);
1094 }
[24486]1095 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
[25379]1096 element->AddInput(WaterfractionDrainageEnum,drainage,finite_element);
[18930]1097
[21721]1098 xDelete<IssmDouble>(waterfractions);
1099 xDelete<IssmDouble>(drainage);
[18930]1100 }
[21721]1101}/*}}}*/
1102void EnthalpyAnalysis::DrainageUpdateWatercolumn(FemModel* femmodel){/*{{{*/
1103
[25539]1104 int k,numnodes, numbasalnodes;
[21721]1105 IssmDouble dt;
1106 int* basalnodeindices=NULL;
1107 Element* element= NULL;
[23066]1108
[21721]1109 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
[18667]1110
[21721]1111 /*depth-integrate the drained water fraction */
1112 femmodel->parameters->SetParam(WaterfractionDrainageEnum,InputToDepthaverageInEnum);
1113 femmodel->parameters->SetParam(WaterfractionDrainageIntegratedEnum,InputToDepthaverageOutEnum);
1114 depthaverage_core(femmodel);
1115 femmodel->parameters->SetParam(WaterfractionDrainageIntegratedEnum,InputToExtrudeEnum);
1116 extrudefrombase_core(femmodel);
1117 /*multiply depth-average by ice thickness*/
[25539]1118 for(Object* & object : femmodel->elements->objects){
1119 element=xDynamicCast<Element*>(object);
[21721]1120 numnodes=element->GetNumberOfNodes();
1121 IssmDouble* drainage_int= xNew<IssmDouble>(numnodes);
1122 IssmDouble* thicknesses= xNew<IssmDouble>(numnodes);
[17014]1123
[21721]1124 element->GetInputListOnNodes(drainage_int,WaterfractionDrainageIntegratedEnum);
1125 element->GetInputListOnNodes(thicknesses,ThicknessEnum);
1126 for(k=0;k<numnodes;k++){
1127 drainage_int[k]*=thicknesses[k];
1128 }
[24486]1129 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
[25379]1130 element->AddInput(WaterfractionDrainageIntegratedEnum, drainage_int,finite_element);
[17166]1131
[21721]1132 xDelete<IssmDouble>(drainage_int);
1133 xDelete<IssmDouble>(thicknesses);
1134 }
[17434]1135
[21721]1136 /*update water column*/
[25539]1137 for(Object* & object : femmodel->elements->objects){
1138 element=xDynamicCast<Element*>(object);
[21721]1139 /* Check if ice in element */
1140 if(!element->IsIceInElement()) continue;
[24240]1141 if(!element->IsOnBase()) continue;
[17166]1142
[21721]1143 numnodes=element->GetNumberOfNodes();
1144 IssmDouble* watercolumn= xNew<IssmDouble>(numnodes);
1145 IssmDouble* drainage_int= xNew<IssmDouble>(numnodes);
1146 element->GetInputListOnNodes(watercolumn,WatercolumnEnum);
1147 element->GetInputListOnNodes(drainage_int,WaterfractionDrainageIntegratedEnum);
[17166]1148
[21721]1149 element->BasalNodeIndices(&numbasalnodes,&basalnodeindices,element->GetElementType());
1150 for(k=0;k<numbasalnodes;k++){
1151 watercolumn[basalnodeindices[k]]+=dt*drainage_int[basalnodeindices[k]];
1152 }
[24486]1153 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
[25379]1154 element->AddInput(WatercolumnEnum, watercolumn,finite_element);
[17166]1155
[21721]1156 xDelete<IssmDouble>(watercolumn);
1157 xDelete<IssmDouble>(drainage_int);
[21779]1158 xDelete<int>(basalnodeindices);
[21721]1159 }
1160}/*}}}*/
1161void EnthalpyAnalysis::DrainageUpdateEnthalpy(FemModel* femmodel){/*{{{*/
[17166]1162
[25539]1163 int k,numnodes;
[21721]1164 IssmDouble dt;
1165 Element* element= NULL;
1166 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
[17166]1167
[25539]1168 for(Object* & object : femmodel->elements->objects){
1169 element=xDynamicCast<Element*>(object);
[21721]1170 numnodes=element->GetNumberOfNodes();
1171 IssmDouble* enthalpies= xNew<IssmDouble>(numnodes);
1172 IssmDouble* pressures= xNew<IssmDouble>(numnodes);
1173 IssmDouble* temperatures= xNew<IssmDouble>(numnodes);
1174 IssmDouble* waterfractions= xNew<IssmDouble>(numnodes);
1175 IssmDouble* drainage= xNew<IssmDouble>(numnodes);
[21653]1176
[21721]1177 element->GetInputListOnNodes(pressures,PressureEnum);
1178 element->GetInputListOnNodes(temperatures,TemperatureEnum);
1179 element->GetInputListOnNodes(waterfractions,WaterfractionEnum);
1180 element->GetInputListOnNodes(drainage,WaterfractionDrainageEnum);
[21653]1181
[21721]1182 for(k=0;k<numnodes;k++){
[24140]1183 if(dt==0.)
1184 waterfractions[k]-=drainage[k];
1185 else
1186 waterfractions[k]-=dt*drainage[k];
[25539]1187
[21721]1188 element->ThermalToEnthalpy(&enthalpies[k], temperatures[k], waterfractions[k], pressures[k]);
1189 }
[24486]1190 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
[25379]1191 element->AddInput(WaterfractionEnum,waterfractions,finite_element);
1192 element->AddInput(EnthalpyEnum,enthalpies,finite_element);
[17166]1193
[21721]1194 xDelete<IssmDouble>(enthalpies);
1195 xDelete<IssmDouble>(pressures);
1196 xDelete<IssmDouble>(temperatures);
1197 xDelete<IssmDouble>(waterfractions);
1198 xDelete<IssmDouble>(drainage);
1199 }
[17002]1200}/*}}}*/
[18930]1201IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
[17014]1202
[23644]1203 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
1204 IssmDouble temperateiceconductivity = element->FindParam(MaterialsTemperateiceconductivityEnum);
1205 IssmDouble thermalconductivity = element->FindParam(MaterialsThermalconductivityEnum);
[17434]1206
[24140]1207 if(enthalpy < PureIceEnthalpy(element,pressure))
[18930]1208 return thermalconductivity/heatcapacity;
[24140]1209 else
[18930]1210 return temperateiceconductivity/heatcapacity;
1211}/*}}}*/
1212IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum){/*{{{*/
[17014]1213
[18930]1214 int iv;
1215 IssmDouble lambda; /* fraction of cold ice */
1216 IssmDouble kappa,kappa_c,kappa_t; /* enthalpy conductivities */
1217 IssmDouble Hc,Ht;
[17014]1218
[18930]1219 /*Get pressures and enthalpies on vertices*/
1220 int numvertices = element->GetNumberOfVertices();
[23697]1221 int effectiveconductivity_averaging;
[18930]1222 IssmDouble* pressures = xNew<IssmDouble>(numvertices);
1223 IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
1224 IssmDouble* PIE = xNew<IssmDouble>(numvertices);
1225 IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
[17014]1226 element->GetInputListOnVertices(pressures,PressureEnum);
[18930]1227 element->GetInputListOnVertices(enthalpies,enthalpy_enum);
[23697]1228 element->FindParam(&effectiveconductivity_averaging,MaterialsEffectiveconductivityAveragingEnum);
1229
[18930]1230 for(iv=0;iv<numvertices;iv++){
1231 PIE[iv] = PureIceEnthalpy(element,pressures[iv]);
1232 dHpmp[iv] = enthalpies[iv]-PIE[iv];
1233 }
[17014]1234
[18930]1235 bool allequalsign = true;
1236 if(dHpmp[0]<0.){
1237 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0.));
[17014]1238 }
[18930]1239 else{
1240 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0.));
[17014]1241 }
1242
[18930]1243 if(allequalsign){
1244 kappa = EnthalpyDiffusionParameter(element,enthalpies[0],pressures[0]);
[17014]1245 }
[18930]1246 else{
1247 kappa_c = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)-1.,0.);
1248 kappa_t = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)+1.,0.);
[23697]1249
[18930]1250 Hc=0.; Ht=0.;
1251 for(iv=0; iv<numvertices;iv++){
1252 if(enthalpies[iv]<PIE[iv])
1253 Hc+=(PIE[iv]-enthalpies[iv]);
1254 else
1255 Ht+=(enthalpies[iv]-PIE[iv]);
1256 }
1257 _assert_((Hc+Ht)>0.);
1258 lambda = Hc/(Hc+Ht);
[23317]1259 _assert_(lambda>=0.);
1260 _assert_(lambda<=1.);
[23697]1261
1262 if(effectiveconductivity_averaging==0){
1263 /* return arithmetic mean (volume average) of thermal conductivities, weighted by fraction of cold/temperate ice */
1264 kappa=kappa_c*lambda+(1.-lambda)*kappa_t;
1265 }
1266 else if(effectiveconductivity_averaging==1){
1267 /* return harmonic mean (reciprocal avarage) of thermal conductivities, weighted by fraction of cold/temperate ice, cf Patankar 1980, pp44 */
[24240]1268 kappa=kappa_c*kappa_t/(lambda*kappa_t+(1.-lambda)*kappa_c);
[23697]1269 }
1270 else if(effectiveconductivity_averaging==2){
1271 /* return geometric mean (power law) of thermal conductivities, weighted by fraction of cold/temperate ice */
[24240]1272 kappa=pow(kappa_c,lambda)*pow(kappa_t,1.-lambda);
[23697]1273 }
1274 else{
1275 _error_("effectiveconductivity_averaging not supported yet");
1276 }
[24240]1277 }
[17014]1278
1279 /*Clean up and return*/
[18930]1280 xDelete<IssmDouble>(PIE);
1281 xDelete<IssmDouble>(dHpmp);
1282 xDelete<IssmDouble>(pressures);
[17014]1283 xDelete<IssmDouble>(enthalpies);
[18930]1284 return kappa;
[17002]1285}/*}}}*/
[18930]1286void EnthalpyAnalysis::GetBasalConstraints(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
[18659]1287
1288 /*Intermediary*/
1289 bool isdynamicbasalspc;
[18612]1290 IssmDouble dt;
1291
1292 /*Check wether dynamic basal boundary conditions are activated */
1293 element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
1294 if(!isdynamicbasalspc) return;
1295
1296 element->FindParam(&dt,TimesteppingTimeStepEnum);
1297 if(dt==0.){
[18659]1298 GetBasalConstraintsSteadystate(vec_spc,element);
[18612]1299 }
1300 else{
[18659]1301 GetBasalConstraintsTransient(vec_spc,element);
[18612]1302 }
1303}/*}}}*/
[18930]1304void EnthalpyAnalysis::GetBasalConstraintsSteadystate(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
[18612]1305
1306 /* Check if ice in element */
1307 if(!element->IsIceInElement()) return;
1308
[24240]1309 /* Only update constraints at the base.
[20213]1310 * Floating ice is not affected by basal BC decision chart. */
[18612]1311 if(!(element->IsOnBase()) || element->IsFloating()) return;
1312
1313 /*Intermediary*/
1314 int numindices, numindicesup, state;
[17027]1315 int *indices = NULL, *indicesup = NULL;
[18612]1316 IssmDouble enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
[17014]1317
[18612]1318 /*Get parameters and inputs: */
[25379]1319 Input* enthalpy_input = element->GetInput(EnthalpyPicardEnum); _assert_(enthalpy_input);
1320 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
1321 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
1322 Input* meltingrate_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(meltingrate_input);
[18612]1323
[17027]1324 /*Fetch indices of basal & surface nodes for this finite element*/
1325 Penta *penta = (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
1326 penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
[18612]1327 penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType()); _assert_(numindices==numindicesup);
[17014]1328
[18612]1329 GaussPenta* gauss=new GaussPenta();
1330 GaussPenta* gaussup=new GaussPenta();
1331 for(int i=0;i<numindices;i++){
1332 gauss->GaussNode(element->GetElementType(),indices[i]);
1333 gaussup->GaussNode(element->GetElementType(),indicesup[i]);
[18930]1334
[18612]1335 enthalpy_input->GetInputValue(&enthalpy,gauss);
1336 enthalpy_input->GetInputValue(&enthalpyup,gaussup);
1337 pressure_input->GetInputValue(&pressure,gauss);
1338 pressure_input->GetInputValue(&pressureup,gaussup);
1339 watercolumn_input->GetInputValue(&watercolumn,gauss);
1340 meltingrate_input->GetInputValue(&meltingrate,gauss);
1341
1342 state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
1343 switch (state) {
1344 case 0:
1345 // cold, dry base: apply basal surface forcing
[18659]1346 vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
[18612]1347 break;
1348 case 1:
[24240]1349 // cold, wet base: keep at pressure melting point
[18659]1350 vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
[18612]1351 break;
1352 case 2:
[24240]1353 // temperate, thin refreezing base:
[20098]1354 vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
[18612]1355 break;
1356 case 3:
1357 // temperate, thin melting base: set spc
[18659]1358 vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
[18612]1359 break;
1360 case 4:
[20098]1361 // temperate, thick melting base:
[18930]1362 vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
[18612]1363 break;
1364 default:
1365 _printf0_(" unknown thermal basal state found!");
1366 }
1367 }
1368
1369 /*Free ressources:*/
1370 xDelete<int>(indices);
1371 xDelete<int>(indicesup);
1372 delete gauss;
1373 delete gaussup;
1374}/*}}}*/
[18930]1375void EnthalpyAnalysis::GetBasalConstraintsTransient(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
[18612]1376
1377 /* Check if ice in element */
1378 if(!element->IsIceInElement()) return;
1379
[24240]1380 /* Only update constraints at the base.
[20213]1381 * Floating ice is not affected by basal BC decision chart.*/
[18612]1382 if(!(element->IsOnBase()) || element->IsFloating()) return;
1383
1384 /*Intermediary*/
1385 int numindices, numindicesup, state;
1386 int *indices = NULL, *indicesup = NULL;
1387 IssmDouble enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
1388
[17027]1389 /*Get parameters and inputs: */
[25379]1390 Input* enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input); //TODO: check EnthalpyPicard?
1391 Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input);
1392 Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
1393 Input* meltingrate_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(meltingrate_input);
[17014]1394
[18612]1395 /*Fetch indices of basal & surface nodes for this finite element*/
1396 Penta *penta = (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
1397 penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
1398 penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType()); _assert_(numindices==numindicesup);
1399
[17027]1400 GaussPenta* gauss=new GaussPenta();
1401 GaussPenta* gaussup=new GaussPenta();
[18930]1402
[17027]1403 for(int i=0;i<numindices;i++){
1404 gauss->GaussNode(element->GetElementType(),indices[i]);
1405 gaussup->GaussNode(element->GetElementType(),indicesup[i]);
[23066]1406
[18612]1407 enthalpy_input->GetInputValue(&enthalpy,gauss);
1408 enthalpy_input->GetInputValue(&enthalpyup,gaussup);
1409 pressure_input->GetInputValue(&pressure,gauss);
1410 pressure_input->GetInputValue(&pressureup,gaussup);
[17027]1411 watercolumn_input->GetInputValue(&watercolumn,gauss);
[18612]1412 meltingrate_input->GetInputValue(&meltingrate,gauss);
1413
1414 state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
[18930]1415
[18612]1416 switch (state) {
1417 case 0:
1418 // cold, dry base: apply basal surface forcing
[18659]1419 vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
[18612]1420 break;
1421 case 1:
[24240]1422 // cold, wet base: keep at pressure melting point
[18659]1423 vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
[18612]1424 break;
1425 case 2:
1426 // temperate, thin refreezing base: release spc
[18659]1427 vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
[18612]1428 break;
1429 case 3:
1430 // temperate, thin melting base: set spc
[18659]1431 vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
[18612]1432 break;
1433 case 4:
[18930]1434 // temperate, thick melting base: set grad H*n=0
1435 vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
[18612]1436 break;
1437 default:
1438 _printf0_(" unknown thermal basal state found!");
[17027]1439 }
[18930]1440
[17027]1441 }
[17014]1442
[17027]1443 /*Free ressources:*/
1444 xDelete<int>(indices);
1445 xDelete<int>(indicesup);
1446 delete gauss;
1447 delete gaussup;
[17002]1448}/*}}}*/
[18930]1449void EnthalpyAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
[24933]1450 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
[18930]1451 * For node i, Bi' can be expressed in the actual coordinate system
[24240]1452 * by:
[18930]1453 * Bi_conduct=[ dh/dx ]
1454 * [ dh/dy ]
1455 * [ dh/dz ]
1456 * where h is the interpolation function for node i.
1457 *
[24933]1458 * We assume B has been allocated already, of size: 3x(1*numnodes)
[18930]1459 */
[17002]1460
[18930]1461 /*Fetch number of nodes for this finite element*/
1462 int numnodes = element->GetNumberOfNodes();
1463
1464 /*Get nodal functions derivatives*/
1465 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
1466 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
1467
1468 /*Build B: */
1469 for(int i=0;i<numnodes;i++){
1470 B[numnodes*0+i] = dbasis[0*numnodes+i];
1471 B[numnodes*1+i] = dbasis[1*numnodes+i];
1472 B[numnodes*2+i] = dbasis[2*numnodes+i];
1473 }
1474
1475 /*Clean-up*/
1476 xDelete<IssmDouble>(dbasis);
1477}/*}}}*/
1478void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
1479 element->GetSolutionFromInputsOneDof(solution,EnthalpyEnum);
1480}/*}}}*/
1481int EnthalpyAnalysis::GetThermalBasalCondition(Element* element, IssmDouble enthalpy, IssmDouble enthalpyup, IssmDouble pressure, IssmDouble pressureup, IssmDouble watercolumn, IssmDouble meltingrate){/*{{{*/
1482
[18612]1483 /* Check if ice in element */
1484 if(!element->IsIceInElement()) return -1;
1485
1486 /* Only update Constraints at the base of grounded ice*/
1487 if(!(element->IsOnBase())) return -1;
1488
1489 /*Intermediary*/
1490 int state=-1;
1491 IssmDouble dt;
1492
1493 /*Get parameters and inputs: */
1494 element->FindParam(&dt,TimesteppingTimeStepEnum);
1495
[18620]1496 if(enthalpy<PureIceEnthalpy(element,pressure)){
1497 if(watercolumn<=0.) state=0; // cold, dry base
1498 else state=1; // cold, wet base (refreezing)
[18612]1499 }
[18620]1500 else{
1501 if(enthalpyup<PureIceEnthalpy(element,pressureup)){
1502 if((dt==0.) && (meltingrate<0.)) state=2; // refreezing temperate base (non-physical, only for steadystate solver)
1503 else state=3; // temperate base, but no temperate layer
[18612]1504 }
[18620]1505 else state=4; // temperate layer with positive thickness
[18612]1506 }
1507
1508 _assert_(state>=0);
1509 return state;
1510}/*}}}*/
[18930]1511IssmDouble EnthalpyAnalysis::GetWetIceConductivity(Element* element, IssmDouble enthalpy, IssmDouble pressure){/*{{{*/
[18612]1512
1513 IssmDouble temperature, waterfraction;
1514 IssmDouble kappa_w = 0.6; // thermal conductivity of water (in W/m/K)
[23644]1515 IssmDouble kappa_i = element->FindParam(MaterialsThermalconductivityEnum);
[18612]1516 element->EnthalpyToThermal(&temperature, &waterfraction, enthalpy, pressure);
1517
1518 return (1.-waterfraction)*kappa_i + waterfraction*kappa_w;
1519}/*}}}*/
[25317]1520void EnthalpyAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
[18930]1521 _error_("Not implemented yet");
1522}/*}}}*/
1523void EnthalpyAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
[18612]1524
[18930]1525 bool converged;
1526 int i,rheology_law;
1527 IssmDouble B_average,s_average,T_average=0.,P_average=0.;
1528 int *doflist = NULL;
1529 IssmDouble *xyz_list = NULL;
[16888]1530
[18930]1531 /*Fetch number of nodes and dof for this finite element*/
1532 int numnodes = element->GetNumberOfNodes();
[16888]1533
[18930]1534 /*Fetch dof list and allocate solution vector*/
[23629]1535 element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
[18930]1536 IssmDouble* values = xNew<IssmDouble>(numnodes);
1537 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
1538 IssmDouble* surface = xNew<IssmDouble>(numnodes);
1539 IssmDouble* B = xNew<IssmDouble>(numnodes);
1540 IssmDouble* temperature = xNew<IssmDouble>(numnodes);
1541 IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
[16888]1542
[18930]1543 /*Use the dof list to index into the solution vector: */
1544 for(i=0;i<numnodes;i++){
1545 values[i]=solution[doflist[i]];
[16888]1546
[18930]1547 /*Check solution*/
1548 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
[20669]1549 if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
[16888]1550 }
1551
[18930]1552 /*Get all inputs and parameters*/
1553 element->GetInputValue(&converged,ConvergedEnum);
1554 element->GetInputListOnNodes(&pressure[0],PressureEnum);
[24486]1555 int finite_element = element->GetElementType(); if(finite_element==P1Enum) finite_element = P1DGEnum;
[18930]1556 if(converged){
1557 for(i=0;i<numnodes;i++){
1558 element->EnthalpyToThermal(&temperature[i],&waterfraction[i],values[i],pressure[i]);
1559 if(waterfraction[i]<0.) _error_("Negative water fraction found in solution vector");
1560 //if(waterfraction[i]>1.) _error_("Water fraction >1 found in solution vector");
1561 }
[25379]1562 element->AddInput(EnthalpyEnum,values,finite_element);
1563 element->AddInput(WaterfractionEnum,waterfraction,finite_element);
1564 element->AddInput(TemperatureEnum,temperature,finite_element);
[18930]1565
[23644]1566 IssmDouble* n = xNew<IssmDouble>(numnodes);
[23645]1567 if(element->material->ObjectEnum()==MatestarEnum){
[23644]1568 for(i=0;i<numnodes;i++) n[i]=3.;
1569 }
1570 else{
1571 element->GetInputListOnNodes(&n[0],MaterialsRheologyNEnum);
1572 }
1573
[18930]1574 /*Update Rheology only if converged (we must make sure that the temperature is below melting point
1575 * otherwise the rheology could be negative*/
[23644]1576 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
[18930]1577 element->GetInputListOnNodes(&surface[0],SurfaceEnum);
1578 switch(rheology_law){
1579 case NoneEnum:
1580 /*Do nothing: B is not temperature dependent*/
1581 break;
[21377]1582 case BuddJackaEnum:
1583 for(i=0;i<numnodes;i++) B[i]=BuddJacka(temperature[i]);
[25379]1584 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
[21377]1585 break;
[18930]1586 case CuffeyEnum:
1587 for(i=0;i<numnodes;i++) B[i]=Cuffey(temperature[i]);
[25379]1588 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
[18930]1589 break;
[20625]1590 case CuffeyTemperateEnum:
[23644]1591 for(i=0;i<numnodes;i++) B[i]=CuffeyTemperate(temperature[i], waterfraction[i],n[i]);
[25379]1592 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
[20625]1593 break;
[18930]1594 case PatersonEnum:
1595 for(i=0;i<numnodes;i++) B[i]=Paterson(temperature[i]);
[25379]1596 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
[18930]1597 break;
[24060]1598 case NyeH2OEnum:
1599 for(i=0;i<numnodes;i++) B[i]=NyeH2O(values[i]);
[25379]1600 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
[24060]1601 break;
1602 case NyeCO2Enum:
1603 for(i=0;i<numnodes;i++) B[i]=NyeCO2(values[i]);
[25379]1604 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
[24060]1605 break;
[23644]1606 case ArrheniusEnum:{
[18930]1607 element->GetVerticesCoordinates(&xyz_list);
[23644]1608 for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],n[i]);
[25379]1609 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
[18930]1610 break;
[23644]1611 }
1612 case LliboutryDuvalEnum:{
[25539]1613 for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],n[i],element->FindParam(MaterialsBetaEnum),element->FindParam(ConstantsReferencetemperatureEnum),element->FindParam(MaterialsHeatcapacityEnum),element->FindParam(MaterialsLatentheatEnum));
1614 element->AddInput(MaterialsRheologyBEnum,&B[0],finite_element);
1615 break;
[23644]1616 }
[18930]1617 default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
1618 }
[23644]1619 xDelete<IssmDouble>(n);
[16888]1620 }
1621 else{
[25379]1622 element->AddInput(EnthalpyPicardEnum,values,finite_element);
[16888]1623 }
1624
[18930]1625 /*Free ressources:*/
1626 xDelete<IssmDouble>(values);
1627 xDelete<IssmDouble>(pressure);
1628 xDelete<IssmDouble>(surface);
1629 xDelete<IssmDouble>(B);
1630 xDelete<IssmDouble>(temperature);
1631 xDelete<IssmDouble>(waterfraction);
1632 xDelete<IssmDouble>(xyz_list);
1633 xDelete<int>(doflist);
1634}/*}}}*/
1635void EnthalpyAnalysis::PostProcessing(FemModel* femmodel){/*{{{*/
[16888]1636
[18930]1637 /*Intermediaries*/
1638 bool computebasalmeltingrates=true;
[24140]1639 bool isdrainicecolumn;
[18930]1640 IssmDouble dt;
1641
1642 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
[24140]1643 femmodel->parameters->FindParam(&isdrainicecolumn,ThermalIsdrainicecolumnEnum);
[18930]1644
[24335]1645 if(isdrainicecolumn){
1646 DrainWaterfraction(femmodel);
1647 }
1648 if(computebasalmeltingrates){
1649 ComputeBasalMeltingrate(femmodel);
1650 }
[18930]1651
[17027]1652}/*}}}*/
[18930]1653IssmDouble EnthalpyAnalysis::PureIceEnthalpy(Element* element,IssmDouble pressure){/*{{{*/
[16888]1654
[23644]1655 IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
1656 IssmDouble referencetemperature = element->FindParam(ConstantsReferencetemperatureEnum);
[16888]1657
1658 return heatcapacity*(TMeltingPoint(element,pressure)-referencetemperature);
1659}/*}}}*/
[18930]1660IssmDouble EnthalpyAnalysis::TMeltingPoint(Element* element,IssmDouble pressure){/*{{{*/
[16888]1661
[23644]1662 IssmDouble meltingpoint = element->FindParam(MaterialsMeltingpointEnum);
1663 IssmDouble beta = element->FindParam(MaterialsBetaEnum);
[16888]1664
1665 return meltingpoint-beta*pressure;
1666}/*}}}*/
[18930]1667void EnthalpyAnalysis::UpdateBasalConstraints(FemModel* femmodel){/*{{{*/
1668
1669 /*Update basal dirichlet BCs for enthalpy: */
1670 Vector<IssmDouble>* spc = NULL;
1671 IssmDouble* serial_spc = NULL;
1672
[23587]1673 spc=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes());
[18930]1674 /*First create a vector to figure out what elements should be constrained*/
[25539]1675 for(Object* & object : femmodel->elements->objects){
1676 Element* element=xDynamicCast<Element*>(object);
[18930]1677 GetBasalConstraints(spc,element);
1678 }
1679
1680 /*Assemble and serialize*/
1681 spc->Assemble();
1682 serial_spc=spc->ToMPISerial();
1683 delete spc;
1684
1685 /*Then update basal constraints nodes accordingly*/
[25539]1686 for(Object* & object : femmodel->elements->objects){
1687 Element* element=xDynamicCast<Element*>(object);
[18930]1688 ApplyBasalConstraints(serial_spc,element);
1689 }
1690
1691 femmodel->UpdateConstraintsx();
1692
1693 /*Delete*/
1694 xDelete<IssmDouble>(serial_spc);
1695}/*}}}*/
1696void EnthalpyAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
[20453]1697 SetActiveNodesLSMx(femmodel);
[18930]1698}/*}}}*/
Note: See TracBrowser for help on using the repository browser.