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

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

CHG: adding Enthalpy GetInputFromSolution
File size: 10.3 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"
6
7/*Model processing*/
[16539]8int EnthalpyAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
[16534]9 return 1;
10}/*}}}*/
[16542]11void EnthalpyAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
[16604]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);
[16539]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}/*}}}*/
[16542]89void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
[16539]90
91 if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
[16542]92 ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
[16539]93 iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
94}/*}}}*/
[16542]95void EnthalpyAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
[16539]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}/*}}}*/
[16542]180void EnthalpyAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
[16539]181
182 /*No loads */
183}/*}}}*/
[16675]184
185/*Numerics*/
186void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
187 element->GetSolutionFromInputsOneDof(solution,EnthalpyEnum);
188}/*}}}*/
[16684]189void EnthalpyAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
[16734]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* temperature = xNew<IssmDouble>(numnodes);
205 IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
206
207 /*Use the dof list to index into the solution vector: */
208 for(i=0;i<numnodes;i++){
209 values[i]=solution[doflist[i]];
210
211 /*Check solution*/
212 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
213 }
214
215 /*Get all inputs and parameters*/
216 element->GetInputValue(&converged,ConvergedEnum);
217 element->GetInputListOnNodes(pressure,PressureEnum);
218 if(converged){
219 for(i=0;i<numnodes;i++){
220 element->EnthalpyToThermal(&temperature[i],&waterfraction[i],values[i],pressure[i]);
221 if(waterfraction[i]<0.) _error_("Negative water fraction found in solution vector");
222 //if(waterfraction[i]>1.) _error_("Water fraction >1 found in solution vector");
223 }
224 element->AddInput(EnthalpyEnum,values,P1Enum);
225 element->AddInput(WaterfractionEnum,waterfraction,P1Enum);
226 element->AddInput(TemperatureEnum,temperature,P1Enum);
227
228 /*Update Rheology only if converged (we must make sure that the temperature is below melting point
229 * otherwise the rheology could be negative*/
230 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
231 switch(rheology_law){
232 case NoneEnum:
233 /*Do nothing: B is not temperature dependent*/
234 break;
235 case PatersonEnum:
236 for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
237 B_average=Paterson(T_average);
238 element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
239 break;
240 case ArrheniusEnum:{
241 Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
242 surface_input->GetInputAverage(&s_average);
243 element->GetVerticesCoordinates(&xyz_list);
244 for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
245 //B_average=Arrhenius(T_average,
246 //s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0),
247 //element->GetMaticeParameter(MaterialsRheologyNEnum));
248 element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
249 break;
250 }
251 case LliboutryDuvalEnum:
252 for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
253 for(i=0;i<numnodes;i++) P_average+=pressure[i]/reCast<IssmDouble>(numnodes);
254 B_average=LliboutryDuval(T_average,P_average,
255 element->GetMaterialParameter(MaterialsRheologyNEnum),
256 element->GetMaterialParameter(MaterialsBetaEnum),
257 element->GetMaterialParameter(ConstantsReferencetemperatureEnum),
258 element->GetMaterialParameter(MaterialsHeatcapacityEnum),
259 element->GetMaterialParameter(MaterialsLatentheatEnum));
260 element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
261 break;
262 default:
263 _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
264 }
265 }
266 else{
267 element->AddInput(EnthalpyPicardEnum,values,P1Enum);
268 }
269
270 /*Free ressources:*/
271 xDelete<IssmDouble>(values);
272 xDelete<IssmDouble>(pressure);
273 xDelete<IssmDouble>(temperature);
274 xDelete<IssmDouble>(waterfraction);
275 xDelete<IssmDouble>(xyz_list);
276 xDelete<int>(doflist);
[16684]277}/*}}}*/
Note: See TracBrowser for help on using the repository browser.