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