Ice Sheet System Model  4.18
Code documentation
Cfsurfacelogvel.cpp
Go to the documentation of this file.
1 
5 /*Headers:*/
6 /*{{{*/
7 #ifdef HAVE_CONFIG_H
8  #include <config.h>
9 #else
10 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
11 #endif
12 
13 #include "./classes.h"
16 #include "../datastructures/datastructures.h"
17 #include "./Elements/Element.h"
18 #include "./Elements/Elements.h"
19 #include "./FemModel.h"
20 #include "../modules/SurfaceAreax/SurfaceAreax.h"
21 #include "../classes/Params/Parameters.h"
22 #include "../classes/gauss/Gauss.h"
24 /*}}}*/
25 
26 /*Cfsurfacelogvel constructors, destructors :*/
28 
29  this->definitionenum = -1;
30  this->name = NULL;
31  this->misfit=0;
32  this->lock=0;
33  this->datatime=0.;
34  this->timepassedflag = false;
35  this->last_time=0.;
36 
37 }
38 /*}}}*/
39 Cfsurfacelogvel::Cfsurfacelogvel(char* in_name, int in_definitionenum, IssmDouble in_datatime, bool in_timepassedflag){/*{{{*/
40 
41  this->definitionenum=in_definitionenum;
42 
43  this->name = xNew<char>(strlen(in_name)+1);
44  xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
45 
46  this->datatime=in_datatime;
47  this->timepassedflag=in_timepassedflag;
48 
49  this->misfit=0;
50  this->lock=0;
51 }
52 /*}}}*/
54  if(this->name)xDelete(this->name);
55  this->misfit=0;
56  this->lock=0;
57 }
58 /*}}}*/
59 /*Object virtual function resolutoin: */
61  Cfsurfacelogvel* mf = new Cfsurfacelogvel(this->name,this->definitionenum,this->datatime,this->timepassedflag);
62  mf->misfit=this->misfit;
63  mf->lock=this->lock;
64  return (Object*) mf;
65 }
66 /*}}}*/
67 void Cfsurfacelogvel::DeepEcho(void){/*{{{*/
68  this->Echo();
69 }
70 /*}}}*/
71 void Cfsurfacelogvel::Echo(void){/*{{{*/
72  _printf_(" Cfsurfacelogvel: " << name << " " << this->definitionenum << "\n");
73  _printf_(" datatime: " << datatime << "\n");
74  _printf_(" timepassedflag: "<<timepassedflag<<"\n");
75 }
76 /*}}}*/
77 int Cfsurfacelogvel::Id(void){/*{{{*/
78  return -1;
79 }
80 /*}}}*/
81 void Cfsurfacelogvel::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/
82  _error_("not implemented yet!");
83 }
84 /*}}}*/
85 int Cfsurfacelogvel::ObjectEnum(void){/*{{{*/
86  return CfsurfacelogvelEnum;
87 }
88 /*}}}*/
89 /*Definition virtual function resolutoin: */
91  return this->definitionenum;
92 }
93 /*}}}*/
94 char* Cfsurfacelogvel::Name(){/*{{{*/
95  char* name2=xNew<char>(strlen(this->name)+1);
96  xMemCpy(name2,this->name,strlen(this->name)+1);
97 
98  return name2;
99 }
100 /*}}}*/
102 
103  /*diverse: */
104  IssmDouble time;
105 
106  /*recover time parameters: */
108  if(time < last_time) timepassedflag = false;
109  last_time = time;
110 
111  int i;
112  IssmDouble J=0.;
113  IssmDouble J_sum=0.;
114 
115  if(datatime<=time && !timepassedflag){
116  for(i=0;i<femmodel->elements->Size();i++){
119  }
120 
121  ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
123  J=J_sum;
124 
125  timepassedflag = true;
126  return J;
127  }
128  else return J;
129  }
130  /*}}}*/
132 
133  int domaintype,numcomponents;
134  IssmDouble Jelem=0.;
135  IssmDouble epsvel=2.220446049250313e-16;
136  IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
137  IssmDouble velocity_mag,obs_velocity_mag;
138  IssmDouble misfit,Jdet;
139  IssmDouble vx,vy,vxobs,vyobs,weight;
140  IssmDouble* xyz_list = NULL;
141 
142  /*Get basal element*/
143  if(!element->IsOnSurface()) return 0.;
144 
145  /*If on water, return 0: */
146  if(!element->IsIceInElement()) return 0.;
147 
148  /*Get problem dimension*/
149  element->FindParam(&domaintype,DomainTypeEnum);
150  switch(domaintype){
151  case Domain2DverticalEnum: numcomponents = 1; break;
152  case Domain3DEnum: numcomponents = 2; break;
153  case Domain2DhorizontalEnum: numcomponents = 2; break;
154  default: _error_("not supported yet");
155  }
156 
157  /*Spawn surface element*/
158  Element* topelement = element->SpawnTopElement();
159 
160  /* Get node coordinates*/
161  topelement->GetVerticesCoordinates(&xyz_list);
162 
163  /*Get model values*/
164  Input2 *vx_input = topelement->GetInput2(VxEnum); _assert_(vx_input);
165  Input2 *vy_input = NULL;
166  if(numcomponents==2){
167  vy_input = topelement->GetInput2(VyEnum); _assert_(vy_input);
168  }
169 
170  /*Retrieve all inputs we will be needing: */
171  DatasetInput2 *datasetinput = topelement->GetDatasetInput2(definitionenum); _assert_(datasetinput);
172 
173  /* Start looping on the number of gaussian points: */
174  Gauss* gauss=topelement->NewGauss(2);
175  for(int ig=gauss->begin();ig<gauss->end();ig++){
176 
177  gauss->GaussPoint(ig);
178 
179  /* Get Jacobian determinant: */
180  topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
181 
182  /*Get all parameters at gaussian point*/
183  datasetinput->GetInputValue(&weight,gauss,WeightsSurfaceObservationEnum);
184  vx_input->GetInputValue(&vx,gauss);
185  datasetinput->GetInputValue(&vxobs,gauss,VxObsEnum);
186  if(numcomponents==2){
187  vy_input->GetInputValue(&vy,gauss);
188  datasetinput->GetInputValue(&vyobs,gauss,VyObsEnum);
189  }
190 
191  /*Compute SurfaceLogVelMisfit:
192  * * [ vel + eps ] 2
193  * * J = 4 \bar{v}^2 | log ( ----------- ) |
194  * * [ vel + eps ]
195  * * obs
196  * */
197  if(numcomponents==1){
198  velocity_mag =fabs(vx)+epsvel;
199  obs_velocity_mag=fabs(vxobs)+epsvel;
200  }
201  else{
202  velocity_mag =sqrt(vx*vx+vy*vy)+epsvel;
203  obs_velocity_mag=sqrt(vxobs*vxobs+vyobs*vyobs)+epsvel;
204  }
205 
206  misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
207 
208  /*Add to cost function*/
209  Jelem+=misfit*weight*Jdet*gauss->weight;
210 
211  }
212 
213  /*clean up and Return: */
214  if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
215  xDelete<IssmDouble>(xyz_list);
216  delete gauss;
217  return Jelem;
218 }/*}}}*/
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
ISSM_MPI_Allreduce
int ISSM_MPI_Allreduce(void *sendbuf, void *recvbuf, int count, ISSM_MPI_Datatype datatype, ISSM_MPI_Op op, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:116
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
Cfsurfacelogvel::last_time
IssmDouble last_time
Definition: Cfsurfacelogvel.h:23
ISSM_MPI_SUM
#define ISSM_MPI_SUM
Definition: issmmpi.h:134
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
Cfsurfacelogvel::~Cfsurfacelogvel
~Cfsurfacelogvel()
Definition: Cfsurfacelogvel.cpp:53
Cfsurfacelogvel::DefinitionEnum
int DefinitionEnum()
Definition: Cfsurfacelogvel.cpp:90
CfsurfacelogvelEnum
@ CfsurfacelogvelEnum
Definition: EnumDefinitions.h:1006
DatasetInput2::GetInputValue
void GetInputValue(IssmDouble *pvalue, Gauss *gauss, int index)
Definition: DatasetInput2.cpp:199
Cfsurfacelogvel::Id
int Id(void)
Definition: Cfsurfacelogvel.cpp:77
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
Cfsurfacelogvel::name
char * name
Definition: Cfsurfacelogvel.h:20
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element::DeleteMaterials
void DeleteMaterials(void)
Definition: Element.cpp:429
Element
Definition: Element.h:41
Cfsurfacelogvel::misfit
IssmDouble misfit
Definition: Cfsurfacelogvel.h:26
Elements.h
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
Element.h
abstract class for Element object This class is a place holder for the Tria and the Penta elements....
Object
Definition: Object.h:13
VyObsEnum
@ VyObsEnum
Definition: EnumDefinitions.h:852
Cfsurfacelogvel::datatime
IssmDouble datatime
Definition: Cfsurfacelogvel.h:21
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
xDelete
void xDelete(T **&aT_pp)
Definition: MemOps.h:97
Element::NewGauss
virtual Gauss * NewGauss(void)=0
ExternalResult.h
abstract class for ExternalResult object
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
Cfsurfacelogvel::Response
IssmDouble Response(FemModel *femmodel)
Definition: Cfsurfacelogvel.cpp:101
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
Cfsurfacelogvel::definitionenum
int definitionenum
Definition: Cfsurfacelogvel.h:19
Cfsurfacelogvel::Name
char * Name()
Definition: Cfsurfacelogvel.cpp:94
FemModel::elements
Elements * elements
Definition: FemModel.h:44
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
Cfsurfacelogvel::copy
Object * copy()
Definition: Cfsurfacelogvel.cpp:60
FemModel
Definition: FemModel.h:31
Cfsurfacelogvel
Definition: Cfsurfacelogvel.h:15
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Cfsurfacelogvel::timepassedflag
bool timepassedflag
Definition: Cfsurfacelogvel.h:22
Element::SpawnTopElement
virtual Element * SpawnTopElement(void)=0
Cfsurfacelogvel::ObjectEnum
int ObjectEnum(void)
Definition: Cfsurfacelogvel.cpp:85
Gauss::begin
virtual int begin(void)=0
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
FemModel.h
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
WeightsSurfaceObservationEnum
@ WeightsSurfaceObservationEnum
Definition: EnumDefinitions.h:864
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
Element::IsOnSurface
bool IsOnSurface()
Definition: Element.cpp:1981
xMemCpy
T * xMemCpy(T *dest, const T *src, unsigned int size)
Definition: MemOps.h:152
Cfsurfacelogvel::Marshall
void Marshall(char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
Definition: Cfsurfacelogvel.cpp:81
Domain2DverticalEnum
@ Domain2DverticalEnum
Definition: EnumDefinitions.h:535
Results.h
Cfsurfacelogvel::DeepEcho
void DeepEcho(void)
Definition: Cfsurfacelogvel.cpp:67
Cfsurfacelogvel::Echo
void Echo(void)
Definition: Cfsurfacelogvel.cpp:71
Gauss
Definition: Gauss.h:8
classes.h
DatasetInput2.h
: header file for datasetinput object
Cfsurfacelogvel::Cfsurfacelogvel
Cfsurfacelogvel()
Definition: Cfsurfacelogvel.cpp:27
Cfsurfacelogvel::Cfsurfacelogvel_Calculation
IssmDouble Cfsurfacelogvel_Calculation(Element *element, int definitionenum)
Definition: Cfsurfacelogvel.cpp:131
VxObsEnum
@ VxObsEnum
Definition: EnumDefinitions.h:848
Cfsurfacelogvel::lock
int lock
Definition: Cfsurfacelogvel.h:25
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16