Ice Sheet System Model  4.18
Code documentation
Functions
SurfaceRelVelMisfitx.h File Reference

header file for inverse methods misfit computation More...

#include "../../classes/classes.h"

Go to the source code of this file.

Functions

void SurfaceRelVelMisfitx (IssmDouble *pJ, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters)
 
IssmDouble SurfaceRelVelMisfit (Element *element)
 

Detailed Description

header file for inverse methods misfit computation

Definition in file SurfaceRelVelMisfitx.h.

Function Documentation

◆ SurfaceRelVelMisfitx()

void SurfaceRelVelMisfitx ( IssmDouble pJ,
Elements elements,
Nodes nodes,
Vertices vertices,
Loads loads,
Materials materials,
Parameters parameters 
)

Definition at line 11 of file SurfaceRelVelMisfitx.cpp.

11  {
12 
13  /*output: */
14  IssmDouble J=0.;
15  IssmDouble J_sum;
16 
17  /*Compute Misfit: */
18  for(int i=0;i<elements->Size();i++){
19  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
20  J+=SurfaceRelVelMisfit(element);
21  }
22 
23  /*Sum all J from all cpus of the cluster:*/
26  J=J_sum;
27 
28  /*Assign output pointers: */
29  *pJ=J;
30 }

◆ SurfaceRelVelMisfit()

IssmDouble SurfaceRelVelMisfit ( Element element)

Definition at line 32 of file SurfaceRelVelMisfitx.cpp.

32  {
33 
34  int domaintype,numcomponents;
35  IssmDouble Jelem=0.;
36  IssmDouble misfit,Jdet,scalex,scaley;
37  IssmDouble epsvel=2.220446049250313e-16;
38  IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
39  IssmDouble vx,vy,vxobs,vyobs,weight;
40  IssmDouble* xyz_list = NULL;
41 
42  /*Get basal element*/
43  if(!element->IsOnSurface()) return 0.;
44 
45  /*If on water, return 0: */
46  if(!element->IsIceInElement()) return 0.;
47 
48  /*Get problem dimension*/
49  element->FindParam(&domaintype,DomainTypeEnum);
50  switch(domaintype){
51  case Domain2DverticalEnum: numcomponents = 1; break;
52  case Domain3DEnum: numcomponents = 2; break;
53  case Domain2DhorizontalEnum: numcomponents = 2; break;
54  default: _error_("not supported yet");
55  }
56 
57  /*Spawn surface element*/
58  Element* topelement = element->SpawnTopElement();
59 
60  /* Get node coordinates*/
61  topelement->GetVerticesCoordinates(&xyz_list);
62 
63  /*Retrieve all inputs we will be needing: */
64  DatasetInput2* weights_input=topelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
65  Input2* vx_input =topelement->GetInput2(VxEnum); _assert_(vx_input);
66  Input2* vxobs_input =topelement->GetInput2(InversionVxObsEnum); _assert_(vxobs_input);
67  Input2* vy_input = NULL;
68  Input2* vyobs_input = NULL;
69  if(numcomponents==2){
70  vy_input =topelement->GetInput2(VyEnum); _assert_(vy_input);
71  vyobs_input =topelement->GetInput2(InversionVyObsEnum); _assert_(vyobs_input);
72  }
73 
74  /* Start looping on the number of gaussian points: */
75  Gauss* gauss=topelement->NewGauss(4);
76  for(int ig=gauss->begin();ig<gauss->end();ig++){
77 
78  gauss->GaussPoint(ig);
79 
80  /* Get Jacobian determinant: */
81  topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
82 
83  /*Get all parameters at gaussian point*/
84  weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum);
85  vx_input->GetInputValue(&vx,gauss);
86  vxobs_input->GetInputValue(&vxobs,gauss);
87  if(numcomponents==2){
88  vy_input->GetInputValue(&vy,gauss);
89  vyobs_input->GetInputValue(&vyobs,gauss);
90  }
91 
92  /*Compute SurfaceRelVelMisfit:
93  *
94  * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
95  * J = --- | ------------- (u - u ) + ------------- (v - v ) |
96  * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
97  * obs obs
98  */
99 
100  if(numcomponents==2){
101  scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
102  scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
103  misfit=0.5*(scalex*pow((vx-vxobs),2)+scaley*pow((vy-vyobs),2));
104  }
105  else{
106  scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
107  misfit=0.5*(scalex*pow((vx-vxobs),2));
108  }
109 
110  /*Add to cost function*/
111  Jelem+=misfit*weight*Jdet*gauss->weight;
112  }
113 
114  /*clean up and Return: */
115  if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
116  xDelete<IssmDouble>(xyz_list);
117  delete gauss;
118  return Jelem;
119 }
DataSet::Size
int Size()
Definition: DataSet.cpp:399
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
DatasetInput2
Definition: DatasetInput2.h:14
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
InversionVyObsEnum
@ InversionVyObsEnum
Definition: EnumDefinitions.h:634
InversionCostFunctionsCoefficientsEnum
@ InversionCostFunctionsCoefficientsEnum
Definition: EnumDefinitions.h:629
ISSM_MPI_SUM
#define ISSM_MPI_SUM
Definition: issmmpi.h:134
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
SurfaceRelVelMisfitEnum
@ SurfaceRelVelMisfitEnum
Definition: EnumDefinitions.h:828
DatasetInput2::GetInputValue
void GetInputValue(IssmDouble *pvalue, Gauss *gauss, int index)
Definition: DatasetInput2.cpp:199
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element::DeleteMaterials
void DeleteMaterials(void)
Definition: Element.cpp:429
Element
Definition: Element.h:41
Element::GetDatasetInput2
virtual DatasetInput2 * GetDatasetInput2(int inputenum)
Definition: Element.h:250
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
ISSM_MPI_DOUBLE
#define ISSM_MPI_DOUBLE
Definition: issmmpi.h:125
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
Element::NewGauss
virtual Gauss * NewGauss(void)=0
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
ISSM_MPI_Bcast
int ISSM_MPI_Bcast(void *buffer, int count, ISSM_MPI_Datatype datatype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:162
Input2
Definition: Input2.h:18
SurfaceRelVelMisfit
IssmDouble SurfaceRelVelMisfit(Element *element)
Definition: SurfaceRelVelMisfitx.cpp:32
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Element::SpawnTopElement
virtual Element * SpawnTopElement(void)=0
Gauss::begin
virtual int begin(void)=0
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
InversionVxObsEnum
@ InversionVxObsEnum
Definition: EnumDefinitions.h:633
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
Element::IsOnSurface
bool IsOnSurface()
Definition: Element.cpp:1981
ISSM_MPI_Reduce
int ISSM_MPI_Reduce(void *sendbuf, void *recvbuf, int count, ISSM_MPI_Datatype datatype, ISSM_MPI_Op op, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:373
Domain2DverticalEnum
@ Domain2DverticalEnum
Definition: EnumDefinitions.h:535
Gauss
Definition: Gauss.h:8