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

Last change on this file since 16542 was 16542, checked in by Mathieu Morlighem, 11 years ago

NEW: create all datasets once for all (nodes, elements, constraints, parameters), do not pass pointers anymore (much simpler)
File size: 5.8 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}/*}}}*/
13void EnthalpyAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
14
15 bool dakota_analysis;
16 bool isenthalpy;
17
18 /*Now, is the model 3d? otherwise, do nothing: */
19 if(iomodel->meshtype==Mesh2DhorizontalEnum)return;
20
21 /*Is enthalpy requested?*/
22 iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
23 if(!isenthalpy) return;
24
25 /*Fetch data needed: */
26 iomodel->FetchData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
27
28 /*Update elements: */
29 int counter=0;
30 for(int i=0;i<iomodel->numberofelements;i++){
31 if(iomodel->my_elements[i]){
32 Element* element=(Element*)elements->GetObjectByOffset(counter);
33 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
34 counter++;
35 }
36 }
37
38 iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
39
40 iomodel->FetchDataToInput(elements,ThicknessEnum);
41 iomodel->FetchDataToInput(elements,SurfaceEnum);
42 iomodel->FetchDataToInput(elements,BedEnum);
43 iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
44 iomodel->FetchDataToInput(elements,FrictionPEnum);
45 iomodel->FetchDataToInput(elements,FrictionQEnum);
46 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
47 iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
48 iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
49 iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
50 iomodel->FetchDataToInput(elements,FlowequationElementEquationEnum);
51 iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
52 iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
53 iomodel->FetchDataToInput(elements,PressureEnum);
54 iomodel->FetchDataToInput(elements,TemperatureEnum);
55 iomodel->FetchDataToInput(elements,WaterfractionEnum);
56 iomodel->FetchDataToInput(elements,EnthalpyEnum);
57 iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum);
58 iomodel->FetchDataToInput(elements,WatercolumnEnum);
59 iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
60 iomodel->FetchDataToInput(elements,VxEnum);
61 iomodel->FetchDataToInput(elements,VyEnum);
62 iomodel->FetchDataToInput(elements,VzEnum);
63 InputUpdateFromConstantx(elements,0.,VxMeshEnum);
64 InputUpdateFromConstantx(elements,0.,VyMeshEnum);
65 InputUpdateFromConstantx(elements,0.,VzMeshEnum);
66 if(dakota_analysis){
67 elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum);
68 elements->InputDuplicate(BasalforcingsMeltingRateEnum,QmuMeltingEnum);
69 elements->InputDuplicate(VxMeshEnum,QmuVxMeshEnum);
70 elements->InputDuplicate(VxMeshEnum,QmuVyMeshEnum);
71 elements->InputDuplicate(VxMeshEnum,QmuVzMeshEnum);
72 }
73
74 /*Free data: */
75 iomodel->DeleteData(3,TemperatureEnum,WaterfractionEnum,PressureEnum);
76}/*}}}*/
77void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
78
79 if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
80 ::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
81 iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
82}/*}}}*/
83void EnthalpyAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
84
85 /*Intermediary*/
86 int count;
87 int M,N;
88 bool spcpresent = false;
89 IssmDouble heatcapacity;
90 IssmDouble referencetemperature;
91
92 /*Output*/
93 IssmDouble *spcvector = NULL;
94 IssmDouble* times=NULL;
95 IssmDouble* values=NULL;
96
97 /*Fetch parameters: */
98 iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
99 iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
100
101 /*return if 2d mesh*/
102 if(iomodel->meshtype==Mesh2DhorizontalEnum) return;
103
104 /*Fetch data: */
105 iomodel->FetchData(&spcvector,&M,&N,ThermalSpctemperatureEnum);
106
107 //FIX ME: SHOULD USE IOMODELCREATECONSTRAINTS
108 /*Transient or static?:*/
109 if(M==iomodel->numberofvertices){
110 /*static: just create Constraints objects*/
111 count=0;
112
113 for(int i=0;i<iomodel->numberofvertices;i++){
114 /*keep only this partition's nodes:*/
115 if((iomodel->my_vertices[i])){
116
117 if (!xIsNan<IssmDouble>(spcvector[i])){
118
119 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum));
120 count++;
121
122 }
123 }
124 }
125 }
126 else if (M==(iomodel->numberofvertices+1)){
127 /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve
128 * various times and values to initialize an SpcTransient object: */
129 count=0;
130
131 /*figure out times: */
132 times=xNew<IssmDouble>(N);
133 for(int j=0;j<N;j++){
134 times[j]=spcvector[(M-1)*N+j];
135 }
136
137 /*Create constraints from x,y,z: */
138 for(int i=0;i<iomodel->numberofvertices;i++){
139
140 /*keep only this partition's nodes:*/
141 if((iomodel->my_vertices[i])){
142
143 /*figure out times and values: */
144 values=xNew<IssmDouble>(N);
145 spcpresent=false;
146 for(int j=0;j<N;j++){
147 values[j]=heatcapacity*(spcvector[i*N+j]-referencetemperature);
148 if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
149 }
150
151 if(spcpresent){
152 constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,EnthalpyAnalysisEnum));
153 count++;
154 }
155 xDelete<IssmDouble>(values);
156 }
157 }
158 }
159 else{
160 _error_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported");
161 }
162
163 /*Free ressources:*/
164 iomodel->DeleteData(spcvector,ThermalSpctemperatureEnum);
165 xDelete<IssmDouble>(times);
166 xDelete<IssmDouble>(values);
167}/*}}}*/
168void EnthalpyAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
169
170 /*No loads */
171}/*}}}*/
Note: See TracBrowser for help on using the repository browser.