Ice Sheet System Model  4.18
Code documentation
GLheightadvectionAnalysis.cpp
Go to the documentation of this file.
2 #include "../toolkits/toolkits.h"
3 #include "../classes/classes.h"
4 #include "../shared/shared.h"
5 #include "../modules/modules.h"
6 
7 /*Model processing*/
9 
10  /*No constraints for now*/
11 }/*}}}*/
13 
14  /*No loads*/
15 }/*}}}*/
16 void GLheightadvectionAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
17 
18  /*First fetch data: */
19  if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
21  iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
22 }/*}}}*/
23 int GLheightadvectionAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
24  return 1;
25 }/*}}}*/
26 void GLheightadvectionAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
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(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
34  counter++;
35  }
36  }
37 
38  if(iomodel->domaintype!=Domain2DhorizontalEnum){
39  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
40  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
41  }
42  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonboundary",MeshVertexonboundaryEnum);
43 }/*}}}*/
44 void GLheightadvectionAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
45 }/*}}}*/
46 
47 /*Finite Element Analysis*/
49  _error_("not implemented");
50 }/*}}}*/
52  /*Default, return NULL*/
53  return NULL;
54 }/*}}}*/
56  _error_("Not implemented");
57 }/*}}}*/
59 
60  /* Spawn basal element */
61  if(!element->IsOnBase()) return NULL;
62  Element* basalelement = element->SpawnBasalElement();
63 
64  /* Check if ice in element */
65  if(!basalelement->IsIceInElement()) return NULL;
66 
67  /*Intermediaries */
68  const IssmPDouble yts = 365*24*3600.;
69  int domaintype,dim;
70  IssmDouble Jdet,D_scalar,onboundary;
71  IssmDouble vel,vx,vy;
72  IssmDouble* xyz_list = NULL;
73  Input2* vx_input = NULL;
74  Input2* vy_input = NULL;
75  Input2* bc_input = NULL;
76 
77  /*Get problem dimension*/
78  basalelement->FindParam(&domaintype,DomainTypeEnum);
79  switch(domaintype){
80  case Domain2DhorizontalEnum: dim = 2; break;
81  case Domain3DEnum: dim = 2; break;
82  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
83  }
84 
85  /*Fetch number of nodes and dof for this finite element*/
86  int numnodes = basalelement->GetNumberOfNodes();
87 
88  /*Initialize Element vector and other vectors*/
89  ElementMatrix* Ke = basalelement->NewElementMatrix();
90  IssmDouble* basis = xNew<IssmDouble>(numnodes);
91  IssmDouble* dbasis = xNew<IssmDouble>(dim*numnodes);
92  IssmDouble* Bprime = xNew<IssmDouble>(dim*numnodes);
93  IssmDouble D[2][2] = {0.};
94 
95  /*Retrieve all inputs and parameters*/
96  basalelement->GetVerticesCoordinates(&xyz_list);
97  switch(domaintype){
99  vx_input=basalelement->GetInput2(VxEnum); _assert_(vx_input);
100  vy_input=basalelement->GetInput2(VyEnum); _assert_(vy_input);
101  bc_input=basalelement->GetInput2(MeshVertexonboundaryEnum); _assert_(bc_input);
102  break;
103  case Domain3DEnum:
104  vx_input=basalelement->GetInput2(VxAverageEnum); _assert_(vx_input);
105  vy_input=basalelement->GetInput2(VyAverageEnum); _assert_(vy_input);
106  bc_input=basalelement->GetInput2(MeshVertexonboundaryEnum); _assert_(bc_input);
107  break;
108  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
109  }
110  IssmDouble h = basalelement->CharacteristicLength();
111 
112  /* Start looping on the number of gaussian points: */
113  Gauss* gauss=basalelement->NewGauss(4);
114  for(int ig=gauss->begin();ig<gauss->end();ig++){
115  gauss->GaussPoint(ig);
116 
117  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
118  basalelement->NodalFunctions(basis,gauss);
119  basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
120  GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
121 
122  /*Get velocity*/
123  vx_input->GetInputValue(&vx,gauss);
124  vy_input->GetInputValue(&vy,gauss);
125 
126  if(false){
127  /*Streamline diffusion*/
128  vel = sqrt(vx*vx+vy*vy);
129  if(vel<10./yts){
130  vx = 0.; vy = 0.;
131  vel = 30./yts*500000.;
132  }
133 
134  for(int i=0;i<numnodes;i++){
135  for(int j=0;j<numnodes;j++){
136  Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
137  (vx*dbasis[0*numnodes+i] + vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j])
138  + vel/500000.*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]));
139  }
140  }
141  }
142  else{
143  D_scalar=gauss->weight*Jdet;
144 
145  bc_input->GetInputValue(&onboundary,gauss);
146  if(onboundary>0.){
147  /*We do not want to advect garbage, make sure only diffusion is applied on boundary*/
148  vx = 0.; vy = 0.;
149  }
150 
151  /*Diffusion */
152  if(sqrt(vx*vx+vy*vy)<1000./31536000.){
153  IssmPDouble kappa = -10.;
154  for(int i=0;i<numnodes;i++){
155  for(int j=0;j<numnodes;j++){
156  Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
157  }
158  }
159  }
160 
161  /*Advection: */
162  for(int i=0;i<numnodes;i++){
163  for(int j=0;j<numnodes;j++){
164  Ke->values[i*numnodes+j] += (D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]))*1e-2;
165  }
166  }
167 
168  /*Artificial diffusivity*/
169  vel=sqrt(vx*vx + vy*vy)+1.e-14;
170  D[0][0]=D_scalar*h/(2.*vel)*fabs(vx*vx); D[0][1]=D_scalar*h/(2.*vel)*fabs(vx*vy);
171  D[1][0]=D_scalar*h/(2.*vel)*fabs(vy*vx); D[1][1]=D_scalar*h/(2.*vel)*fabs(vy*vy);
172  TripleMultiply(Bprime,dim,numnodes,1,
173  &D[0][0],2,2,0,
174  Bprime,dim,numnodes,0,
175  &Ke->values[0],1);
176  }
177  }
178 
179  /*Clean up and return*/
180  xDelete<IssmDouble>(xyz_list);
181  xDelete<IssmDouble>(basis);
182  xDelete<IssmDouble>(dbasis);
183  xDelete<IssmDouble>(Bprime);
184  delete gauss;
185  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
186  return Ke;
187 
188 }/*}}}*/
190  return NULL;
191 }/*}}}*/
193  _error_("not implemented yet");
194 }/*}}}*/
195 void GLheightadvectionAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
196  _error_("Not implemented yet");
197 }/*}}}*/
199 
200  int domaintype;
201  element->FindParam(&domaintype,DomainTypeEnum);
202  switch(domaintype){
205  break;
206  case Domain3DEnum:
208  break;
209  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
210  }
211 }/*}}}*/
213 
214  /*Deal with ocean constraint*/
216 
217  /*Constrain all nodes that are grounded and unconstrain the ones that float*/
218  for(int i=0;i<femmodel->elements->Size();i++){
219  Element *element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
220  int numnodes = element->GetNumberOfNodes();
221  IssmDouble *mask = xNew<IssmDouble>(numnodes);
222  IssmDouble *bed = xNew<IssmDouble>(numnodes);
223  IssmDouble *ls_active = xNew<IssmDouble>(numnodes);
224 
225  element->GetInputListOnNodes(&mask[0],MaskOceanLevelsetEnum);
226  element->GetInputListOnNodes(&bed[0],BaseEnum);
227  element->GetInputListOnNodes(&ls_active[0],IceMaskNodeActivationEnum);
228 
229  for(int in=0;in<numnodes;in++){
230  Node* node=element->GetNode(in);
231  if(mask[in]<0. && ls_active[in]==1.){
232  node->Activate();
233  }
234  else{
235  node->Deactivate();
236  node->ApplyConstraint(0,bed[in]);
237  }
238  }
239  xDelete<IssmDouble>(mask);
240  xDelete<IssmDouble>(bed);
241  xDelete<IssmDouble>(ls_active);
242  }
243 
244  return;
245 }/*}}}*/
246 
247 void GLheightadvectionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
248  /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*2.
249  * For node i, Bi can be expressed in the actual coordinate system
250  * by:
251  * Bi=[ N ]
252  * [ N ]
253  * where N is the finiteelement function for node i.
254  *
255  * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
256  */
257 
258  /*Fetch number of nodes for this finite element*/
259  int numnodes = element->GetNumberOfNodes();
260 
261  /*Get nodal functions*/
262  IssmDouble* basis=xNew<IssmDouble>(numnodes);
263  element->NodalFunctions(basis,gauss);
264 
265  /*Build B: */
266  for(int i=0;i<numnodes;i++){
267  for(int j=0;j<dim;j++){
268  B[numnodes*j+i] = basis[i];
269  }
270  }
271 
272  /*Clean-up*/
273  xDelete<IssmDouble>(basis);
274 }/*}}}*/
275 void GLheightadvectionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
276  /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
277  * For node i, Bi' can be expressed in the actual coordinate system
278  * by:
279  * Bi_prime=[ dN/dx ]
280  * [ dN/dy ]
281  * where N is the finiteelement function for node i.
282  *
283  * We assume B' has been allocated already, of size: 3x(2*numnodes)
284  */
285 
286  /*Fetch number of nodes for this finite element*/
287  int numnodes = element->GetNumberOfNodes();
288 
289  /*Get nodal functions derivatives*/
290  IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
291  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
292 
293  /*Build B': */
294  for(int i=0;i<numnodes;i++){
295  for(int j=0;j<dim;j++){
296  Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
297  }
298  }
299 
300  /*Clean-up*/
301  xDelete<IssmDouble>(dbasis);
302 
303 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
GLheightadvectionAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: GLheightadvectionAnalysis.cpp:51
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
Element::GetInputListOnNodes
void GetInputListOnNodes(IssmDouble *pvalue, int enumtype)
Definition: Element.cpp:1106
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Element::IsOnBase
bool IsOnBase()
Definition: Element.cpp:1984
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
GLheightadvectionAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: GLheightadvectionAnalysis.cpp:44
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
GroundinglineHeightEnum
@ GroundinglineHeightEnum
Definition: EnumDefinitions.h:596
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
MeshVertexonboundaryEnum
@ MeshVertexonboundaryEnum
Definition: EnumDefinitions.h:654
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
Element::SpawnBasalElement
virtual Element * SpawnBasalElement(void)=0
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
Elements
Declaration of Elements class.
Definition: Elements.h:17
TripleMultiply
int TripleMultiply(IssmDouble *a, int nrowa, int ncola, int itrna, IssmDouble *b, int nrowb, int ncolb, int itrnb, IssmDouble *c, int nrowc, int ncolc, int itrnc, IssmDouble *d, int iaddd)
Definition: MatrixUtils.cpp:20
IceMaskNodeActivationEnum
@ IceMaskNodeActivationEnum
Definition: EnumDefinitions.h:627
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
GLheightadvectionAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: GLheightadvectionAnalysis.cpp:55
GLheightadvectionAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: GLheightadvectionAnalysis.cpp:26
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element::DeleteMaterials
void DeleteMaterials(void)
Definition: Element.cpp:429
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
GLheightadvectionAnalysis.h
: header file for generic external result object
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
GLheightadvectionAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: GLheightadvectionAnalysis.cpp:198
VxAverageEnum
@ VxAverageEnum
Definition: EnumDefinitions.h:845
Element::NewGauss
virtual Gauss * NewGauss(void)=0
GLheightadvectionAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: GLheightadvectionAnalysis.cpp:8
Element::InputUpdateFromSolutionOneDof
virtual void InputUpdateFromSolutionOneDof(IssmDouble *solution, int inputenum)=0
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
GLheightadvectionAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: GLheightadvectionAnalysis.cpp:192
GLheightadvectionAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: GLheightadvectionAnalysis.cpp:16
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
Element::GetNode
Node * GetNode(int nodeindex)
Definition: Element.cpp:1207
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
GLheightadvectionAnalysis::Core
void Core(FemModel *femmodel)
Definition: GLheightadvectionAnalysis.cpp:48
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
FemModel::elements
Elements * elements
Definition: FemModel.h:44
GLheightadvectionAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: GLheightadvectionAnalysis.cpp:23
Input2
Definition: Input2.h:18
GLheightadvectionAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: GLheightadvectionAnalysis.cpp:12
FemModel
Definition: FemModel.h:31
GLheightadvectionAnalysis::GetBprime
void GetBprime(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
Definition: GLheightadvectionAnalysis.cpp:275
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
Loads
Declaration of Loads class.
Definition: Loads.h:16
Node
Definition: Node.h:23
GLheightadvectionAnalysisEnum
@ GLheightadvectionAnalysisEnum
Definition: EnumDefinitions.h:1077
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Node::Deactivate
void Deactivate(void)
Definition: Node.cpp:681
SetActiveNodesLSMx
void SetActiveNodesLSMx(FemModel *femmodel)
Definition: SetActiveNodesLSMx.cpp:12
Gauss::begin
virtual int begin(void)=0
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
Element::CharacteristicLength
virtual IssmDouble CharacteristicLength(void)=0
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
MeshVertexonsurfaceEnum
@ MeshVertexonsurfaceEnum
Definition: EnumDefinitions.h:655
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
Element::Update
virtual void Update(Inputs2 *inputs2, int index, IoModel *iomodel, int analysis_counter, int analysis_type, int finite_element)=0
GLheightadvectionAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: GLheightadvectionAnalysis.cpp:189
VyAverageEnum
@ VyAverageEnum
Definition: EnumDefinitions.h:849
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
Node::ApplyConstraint
void ApplyConstraint(int dof, IssmDouble value)
Definition: Node.cpp:646
ElementVector
Definition: ElementVector.h:20
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
IoModel
Definition: IoModel.h:48
GLheightadvectionAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: GLheightadvectionAnalysis.cpp:58
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
Element::InputUpdateFromSolutionOneDofCollapsed
virtual void InputUpdateFromSolutionOneDofCollapsed(IssmDouble *solution, int inputenum)=0
Gauss
Definition: Gauss.h:8
GLheightadvectionAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: GLheightadvectionAnalysis.cpp:212
GLheightadvectionAnalysis::GetB
void GetB(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
Definition: GLheightadvectionAnalysis.cpp:247
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
Node::Activate
void Activate(void)
Definition: Node.cpp:632
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
GLheightadvectionAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: GLheightadvectionAnalysis.cpp:195
Element::NewElementMatrix
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2497
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16