Skip to content

Commit c24a91f

Browse files
MathieuMorlighemMathieuMorlighem
committed
Merge branch 'main' of github.com:ISSMteam/ISSM
2 parents 58c5d9f + 1744867 commit c24a91f

File tree

20 files changed

+338
-101
lines changed

20 files changed

+338
-101
lines changed

‎src/c/analyses/AgeAnalysis.cpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -245,7 +245,7 @@ ElementVector* AgeAnalysis::CreatePVector(Element* element){/*{{{*/
245245
/*Intermediaries*/
246246
int stabilization;
247247
IssmDouble Jdet,dt,factor;
248-
IssmDouble temperature;
248+
IssmDouble age;
249249
IssmDouble tau_parameter,diameter,hx,hy,hz;
250250
IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
251251
IssmDouble u,v,w;
@@ -267,8 +267,8 @@ ElementVector* AgeAnalysis::CreatePVector(Element* element){/*{{{*/
267267
Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
268268
Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
269269
Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
270-
Input* temperature_input = NULL;
271-
if(reCast<bool,IssmDouble>(dt)){temperature_input = element->GetInput(TemperatureEnum); _assert_(temperature_input);}
270+
Input* age_input = NULL;
271+
if(reCast<bool,IssmDouble>(dt)){age_input = element->GetInput(TemperatureEnum); _assert_(age_input);}
272272

273273
/* Start looping on the number of gaussian points: */
274274
Gauss* gauss=element->NewGauss(4);
@@ -277,15 +277,15 @@ ElementVector* AgeAnalysis::CreatePVector(Element* element){/*{{{*/
277277
element->JacobianDeterminant(&Jdet,xyz_list,gauss);
278278
element->NodalFunctions(basis,gauss);
279279

280+
/*The RHS is 1*/
280281
scalar_def=1.*Jdet*gauss->weight;
281282
if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
282-
283283
for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
284284

285285
/* Build transient now */
286286
if(reCast<bool,IssmDouble>(dt)){
287-
temperature_input->GetInputValue(&temperature, gauss);
288-
scalar_transient=temperature*Jdet*gauss->weight;
287+
age_input->GetInputValue(&age, gauss);
288+
scalar_transient=age*Jdet*gauss->weight;
289289
for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
290290
}
291291

‎src/c/analyses/EsaAnalysis.cpp‎

Lines changed: 25 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -39,11 +39,14 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
3939
int nl;
4040
IssmDouble* love_h=NULL;
4141
IssmDouble* love_l=NULL;
42+
IssmDouble* love_k=NULL;
4243

4344
IssmDouble* U_elastic = NULL;
4445
IssmDouble* U_elastic_local = NULL;
4546
IssmDouble* H_elastic = NULL;
4647
IssmDouble* H_elastic_local = NULL;
48+
IssmDouble* G_elastic = NULL;
49+
IssmDouble* G_elastic_local = NULL;
4750
int M,m,lower_row,upper_row;
4851
IssmDouble degacc=.01;
4952
IssmDouble planetradius=0;
@@ -74,26 +77,30 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
7477
/*love numbers: */
7578
iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
7679
iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
80+
iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
7781

7882
/*compute elastic green function for a range of angles*/
7983
iomodel->FetchData(&degacc,"md.esa.degacc");
8084
M=reCast<int,IssmDouble>(180./degacc+1.);
8185
U_elastic=xNew<IssmDouble>(M);
8286
H_elastic=xNew<IssmDouble>(M);
87+
G_elastic=xNew<IssmDouble>(M);
8388

8489
/*compute combined legendre + love number (elastic green function:*/
8590
m=DetermineLocalSize(M,IssmComm::GetComm());
8691
GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
8792
U_elastic_local=xNew<IssmDouble>(m);
8893
H_elastic_local=xNew<IssmDouble>(m);
94+
G_elastic_local=xNew<IssmDouble>(m);
8995

90-
/*compute U_elastic_local and H_elastic_local {{{ */
96+
/*compute U_elastic_local, H_elastic_local, and G_elastic_local {{{ */
9197
for(int i=lower_row;i<upper_row;i++){
9298
IssmDouble alpha,x;
9399
alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
94100

95-
U_elastic_local[i-lower_row]= (love_h[nl-1])/2.0/sin(alpha/2.0);
96-
H_elastic_local[i-lower_row]= 0;
101+
U_elastic_local[i-lower_row]= 0.5*love_h[nl-1]/sin(alpha/2.0);
102+
H_elastic_local[i-lower_row]= -love_l[nl-1]*(nl-1) * cos(alpha/2)*(1 + 2*sin(alpha/2)) / (2*sin(alpha/2)*(1 + sin(alpha/2)));
103+
G_elastic_local[i-lower_row]= -(-0.25 + love_h[nl-1] - 0.5*love_k[nl-1]*(nl-1)) / sin(alpha/2.0); // negative sign is imposed to mean g_deformed_earth minus g_initial_undeformed_earth. Ferrell defined it as the difference in g between the undeformed initial Earth and deformed Earth.
97104
//IssmDouble Pn,Pn1,Pn2;
98105
//IssmDouble Pn_p,Pn_p1,Pn_p2;
99106
IssmDouble Pn = 0.;
@@ -105,8 +112,12 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
105112

106113
for (int n=0;n<nl;n++) {
107114
IssmDouble deltalove_U;
108-
109-
deltalove_U = (love_h[n]-love_h[nl-1]);
115+
IssmDouble deltalove_H;
116+
IssmDouble deltalove_G;
117+
118+
deltalove_U = love_h[n]-love_h[nl-1];
119+
deltalove_H = love_l[n] - (love_l[nl-1]*(nl-1)/(n+1e-12));
120+
deltalove_G = 2*(love_h[n]-love_h[nl-1]) - (n*love_k[n]-love_k[nl-1]*(nl-1)) - love_k[n];
110121

111122
/*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
112123
if(n==0){
@@ -125,7 +136,9 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
125136
Pn_p2=Pn_p1; Pn_p1=Pn_p;
126137

127138
U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
128-
H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
139+
H_elastic_local[i-lower_row] -= sin(alpha)*deltalove_H*Pn_p; // horizontal displacements
140+
G_elastic_local[i-lower_row] -= deltalove_G*Pn; // change in gravitational acceleration => negative sign is imposed to mean g_deformed_earth minus g_initial_undeformed_earth. Ferrell defined it as the difference in g between the undeformed initial Earth and deformed Earth.
141+
//IssmDouble Pn,Pn1,Pn2;
129142
}
130143
}
131144
/* }}} */
@@ -143,6 +156,7 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
143156
/*All gather:*/
144157
ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
145158
ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
159+
ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
146160
/*Free resources: */
147161
xDelete<int>(recvcounts);
148162
xDelete<int>(displs);
@@ -154,14 +168,19 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
154168
parameters->AddObject(new DoubleVecParam(EsaUElasticEnum,U_elastic,M));
155169
H_elastic[0]=H_elastic[1];
156170
parameters->AddObject(new DoubleVecParam(EsaHElasticEnum,H_elastic,M));
171+
G_elastic[0]=G_elastic[1];
172+
parameters->AddObject(new DoubleVecParam(EsaGElasticEnum,G_elastic,M));
157173

158174
/*Free resources: */
159175
xDelete<IssmDouble>(love_h);
160176
xDelete<IssmDouble>(love_l);
177+
xDelete<IssmDouble>(love_k);
161178
xDelete<IssmDouble>(U_elastic);
162179
xDelete<IssmDouble>(U_elastic_local);
163180
xDelete<IssmDouble>(H_elastic);
164181
xDelete<IssmDouble>(H_elastic_local);
182+
xDelete<IssmDouble>(G_elastic);
183+
xDelete<IssmDouble>(G_elastic_local);
165184

166185
/*Transitions: */
167186
iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.esa.transitions");

‎src/c/classes/Elements/Element.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -406,8 +406,8 @@ class Element: public Object{
406406
virtual void WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue){_error_("not implemented yet");};
407407

408408
#ifdef _HAVE_ESA_
409-
virtual void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy)=0;
410-
virtual void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz)=0;
409+
virtual void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast, Vector<IssmDouble>* pGravity, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy)=0;
410+
virtual void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pGravity, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz)=0;
411411
#endif
412412
#ifdef _HAVE_SEALEVELCHANGE_
413413
virtual IssmDouble GetArea3D(void)=0;

‎src/c/classes/Elements/Penta.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -223,8 +223,8 @@ class Penta: public Element,public ElementHook,public PentaRef{
223223
#endif
224224

225225
#ifdef _HAVE_ESA_
226-
void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
227-
void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
226+
void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pGravity,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
227+
void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pGravity,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
228228
#endif
229229
#ifdef _HAVE_SEALEVELCHANGE_
230230
void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};

‎src/c/classes/Elements/Seg.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -174,8 +174,8 @@ class Seg: public Element,public ElementHook,public SegRef{
174174
IssmDouble GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]){_error_("not implemented yet");};
175175

176176
#ifdef _HAVE_ESA_
177-
void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
178-
void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
177+
void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pGravity,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
178+
void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pGravity,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
179179
#endif
180180
#ifdef _HAVE_SEALEVELCHANGE_
181181
void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};

‎src/c/classes/Elements/Tetra.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -181,8 +181,8 @@ class Tetra: public Element,public ElementHook,public TetraRef{
181181
void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
182182

183183
#ifdef _HAVE_ESA_
184-
void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
185-
void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
184+
void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pGravity,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
185+
void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pGravity,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
186186
#endif
187187
#ifdef _HAVE_SEALEVELCHANGE_
188188
void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};

0 commit comments

Comments
 (0)