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

Last change on this file since 16748 was 16748, checked in by seroussi, 11 years ago

BUG: fixed enthalpy analysis solution
File size: 9.9 KB
Line 
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
7/*Model processing*/
8int EnthalpyAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
9 return 1;
10}/*}}}*/
11void EnthalpyAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
12
13 int numoutputs;
14 char** requestedoutputs = NULL;
15
16 parameters->AddObject(iomodel->CopyConstantObject(ThermalStabilizationEnum));
17 parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
18 parameters->AddObject(iomodel->CopyConstantObject(ThermalIsdynamicbasalspcEnum));
19
20 iomodel->FetchData(&requestedoutputs,&numoutputs,ThermalRequestedOutputsEnum);
21 parameters->AddObject(new IntParam(ThermalNumRequestedOutputsEnum,numoutputs));
22 if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs));
23 iomodel->DeleteData(&requestedoutputs,numoutputs,ThermalRequestedOutputsEnum);
24}/*}}}*/
25void EnthalpyAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
26
27 bool dakota_analysis;
28 bool isenthalpy;
29
30 /*Now, is the model 3d? otherwise, do nothing: */
31 if(iomodel->meshtype==Mesh2DhorizontalEnum)return;
32
33 /*Is enthalpy requested?*/
34 iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
35 if(!isenthalpy) return;
36
37 /*Fetch data needed: */
38 iomodel->FetchData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
39
40 /*Update elements: */
41 int counter=0;
42 for(int i=0;i<iomodel->numberofelements;i++){
43 if(iomodel->my_elements[i]){
44 Element* element=(Element*)elements->GetObjectByOffset(counter);
45 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
46 counter++;
47 }
48 }
49
50 iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
51
52 iomodel->FetchDataToInput(elements,ThicknessEnum);
53 iomodel->FetchDataToInput(elements,SurfaceEnum);
54 iomodel->FetchDataToInput(elements,BedEnum);
55 iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
56 iomodel->FetchDataToInput(elements,FrictionPEnum);
57 iomodel->FetchDataToInput(elements,FrictionQEnum);
58 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
59 iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
60 iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
61 iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
62 iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum);
63 iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
64 iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
65 iomodel->FetchDataToInput(elements,PressureEnum);
66 iomodel->FetchDataToInput(elements,TemperatureEnum);
67 iomodel->FetchDataToInput(elements,WaterfractionEnum);
68 iomodel->FetchDataToInput(elements,EnthalpyEnum);
69 iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum);
70 iomodel->FetchDataToInput(elements,WatercolumnEnum);
71 iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
72 iomodel->FetchDataToInput(elements,VxEnum);
73 iomodel->FetchDataToInput(elements,VyEnum);
74 iomodel->FetchDataToInput(elements,VzEnum);
75 InputUpdateFromConstantx(elements,0.,VxMeshEnum);
76 InputUpdateFromConstantx(elements,0.,VyMeshEnum);
77 InputUpdateFromConstantx(elements,0.,VzMeshEnum);
78 if(dakota_analysis){
79 elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum);
80 elements->InputDuplicate(BasalforcingsMeltingRateEnum,QmuMeltingEnum);
81 elements->InputDuplicate(VxMeshEnum,QmuVxMeshEnum);
82 elements->InputDuplicate(VxMeshEnum,QmuVyMeshEnum);
83 elements->InputDuplicate(VxMeshEnum,QmuVzMeshEnum);
84 }
85
86 /*Free data: */
87 iomodel->DeleteData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
88}/*}}}*/
89void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
90
91 if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
92 ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
93 iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
94}/*}}}*/
95void EnthalpyAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
96
97 /*Intermediary*/
98 int count;
99 int M,N;
100 bool spcpresent = false;
101 IssmDouble heatcapacity;
102 IssmDouble referencetemperature;
103
104 /*Output*/
105 IssmDouble *spcvector = NULL;
106 IssmDouble* times=NULL;
107 IssmDouble* values=NULL;
108
109 /*Fetch parameters: */
110 iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
111 iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
112
113 /*return if 2d mesh*/
114 if(iomodel->meshtype==Mesh2DhorizontalEnum) return;
115
116 /*Fetch data: */
117 iomodel->FetchData(&spcvector,&M,&N,ThermalSpctemperatureEnum);
118
119 //FIX ME: SHOULD USE IOMODELCREATECONSTRAINTS
120 /*Transient or static?:*/
121 if(M==iomodel->numberofvertices){
122 /*static: just create Constraints objects*/
123 count=0;
124
125 for(int i=0;i<iomodel->numberofvertices;i++){
126 /*keep only this partition's nodes:*/
127 if((iomodel->my_vertices[i])){
128
129 if (!xIsNan<IssmDouble>(spcvector[i])){
130
131 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum));
132 count++;
133
134 }
135 }
136 }
137 }
138 else if (M==(iomodel->numberofvertices+1)){
139 /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve
140 * various times and values to initialize an SpcTransient object: */
141 count=0;
142
143 /*figure out times: */
144 times=xNew<IssmDouble>(N);
145 for(int j=0;j<N;j++){
146 times[j]=spcvector[(M-1)*N+j];
147 }
148
149 /*Create constraints from x,y,z: */
150 for(int i=0;i<iomodel->numberofvertices;i++){
151
152 /*keep only this partition's nodes:*/
153 if((iomodel->my_vertices[i])){
154
155 /*figure out times and values: */
156 values=xNew<IssmDouble>(N);
157 spcpresent=false;
158 for(int j=0;j<N;j++){
159 values[j]=heatcapacity*(spcvector[i*N+j]-referencetemperature);
160 if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
161 }
162
163 if(spcpresent){
164 constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,EnthalpyAnalysisEnum));
165 count++;
166 }
167 xDelete<IssmDouble>(values);
168 }
169 }
170 }
171 else{
172 _error_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported");
173 }
174
175 /*Free ressources:*/
176 iomodel->DeleteData(spcvector,ThermalSpctemperatureEnum);
177 xDelete<IssmDouble>(times);
178 xDelete<IssmDouble>(values);
179}/*}}}*/
180void EnthalpyAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
181
182 /*No loads */
183}/*}}}*/
184
185/*Numerics*/
186void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
187 element->GetSolutionFromInputsOneDof(solution,EnthalpyEnum);
188}/*}}}*/
189void EnthalpyAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
190
191 bool converged;
192 int i,rheology_law;
193 IssmDouble B_average,s_average,T_average=0.,P_average=0.;
194 int *doflist = NULL;
195 IssmDouble *xyz_list = NULL;
196
197 /*Fetch number of nodes and dof for this finite element*/
198 int numnodes = element->GetNumberOfNodes();
199
200 /*Fetch dof list and allocate solution vector*/
201 element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
202 IssmDouble* values = xNew<IssmDouble>(numnodes);
203 IssmDouble* pressure = xNew<IssmDouble>(numnodes);
204 IssmDouble* surface = xNew<IssmDouble>(numnodes);
205 IssmDouble* B = xNew<IssmDouble>(numnodes);
206 IssmDouble* temperature = xNew<IssmDouble>(numnodes);
207 IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
208
209 /*Use the dof list to index into the solution vector: */
210 for(i=0;i<numnodes;i++){
211 values[i]=solution[doflist[i]];
212
213 /*Check solution*/
214 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
215 }
216
217 /*Get all inputs and parameters*/
218 element->GetInputValue(&converged,ConvergedEnum);
219 element->GetInputListOnNodes(&pressure[0],PressureEnum);
220 if(converged){
221 for(i=0;i<numnodes;i++){
222 element->EnthalpyToThermal(&temperature[i],&waterfraction[i],values[i],pressure[i]);
223 if(waterfraction[i]<0.) _error_("Negative water fraction found in solution vector");
224 //if(waterfraction[i]>1.) _error_("Water fraction >1 found in solution vector");
225 }
226 element->AddInput(EnthalpyEnum,values,P1Enum);
227 element->AddInput(WaterfractionEnum,waterfraction,P1Enum);
228 element->AddInput(TemperatureEnum,temperature,P1Enum);
229
230 /*Update Rheology only if converged (we must make sure that the temperature is below melting point
231 * otherwise the rheology could be negative*/
232 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
233 element->GetInputListOnNodes(&surface[0],SurfaceEnum);
234 switch(rheology_law){
235 case NoneEnum:
236 /*Do nothing: B is not temperature dependent*/
237 break;
238 case PatersonEnum:
239 for(i=0;i<numnodes;i++) B[i]=Paterson(temperature[i]);
240 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
241 break;
242 case ArrheniusEnum:
243 element->GetVerticesCoordinates(&xyz_list);
244 for(i=0;i<numnodes;i++) B[i]=Arrhenius(temperature[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
245 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
246 break;
247 case LliboutryDuvalEnum:
248 //for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],material->GetN(),matpar->GetBeta(),matpar->GetReferenceTemperature(),matpar->GetHeatCapacity(),matpar->GetLatentHeat()); for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]);
249 for(i=0;i<numnodes;i++) B[i]=Paterson(temperature[i]);
250 element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
251 break;
252 default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
253 }
254 }
255 else{
256 element->AddInput(EnthalpyPicardEnum,values,P1Enum);
257 }
258
259 /*Free ressources:*/
260 xDelete<IssmDouble>(values);
261 xDelete<IssmDouble>(pressure);
262 xDelete<IssmDouble>(surface);
263 xDelete<IssmDouble>(B);
264 xDelete<IssmDouble>(temperature);
265 xDelete<IssmDouble>(waterfraction);
266 xDelete<IssmDouble>(xyz_list);
267 xDelete<int>(doflist);
268}/*}}}*/
Note: See TracBrowser for help on using the repository browser.