Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 21260) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 21261) @@ -466,4 +466,10 @@ endif #}}} +#Esa sources {{{ +if ESA +issm_sources += ./cores/esa_core.cpp\ + ./analyses/EsaAnalysis.cpp +endif +#}}} #Slr sources {{{ if SEALEVELRISE Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 21260) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 21261) @@ -107,4 +107,7 @@ case GiaAnalysisEnum : return new GiaAnalysis(); #endif + #ifdef _HAVE_ESA_ + case EsaAnalysisEnum : return new EsaAnalysis(); + #endif #ifdef _HAVE_MESHDEFORMATION_ case MeshdeformationAnalysisEnum : return new MeshdeformationAnalysis(); Index: /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 21261) +++ /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 21261) @@ -0,0 +1,199 @@ +#include "./EsaAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Model processing*/ +void EsaAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + /*No constraints*/ +}/*}}}*/ +void EsaAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + /*No loads*/ +}/*}}}*/ +void EsaAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + ::CreateNodes(nodes,iomodel,EsaAnalysisEnum,P1Enum); +}/*}}}*/ +int EsaAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ + return 1; +}/*}}}*/ +void EsaAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum); + counter++; + } + } + + /*Create inputs: */ + iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); + iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum); + iomodel->FetchDataToInput(elements,"md.esa.deltathickness",EsaDeltathicknessEnum); + +}/*}}}*/ +void EsaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + int nl; + IssmDouble* love_h=NULL; + IssmDouble* love_l=NULL; + + IssmDouble* U_elastic = NULL; + IssmDouble* U_elastic_local = NULL; + IssmDouble* H_elastic = NULL; + IssmDouble* H_elastic_local = NULL; + int M,m,lower_row,upper_row; + IssmDouble degacc=.01; + + int numoutputs; + char** requestedoutputs = NULL; + + /*transition vectors: */ + IssmDouble **transitions = NULL; + int *transitions_M = NULL; + int *transitions_N = NULL; + int ntransitions; + + /*love numbers: */ + iomodel->FetchData(&love_h,&nl,NULL,"md.esa.love_h"); + iomodel->FetchData(&love_l,&nl,NULL,"md.esa.love_l"); + + /*compute elastic green function for a range of angles*/ + iomodel->FetchData(°acc,"md.esa.degacc"); + M=reCast(180./degacc+1.); + U_elastic=xNew(M); + H_elastic=xNew(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(m); + H_elastic_local=xNew(m); + + /*compute U_elastic_local and H_elastic_local {{{ */ + for(int i=lower_row;i(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; + IssmDouble Pn,Pn1,Pn2; + IssmDouble Pn_p,Pn_p1,Pn_p2; + for (int n=0;n(IssmComm::GetSize()); + int* displs=xNew(IssmComm::GetSize()); + + //recvcounts: + ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); + + /*displs: */ + ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); + + /*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()); + /*free ressources: */ + xDelete(recvcounts); + xDelete(displs); + + /*}}}*/ + + /*Avoid singularity at 0: */ + U_elastic[0]=U_elastic[1]; + parameters->AddObject(new DoubleVecParam(EsaUElasticEnum,U_elastic,M)); + H_elastic[0]=H_elastic[1]; + parameters->AddObject(new DoubleVecParam(EsaHElasticEnum,H_elastic,M)); + + /*free ressources: */ + xDelete(love_h); + xDelete(love_l); + xDelete(U_elastic); + xDelete(U_elastic_local); + xDelete(H_elastic); + xDelete(H_elastic_local); + + /*Transitions: */ + iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.esa.transitions"); + if(transitions){ + parameters->AddObject(new DoubleMatArrayParam(EsaTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N)); + + for(int i=0;i(transition); + } + xDelete(transitions); + xDelete(transitions_M); + xDelete(transitions_N); + } + + /*Requested outputs*/ + iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.esa.requested_outputs"); + parameters->AddObject(new IntParam(EsaNumRequestedOutputsEnum,numoutputs)); + if(numoutputs)parameters->AddObject(new StringArrayParam(EsaRequestedOutputsEnum,requestedoutputs,numoutputs)); + iomodel->DeleteData(&requestedoutputs,numoutputs,"md.esa.requested_outputs"); + +}/*}}}*/ + +/*Finite Element Analysis*/ +void EsaAnalysis::Core(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementVector* EsaAnalysis::CreateDVector(Element* element){/*{{{*/ + /*Default, return NULL*/ + return NULL; +}/*}}}*/ +ElementMatrix* EsaAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ +_error_("Not implemented"); +}/*}}}*/ +ElementMatrix* EsaAnalysis::CreateKMatrix(Element* element){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +ElementVector* EsaAnalysis::CreatePVector(Element* element){/*{{{*/ +_error_("not implemented yet"); +}/*}}}*/ +void EsaAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void EsaAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ + _error_("Not implemented yet"); +}/*}}}*/ +void EsaAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + /*Default, do nothing*/ + return; + +}/*}}}*/ +void EsaAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + /*Default, do nothing*/ + return; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/EsaAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/EsaAnalysis.h (revision 21261) +++ /issm/trunk-jpl/src/c/analyses/EsaAnalysis.h (revision 21261) @@ -0,0 +1,33 @@ +/*! \file EsaAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _EsaAnalysis_ +#define _EsaAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class EsaAnalysis: public Analysis{ + + public: + /*Model processing*/ + void CreateConstraints(Constraints* constraints,IoModel* iomodel); + void CreateLoads(Loads* loads, IoModel* iomodel); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); + + /*Finite element Analysis*/ + void Core(FemModel* femmodel); + ElementVector* CreateDVector(Element* element); + ElementMatrix* CreateJacobianMatrix(Element* element); + ElementMatrix* CreateKMatrix(Element* element); + ElementVector* CreatePVector(Element* element); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); +}; +#endif Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 21260) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 21261) @@ -24,4 +24,5 @@ #include "./FreeSurfaceTopAnalysis.h" #include "./GiaAnalysis.h" +#include "./EsaAnalysis.h" #include "./HydrologyDCEfficientAnalysis.h" #include "./HydrologyDCInefficientAnalysis.h" Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21260) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21261) @@ -1527,4 +1527,8 @@ name==SealevelEustaticEnum || name==SealevelriseDeltathicknessEnum || + name==EsaUmotionEnum || + name==EsaNmotionEnum || + name==EsaEmotionEnum || + name==EsaDeltathicknessEnum || name==GiaWEnum || name==GiadWdtEnum || Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 21260) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 21261) @@ -296,4 +296,7 @@ virtual void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y)=0; #endif + #ifdef _HAVE_ESA_ + virtual void EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0; + #endif #ifdef _HAVE_SEALEVELRISE_ virtual IssmDouble GetArea3D(void)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 21260) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 21261) @@ -183,4 +183,7 @@ void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y); #endif + #ifdef _HAVE_ESA_ + void EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + #endif #ifdef _HAVE_SEALEVELRISE_ IssmDouble OceanArea(void){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 21260) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 21261) @@ -168,5 +168,7 @@ void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; #endif - +#ifdef _HAVE_ESA_ + void EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; +#endif #ifdef _HAVE_SEALEVELRISE_ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21260) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21261) @@ -175,4 +175,7 @@ void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; #endif +#ifdef _HAVE_ESA_ + void EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; +#endif #ifdef _HAVE_SEALEVELRISE_ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21260) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21261) @@ -3580,5 +3580,162 @@ /*}}}*/ #endif - +#ifdef _HAVE_ESA_ +void Tria::EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/ + + /*diverse:*/ + int gsize; + bool spherical=true; + IssmDouble llr_list[NUMVERTICES][3]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble area; + IssmDouble I; //ice/water loading + IssmDouble late,longe,re; + IssmDouble lati,longi,ri; + IssmDouble rho_ice,rho_water,rho_earth; + IssmDouble minlong=400; + IssmDouble maxlong=-20; + + /*precomputed elastic green functions:*/ + IssmDouble* U_elastic_precomputed = NULL; + IssmDouble* H_elastic_precomputed = NULL; + int M; + + /*computation of Green functions:*/ + IssmDouble* U_elastic= NULL; + IssmDouble* N_elastic= NULL; + IssmDouble* E_elastic= NULL; + + /*optimization:*/ + bool store_green_functions=false; + + /*early return if we are not on the (ice) loading point: */ + //if(!(this->inputs->Max(MaskIceLevelsetEnum)<0)) return; + + /*recover material parameters: */ + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + rho_earth=matpar->GetMaterialParameter(MaterialsEarthDensityEnum); + + /*how many dofs are we working with here? */ + this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /*compute area of element:*/ + area=GetAreaSpherical(); + + /*element centroid (spherical): */ + /* Where is the centroid of this element?:{{{*/ + ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); + + minlong=400; maxlong=-20; + for (int i=0;imaxlong)maxlong=llr_list[i][1]; + if(llr_list[i][1]180){ + if (llr_list[0][1]==0)llr_list[0][1]=360; + if (llr_list[1][1]==0)llr_list[1][1]=360; + if (llr_list[2][1]==0)llr_list[2][1]=360; + } + + // correction at the north pole + if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; + if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; + if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; + + //correction at the south pole + if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; + if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; + if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; + + late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0; + longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0; + + late=90-late; + if(longe>180)longe=(longe-180)-180; + + late=late/180*PI; + longe=longe/180*PI; + /*}}}*/ + + /*figure out gravity center of our element (Cartesian): */ + IssmDouble x_element, y_element, z_element; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; + y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; + z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; + + /*recover elastic Green's functions for displacement:*/ + DoubleVecParam* U_parameter = static_cast(this->parameters->FindParamObject(EsaUElasticEnum)); _assert_(U_parameter); + DoubleVecParam* H_parameter = static_cast(this->parameters->FindParamObject(EsaHElasticEnum)); _assert_(H_parameter); + U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M); + H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M); + + /*Compute ice thickness change: */ + Input* deltathickness_input=inputs->GetInput(EsaDeltathicknessEnum); + if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!"); + deltathickness_input->GetInputAverage(&I); + + /*initialize: */ + U_elastic=xNewZeroInit(gsize); + N_elastic=xNewZeroInit(gsize); + E_elastic=xNewZeroInit(gsize); + + int* indices=xNew(gsize); + IssmDouble* U_values=xNewZeroInit(gsize); + IssmDouble* N_values=xNewZeroInit(gsize); + IssmDouble* E_values=xNewZeroInit(gsize); + IssmDouble alpha; + IssmDouble delPhi,delLambda; + IssmDouble dx, dy, dz, x, y, z; + IssmDouble N_azim, E_azim; + + for(int i=0;i(alpha/PI*(M-1)); + U_elastic[i] += U_elastic_precomputed[index]; + N_elastic[i] += H_elastic_precomputed[index]*N_azim; + E_elastic[i] += H_elastic_precomputed[index]*E_azim; + + /*Add all components to the pUp solution vectors:*/ + U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i]; + N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i]; + E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i]; + } + pUp->SetValues(gsize,indices,U_values,ADD_VAL); + pNorth->SetValues(gsize,indices,N_values,ADD_VAL); + pEast->SetValues(gsize,indices,E_values,ADD_VAL); + + /*free ressources:*/ + xDelete(indices); + xDelete(U_values); xDelete(N_values); xDelete(E_values); + xDelete(U_elastic); xDelete(N_elastic); xDelete(E_elastic); + + return; +} +/*}}}*/ +#endif #ifdef _HAVE_SEALEVELRISE_ IssmDouble Tria::OceanArea(void){ /*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 21260) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 21261) @@ -144,4 +144,7 @@ void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y); #endif + #ifdef _HAVE_ESA_ + void EsaGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea); + #endif #ifdef _HAVE_SEALEVELRISE_ IssmDouble OceanArea(void); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21260) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21261) @@ -625,4 +625,8 @@ break; + case EsaSolutionEnum: + analyses_temp[numanalyses++]=EsaAnalysisEnum; + break; + case SealevelriseSolutionEnum: analyses_temp[numanalyses++]=SealevelriseAnalysisEnum; @@ -638,5 +642,5 @@ case TransientSolutionEnum:{ - bool isSIA,isFS,isthermal,isenthalpy,ismasstransport,isgroundingline,isstressbalance,ismovingfront,ishydrology,isdamage,issmb,isslr,isgia; + bool isSIA,isFS,isthermal,isenthalpy,ismasstransport,isgroundingline,isstressbalance,ismovingfront,ishydrology,isdamage,issmb,isslr,isesa,isgia; iomodel->FindConstant(&isSIA,"md.flowequation.isSIA"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); @@ -651,4 +655,5 @@ iomodel->FindConstant(&issmb,"md.transient.issmb"); iomodel->FindConstant(&isslr,"md.transient.isslr"); + iomodel->FindConstant(&isesa,"md.transient.isesa"); iomodel->FindConstant(&isgia,"md.transient.isgia"); if(isstressbalance){ @@ -694,4 +699,7 @@ if(isslr){ analyses_temp[numanalyses++]=SealevelriseAnalysisEnum; + } + if(isesa){ + analyses_temp[numanalyses++]=EsaAnalysisEnum; } if(isgia){ @@ -2291,4 +2299,55 @@ /*}}}*/ #endif +#ifdef _HAVE_ESA_ +void FemModel::EsaGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ + + IssmDouble eartharea=0; + IssmDouble eartharea_cpu=0; + + int ns,nsmax; + + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ + ns = elements->Size(); + + /*First, figure out the surface area of Earth: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + eartharea_cpu += element->GetAreaSpherical(); + } + ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Figure out max of ns: */ + ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Call the esa geodetic core: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + element->EsaGeodetic(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea); + } + if(i%100==0){ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + } + } + + /*One last time: */ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + + /*Free ressources:*/ + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + xDelete(xx); + xDelete(yy); + xDelete(zz); +} +/*}}}*/ +#endif #ifdef _HAVE_SEALEVELRISE_ void FemModel::SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21260) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21261) @@ -113,4 +113,7 @@ void Deflection(Vector* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y); #endif + #ifdef _HAVE_ESA_ + void EsaGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); + #endif #ifdef _HAVE_SEALEVELRISE_ void SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); Index: /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 21260) +++ /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 21261) @@ -62,4 +62,7 @@ solutioncore=&sealevelrise_core; break; + case EsaSolutionEnum: + solutioncore=&esa_core; + break; case GiaSolutionEnum: #if _HAVE_GIA_ Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 21260) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 21261) @@ -47,4 +47,5 @@ void dummy_core(FemModel* femmodel); void gia_core(FemModel* femmodel); +void esa_core(FemModel* femmodel); void smb_core(FemModel* femmodel); void damage_core(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/cores/esa_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/esa_core.cpp (revision 21261) +++ /issm/trunk-jpl/src/c/cores/esa_core.cpp (revision 21261) @@ -0,0 +1,104 @@ +/*!\file: esa_core.cpp + * \brief: core of the ESA solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void esa_core(FemModel* femmodel){ /*{{{*/ + + Vector *U_radial = NULL; + Vector *U_north = NULL; + Vector *U_east = NULL; + bool save_results,isesa,iscoupler; + int configuration_type; + int solution_type; + int numoutputs = 0; + char **requested_outputs = NULL; + + /*additional parameters: */ + int gsize; + bool spherical=true; + IssmDouble *latitude = NULL; + IssmDouble *longitude = NULL; + IssmDouble *radius = NULL; + IssmDouble *xx = NULL; + IssmDouble *yy = NULL; + IssmDouble *zz = NULL; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&isesa,TransientIsesaEnum); + femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); + + /*first, recover lat,long and radius vectors from vertices: */ + VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + + /*recover x,y,z vectors from vertices: */ + VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + + /*several cases here, depending on value of iscoupler and isesa: + solution_type == EsaSolutionEnum) we are running elastic adjustment core (no coupler) + ( !iscoupler & !isesa) we are not interested in being here :) + ( !iscoupler & isesa) we are running in uncoupled mode + ( iscoupler & isesa) we are running in coupled mode, and better be earth + ( iscoupler & !isesa) we are running in coupled mode, and better be an ice cap + */ + + if(solution_type==EsaSolutionEnum){ + isesa=1; + iscoupler=0; + } + + /*early return: */ + if( !iscoupler & !isesa) return; //we are not interested in being here :) + + /*In what follows we assume we are all running esa, either in coupled, or uncoupled mode:*/ + if(VerboseSolution()) _printf0_(" computing elastic adjustment\n"); + + /*set configuration: */ + if(isesa)femmodel->SetCurrentConfiguration(EsaAnalysisEnum); + + if(isesa){ + + /*compute components of 3-D crustal motion: */ + /*Initialize:*/ + U_radial = new Vector(gsize); + U_north = new Vector(gsize); + U_east = new Vector(gsize); + + /*call the geodetic main modlule:*/ + femmodel->EsaGeodetic(U_radial,U_north,U_east,latitude,longitude,radius,xx,yy,zz); + + /*get results into elements:*/ + InputUpdateFromVectorx(femmodel,U_radial,EsaUmotionEnum,VertexSIdEnum); // radial displacement + InputUpdateFromVectorx(femmodel,U_north,EsaNmotionEnum,VertexSIdEnum); // north motion + InputUpdateFromVectorx(femmodel,U_east,EsaEmotionEnum,VertexSIdEnum); // east motion + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + femmodel->parameters->FindParam(&requested_outputs,&numoutputs,EsaRequestedOutputsEnum); + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + } + + if(solution_type==EsaSolutionEnum)femmodel->RequestedDependentsx(); + + /*Free ressources:*/ + delete U_radial; + delete U_north; + delete U_east; + if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} + } + +} +/*}}}*/ + Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 21260) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 21261) @@ -77,4 +77,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.transient.isgroundingline",TransientIsgroundinglineEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.isgia",TransientIsgiaEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.transient.isesa",TransientIsesaEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.isdamageevolution",TransientIsdamageevolutionEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.ishydrology",TransientIshydrologyEnum)); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21260) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21261) @@ -343,4 +343,5 @@ TransientIsthermalEnum, TransientIsgiaEnum, + TransientIsesaEnum, TransientIsdamageevolutionEnum, TransientIshydrologyEnum, @@ -779,4 +780,15 @@ SealevelriseNumRequestedOutputsEnum, /*}}}*/ + /*Elastic adjustment (ESA){{{*/ + EsaUmotionEnum, + EsaNmotionEnum, + EsaEmotionEnum, + EsaDeltathicknessEnum, + EsaUElasticEnum, + EsaHElasticEnum, + EsaTransitionsEnum, + EsaRequestedOutputsEnum, + EsaNumRequestedOutputsEnum, + /*}}}*/ ParametersENDEnum, /*Coordinate Systems{{{*/ @@ -898,4 +910,6 @@ GiaSolutionEnum, GiaAnalysisEnum, + EsaSolutionEnum, + EsaAnalysisEnum, MeshdeformationAnalysisEnum, LevelsetAnalysisEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21260) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21261) @@ -349,4 +349,5 @@ case TransientIsthermalEnum : return "TransientIsthermal"; case TransientIsgiaEnum : return "TransientIsgia"; + case TransientIsesaEnum : return "TransientIsesa"; case TransientIsdamageevolutionEnum : return "TransientIsdamageevolution"; case TransientIshydrologyEnum : return "TransientIshydrology"; @@ -759,4 +760,13 @@ case SealevelriseRequestedOutputsEnum : return "SealevelriseRequestedOutputs"; case SealevelriseNumRequestedOutputsEnum : return "SealevelriseNumRequestedOutputs"; + case EsaUmotionEnum : return "EsaUmotion"; + case EsaNmotionEnum : return "EsaNmotion"; + case EsaEmotionEnum : return "EsaEmotion"; + case EsaDeltathicknessEnum : return "EsaDeltathickness"; + case EsaUElasticEnum : return "EsaUElastic"; + case EsaHElasticEnum : return "EsaHElastic"; + case EsaTransitionsEnum : return "EsaTransitions"; + case EsaRequestedOutputsEnum : return "EsaRequestedOutputs"; + case EsaNumRequestedOutputsEnum : return "EsaNumRequestedOutputs"; case ParametersENDEnum : return "ParametersEND"; case XYEnum : return "XY"; @@ -871,4 +881,6 @@ case GiaSolutionEnum : return "GiaSolution"; case GiaAnalysisEnum : return "GiaAnalysis"; + case EsaSolutionEnum : return "EsaSolution"; + case EsaAnalysisEnum : return "EsaAnalysis"; case MeshdeformationAnalysisEnum : return "MeshdeformationAnalysis"; case LevelsetAnalysisEnum : return "LevelsetAnalysis"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21260) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21261) @@ -355,4 +355,5 @@ else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum; else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; + else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum; else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum; else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; - else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; + if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; + else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum; else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; else if (strcmp(name,"Vx")==0) return VxEnum; - else if (strcmp(name,"VxPicard")==0) return VxPicardEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"VyAverage")==0) return VyAverageEnum; + if (strcmp(name,"VxPicard")==0) return VxPicardEnum; + else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; else if (strcmp(name,"Vy")==0) return VyEnum; else if (strcmp(name,"VyPicard")==0) return VyPicardEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; - else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; else stage=6; } if(stage==6){ - if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; + if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; + else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; @@ -751,9 +752,9 @@ else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum; - else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; + if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; + else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; @@ -777,4 +778,13 @@ else if (strcmp(name,"SealevelriseRequestedOutputs")==0) return SealevelriseRequestedOutputsEnum; else if (strcmp(name,"SealevelriseNumRequestedOutputs")==0) return SealevelriseNumRequestedOutputsEnum; + else if (strcmp(name,"EsaUmotion")==0) return EsaUmotionEnum; + else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum; + else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum; + else if (strcmp(name,"EsaDeltathickness")==0) return EsaDeltathicknessEnum; + else if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum; + else if (strcmp(name,"EsaHElastic")==0) return EsaHElasticEnum; + else if (strcmp(name,"EsaTransitions")==0) return EsaTransitionsEnum; + else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum; + else if (strcmp(name,"EsaNumRequestedOutputs")==0) return EsaNumRequestedOutputsEnum; else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; else if (strcmp(name,"XY")==0) return XYEnum; @@ -865,5 +875,8 @@ else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum; else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; - else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum; else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum; else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; @@ -875,8 +888,5 @@ else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; + else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; @@ -892,4 +902,6 @@ else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum; else if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum; + else if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum; + else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum; else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum; else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;