Skip to content

Commit

Permalink
NEW: retrieving change in gravitational acceleration from ESA solver
Browse files Browse the repository at this point in the history
  • Loading branch information
@AdhikariJPL
AdhikariJPL committed Oct 1, 2024
1 parent 40790d0 commit 7bb99ab
Show file tree
Hide file tree
Showing 16 changed files with 150 additions and 82 deletions.
31 changes: 25 additions & 6 deletions src/c/analyses/EsaAnalysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,14 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
int nl;
IssmDouble* love_h=NULL;
IssmDouble* love_l=NULL;
IssmDouble* love_k=NULL;

IssmDouble* U_elastic = NULL;
IssmDouble* U_elastic_local = NULL;
IssmDouble* H_elastic = NULL;
IssmDouble* H_elastic_local = NULL;
IssmDouble* G_elastic = NULL;
IssmDouble* G_elastic_local = NULL;
int M,m,lower_row,upper_row;
IssmDouble degacc=.01;
IssmDouble planetradius=0;
Expand Down Expand Up @@ -74,26 +77,30 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
/*love numbers: */
iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");

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

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

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

U_elastic_local[i-lower_row]= (love_h[nl-1])/2.0/sin(alpha/2.0);
H_elastic_local[i-lower_row]= 0;
U_elastic_local[i-lower_row]= 0.5*love_h[nl-1]/sin(alpha/2.0);
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)));
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.
//IssmDouble Pn,Pn1,Pn2;
//IssmDouble Pn_p,Pn_p1,Pn_p2;
IssmDouble Pn = 0.;
Expand All @@ -105,8 +112,12 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s

for (int n=0;n<nl;n++) {
IssmDouble deltalove_U;

deltalove_U = (love_h[n]-love_h[nl-1]);
IssmDouble deltalove_H;
IssmDouble deltalove_G;

deltalove_U = love_h[n]-love_h[nl-1];
deltalove_H = love_l[n] - (love_l[nl-1]*(nl-1)/(n+1e-12));
deltalove_G = 2*(love_h[n]-love_h[nl-1]) - (n*love_k[n]-love_k[nl-1]*(nl-1)) - love_k[n];

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

U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
H_elastic_local[i-lower_row] -= sin(alpha)*deltalove_H*Pn_p; // horizontal displacements
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.
//IssmDouble Pn,Pn1,Pn2;
}
}
/* }}} */
Expand All @@ -143,6 +156,7 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
/*All gather:*/
ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
/*Free resources: */
xDelete<int>(recvcounts);
xDelete<int>(displs);
Expand All @@ -154,14 +168,19 @@ void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
parameters->AddObject(new DoubleVecParam(EsaUElasticEnum,U_elastic,M));
H_elastic[0]=H_elastic[1];
parameters->AddObject(new DoubleVecParam(EsaHElasticEnum,H_elastic,M));
G_elastic[0]=G_elastic[1];
parameters->AddObject(new DoubleVecParam(EsaGElasticEnum,G_elastic,M));

/*Free resources: */
xDelete<IssmDouble>(love_h);
xDelete<IssmDouble>(love_l);
xDelete<IssmDouble>(love_k);
xDelete<IssmDouble>(U_elastic);
xDelete<IssmDouble>(U_elastic_local);
xDelete<IssmDouble>(H_elastic);
xDelete<IssmDouble>(H_elastic_local);
xDelete<IssmDouble>(G_elastic);
xDelete<IssmDouble>(G_elastic_local);

/*Transitions: */
iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.esa.transitions");
Expand Down
4 changes: 2 additions & 2 deletions src/c/classes/Elements/Element.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -406,8 +406,8 @@ class Element: public Object{
virtual void WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue){_error_("not implemented yet");};

#ifdef _HAVE_ESA_
virtual void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy)=0;
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;
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;
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;
#endif
#ifdef _HAVE_SEALEVELCHANGE_
virtual IssmDouble GetArea3D(void)=0;
Expand Down
4 changes: 2 additions & 2 deletions src/c/classes/Elements/Penta.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,8 +223,8 @@ class Penta: public Element,public ElementHook,public PentaRef{
#endif

#ifdef _HAVE_ESA_
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!");};
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!");};
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!");};
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!");};
#endif
#ifdef _HAVE_SEALEVELCHANGE_
void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
Expand Down
4 changes: 2 additions & 2 deletions src/c/classes/Elements/Seg.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,8 +174,8 @@ class Seg: public Element,public ElementHook,public SegRef{
IssmDouble GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]){_error_("not implemented yet");};

#ifdef _HAVE_ESA_
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!");};
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!");};
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!");};
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!");};
#endif
#ifdef _HAVE_SEALEVELCHANGE_
void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
Expand Down
4 changes: 2 additions & 2 deletions src/c/classes/Elements/Tetra.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,8 +181,8 @@ class Tetra: public Element,public ElementHook,public TetraRef{
void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);

#ifdef _HAVE_ESA_
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!");};
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!");};
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!");};
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!");};
#endif
#ifdef _HAVE_SEALEVELCHANGE_
void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
Expand Down
Loading

0 comments on commit 7bb99ab

Please sign in to comment.