Skip to content

Commit

Permalink
NEW: adding continuous design QMU capabilities and Mmemasstransport
Browse files Browse the repository at this point in the history
analysis. The Mmemasstransport analysis allows for dakota runs with mass
tranposrt thickness times series that are an ensemble -> very nice to
compute sea level change fingerprints in a MME manner.
  • Loading branch information
larour committed Nov 29, 2024
1 parent 84b0849 commit 2cbbce0
Show file tree
Hide file tree
Showing 30 changed files with 752 additions and 85 deletions.
14 changes: 14 additions & 0 deletions m4/analyses.m4
Original file line number Diff line number Diff line change
Expand Up @@ -457,6 +457,20 @@ fi
AM_CONDITIONAL([MASSTRANSPORT], [test x$HAVE_MASSTRANSPORT = xyes])
AC_MSG_RESULT($HAVE_MASSTRANSPORT)
dnl }}}
dnl with-Mmemasstransport{{{
AC_ARG_WITH([Mmemasstransport],
AS_HELP_STRING([--with-Mmemasstransport = YES], [compile with Mmemasstransport capabilities (default is yes)]),
[MMEMASSTRANSPORT=$withval],[MMEMASSTRANSPORT=yes])
AC_MSG_CHECKING(for Mmemasstransport capability compilation)
HAVE_MMEMASSTRANSPORT=no
if test "x$MMEMASSTRANSPORT" = "xyes"; then
HAVE_MMEMASSTRANSPORT=yes
AC_DEFINE([_HAVE_MMEMASSTRANSPORT_],[1],[with Mmemasstransport capability])
fi
AM_CONDITIONAL([MMEMASSTRANSPORT], [test x$HAVE_MMEMASSTRANSPORT = xyes])
AC_MSG_RESULT($HAVE_MMEMASSTRANSPORT)
dnl }}}
dnl with-Melting{{{
AC_ARG_WITH([Melting],
AS_HELP_STRING([--with-Melting = YES], [compile with Melting capabilities (default is yes)]),
Expand Down
6 changes: 5 additions & 1 deletion src/c/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -197,6 +197,7 @@ issm_sources += \
./toolkits/mpi/commops/GetOwnershipBoundariesFromRange.cpp \
./toolkits/ToolkitOptions.cpp \
./modules/MmeToInputFromIdx/MmeToInputFromIdx.cpp\
./modules/MmeToInputx/MmeToInputx.cpp\
./modules/ModelProcessorx/ModelProcessorx.cpp \
./modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp \
./modules/ModelProcessorx/NodesPartitioning.cpp \
Expand Down Expand Up @@ -248,6 +249,7 @@ issm_sources += \
./modules/Solverx/Solverx.cpp \
./modules/StochasticForcingx/StochasticForcingx.cpp \
./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp \
./modules/UpdateMmesx/UpdateMmesx.cpp \
./cores/ProcessArguments.cpp \
./cores/ResetBoundaryConditions.cpp \
./cores/WrapperCorePointerFromSolutionEnum.cpp \
Expand All @@ -274,6 +276,7 @@ issm_sources += \
./cores/transient_core.cpp \
./cores/steadystate_core.cpp \
./cores/masstransport_core.cpp \
./cores/mmemasstransport_core.cpp \
./cores/oceantransport_core.cpp \
./cores/depthaverage_core.cpp \
./cores/extrudefrombase_core.cpp \
Expand Down Expand Up @@ -497,7 +500,8 @@ if MELTING
issm_sources += ./analyses/MeltingAnalysis.cpp
endif
if MASSTRANSPORT
issm_sources += ./analyses/MasstransportAnalysis.cpp
issm_sources += ./analyses/MasstransportAnalysis.cpp \
./analyses/MmemasstransportAnalysis.cpp
endif
if OCEANTRANSPORT
issm_sources += ./analyses/OceantransportAnalysis.cpp
Expand Down
3 changes: 3 additions & 0 deletions src/c/analyses/EnumToAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,9 @@ Analysis* EnumToAnalysis(int analysis_enum){
#ifdef _HAVE_MASSTRANSPORT_
case MasstransportAnalysisEnum : return new MasstransportAnalysis();
#endif
#ifdef _HAVE_MMEMASSTRANSPORT_
case MmemasstransportAnalysisEnum : return new MmemasstransportAnalysis();
#endif
#ifdef _HAVE_MELTING_
case MeltingAnalysisEnum : return new MeltingAnalysis();
#endif
Expand Down
143 changes: 143 additions & 0 deletions src/c/analyses/MmemasstransportAnalysis.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
#include "./MmemasstransportAnalysis.h"
#include <math.h>
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../classes/Inputs/TransientInput.h"
#include "../classes/Inputs/TriaInput.h"
#include "../classes/gauss/Gauss.h"
#include "../shared/shared.h"
#include "../modules/modules.h"

/*Model processing*/
void MmemasstransportAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
/*No constraints*/
}/*}}}*/
void MmemasstransportAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
/*No loads*/
}/*}}}*/
void MmemasstransportAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
::CreateNodes(nodes,iomodel,MmemasstransportAnalysisEnum,P1Enum);
}/*}}}*/
int MmemasstransportAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
return 3;
}/*}}}*/
void MmemasstransportAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/

int nature=0;
bool isdakota=0;

/*Update elements: */
int counter=0;
for(int i=0;i<iomodel->numberofelements;i++){
if(iomodel->my_elements[i]){
Element* element=(Element*)elements->GetObjectByOffset(counter);
element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum);
counter++;
}
}

/*Plug inputs into element:*/
iomodel->FetchDataToInput(inputs,elements,"md.geometry.thickness",ThicknessEnum);
iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0);
iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MmemasstransportMaskIceLevelsetEnum);
iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MmemasstransportMaskOceanLevelsetEnum);
iomodel->FetchDataToInput(inputs,elements,"md.mmemasstransport.thickness", MmemasstransportThicknessEnum);

/*Initialize sea level cumulated sea level loads :*/
iomodel->ConstantToInput(inputs,elements,0.,AccumulatedDeltaIceThicknessEnum,P0Enum);
iomodel->ConstantToInput(inputs,elements,0.,OldAccumulatedDeltaIceThicknessEnum,P0Enum);
iomodel->ConstantToInput(inputs,elements,0.,DeltaIceThicknessEnum,P0Enum);

}/*}}}*/
void MmemasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/

int numoutputs;
char** requestedoutputs = NULL;

int nids,npart,nel;
IssmDouble* ids=NULL;
IssmDouble* partition = NULL;

iomodel->FetchData(&nel,"md.mesh.numberofelements");
iomodel->FetchData(&ids,&nids,NULL,"md.mmemasstransport.ids");
//_printf_("nids: " << nids << "\n"); for(int i=0;i<nids;i++)_printf_(ids[i] << "|"); _printf_("\n");
parameters->AddObject(new DoubleMatParam(MmemasstransportModelidsEnum,ids,nids,1));
iomodel->FetchData(&partition,&npart,NULL,"md.mmemasstransport.partition");
if (npart!=nel)_error_("MmemasstransportAnalysis:UpdateParameters: partition vector should be distributed over elements, not vertices!");
parameters->AddObject(new DoubleMatParam(MmemasstransportPartitionEnum,partition,nel,1));

xDelete<IssmDouble>(ids);
xDelete<IssmDouble>(partition);

/*Requested outputs*/
iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.mmemasstransport.requested_outputs");
parameters->AddObject(new IntParam(MmemasstransportNumRequestedOutputsEnum,numoutputs));
if(numoutputs)parameters->AddObject(new StringArrayParam(MmemasstransportRequestedOutputsEnum,requestedoutputs,numoutputs));
iomodel->DeleteData(&requestedoutputs,numoutputs,"md.mmemasstransport.requested_outputs");

}/*}}}*/

/*Finite Element Analysis*/
void MmemasstransportAnalysis::Core(FemModel* femmodel){/*{{{*/
_error_("not implemented");
}/*}}}*/
void MmemasstransportAnalysis::PreCore(FemModel* femmodel){/*{{{*/
_error_("not implemented");
}/*}}}*/
ElementVector* MmemasstransportAnalysis::CreateDVector(Element* element){/*{{{*/
/*Default, return NULL*/
return NULL;
}/*}}}*/
ElementMatrix* MmemasstransportAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
_error_("Not implemented");
}/*}}}*/
ElementMatrix* MmemasstransportAnalysis::CreateKMatrix(Element* element){/*{{{*/
_error_("not implemented yet");
}/*}}}*/
ElementVector* MmemasstransportAnalysis::CreatePVector(Element* element){/*{{{*/
_error_("not implemented yet");
}/*}}}*/
void MmemasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

/*do nothing:*/
return;
}/*}}}*/
void MmemasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
_error_("Not implemented yet");
}/*}}}*/
void MmemasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

/*Fetch number of nodes and dof for this finite element*/
IssmDouble time;
IssmDouble ice[3];
IssmDouble ocean[3];
IssmDouble height;
int numnodes = element->GetNumberOfNodes();

element->parameters->FindParam(&time,TimeEnum);

TriaInput* h_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportThicknessEnum,time)); _assert_(h_input);
TriaInput* i_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportMaskIceLevelsetEnum,time)); _assert_(i_input);
TriaInput* o_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportMaskOceanLevelsetEnum,time)); _assert_(o_input);

Gauss* gauss=element->NewGauss();
for(int iv=0;iv<3;iv++){
gauss->GaussVertex(iv);
i_input->GetInputValue(&ice[iv],gauss);
o_input->GetInputValue(&ocean[iv],gauss);
}
h_input->GetInputAverage(&height);
delete gauss;

/*Add thickness ice and ocean levelsets as inputs to the tria element: */
element->AddInput(ThicknessEnum,&height,P0Enum); //very important, do not change the type of ThicknessEnum to P1 when it should be P0!
element->AddInput(MaskIceLevelsetEnum,ice,P1Enum);
element->AddInput(MaskOceanLevelsetEnum,ocean,P1Enum);

}/*}}}*/
void MmemasstransportAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
/*Default, do nothing*/
return;
}/*}}}*/
34 changes: 34 additions & 0 deletions src/c/analyses/MmemasstransportAnalysis.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
/*! \file MmemasstransportAnalysis.h
* \brief: header file for generic external result object
*/

#ifndef _MmemasstransportAnalysis_
#define _MmemasstransportAnalysis_

/*Headers*/
#include "./Analysis.h"

class MmemasstransportAnalysis: public Analysis{

public:
/*Model processing*/
void CreateConstraints(Constraints* constraints,IoModel* iomodel);
void CreateLoads(Loads* loads, IoModel* iomodel);
void CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr=false);
int DofsPerNode(int** doflist,int domaintype,int approximation);
void UpdateElements(Elements* elements,Inputs* inputs,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);
void PreCore(FemModel* femmodel);
ElementVector* CreateDVector(Element* element);
ElementMatrix* CreateJacobianMatrix(Element* element);
ElementMatrix* CreateKMatrix(Element* element);
ElementVector* CreatePVector(Element* element);
void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index);
void InputUpdateFromSolution(IssmDouble* solution,Element* element);
void UpdateConstraints(FemModel* femmodel);
};
#endif
6 changes: 3 additions & 3 deletions src/c/analyses/OceantransportAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -49,9 +49,9 @@ void OceantransportAnalysis::UpdateElements(Elements* elements,Inputs* inputs,Io
iomodel->FetchData(&modelid,"md.dsl.modelid");

/*replace dataset of forcings with only one, the modelid'th:*/
MmeToInputFromIdx(inputs,elements,modelid,OceantransportSpcbottompressureEnum, P1Enum);
MmeToInputFromIdx(inputs,elements,modelid,OceantransportSpcdslEnum, P1Enum);
MmeToInputFromIdx(inputs,elements,modelid,OceantransportSpcstrEnum, P0Enum);
MmeToInputFromIdx(inputs,elements,NULL,modelid,OceantransportSpcbottompressureEnum, P1Enum);
MmeToInputFromIdx(inputs,elements,NULL,modelid,OceantransportSpcdslEnum, P1Enum);
MmeToInputFromIdx(inputs,elements,NULL,modelid,OceantransportSpcstrEnum, P0Enum);
}

/*Initialize sea level cumulated sea level loads :*/
Expand Down
21 changes: 20 additions & 1 deletion src/c/analyses/SealevelchangeAnalysis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ void SealevelchangeAnalysis::UpdateElements(Elements* elements,Inputs* inputs,Io
iomodel->ConstantToInput(inputs,elements,0.,BedNorthEnum,P1Enum);
iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum);

/*Initialize loads:*/
/*Initialize loads: no! this should be done by the corresponding mass transports!*/
iomodel->ConstantToInput(inputs,elements,0.,DeltaTwsEnum,P1Enum);
iomodel->ConstantToInput(inputs,elements,0.,DeltaIceThicknessEnum,P1Enum);
iomodel->ConstantToInput(inputs,elements,0.,DeltaBottomPressureEnum,P1Enum);
Expand Down Expand Up @@ -251,6 +251,9 @@ void SealevelchangeAnalysis::UpdateParameters(Parameters* parameters,IoModel* io
} /*}}}*/
}

/*Indicate we have not yet run the Geometry Core module: */
parameters->AddObject(new BoolParam(SealevelchangeGeometryDoneEnum,false));

parameters->FindParam(&grdmodel,GrdModelEnum);
parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
if(grdmodel==ElasticEnum && isgrd){
Expand Down Expand Up @@ -669,6 +672,22 @@ void SealevelchangeAnalysis::Core(FemModel* femmodel){/*{{{*/
}/*}}}*/
void SealevelchangeAnalysis::PreCore(FemModel* femmodel){/*{{{*/

int isuq=0;
int modelid=0;

/*Resolve Mmes using the modelid, if necessary: meaning if we are running a transient model and that UQ computations have not been triggered:*/
femmodel->parameters->FindParam(&isuq,QmuIsdakotaEnum);
if (!isuq && femmodel->inputs->GetInputObjectEnum(SolidearthExternalDisplacementEastRateEnum)==DatasetInputEnum){
femmodel->parameters->FindParam(&modelid,SolidearthExternalModelidEnum);

/*replace dataset of forcings with only one, the modelid'th:*/
MmeToInputFromIdx(femmodel->inputs,femmodel->elements,femmodel->parameters,modelid-1,SolidearthExternalDisplacementNorthRateEnum, P1Enum);
MmeToInputFromIdx(femmodel->inputs,femmodel->elements,femmodel->parameters,modelid-1,SolidearthExternalDisplacementEastRateEnum, P1Enum);
MmeToInputFromIdx(femmodel->inputs,femmodel->elements,femmodel->parameters,modelid-1,SolidearthExternalDisplacementUpRateEnum, P1Enum);
MmeToInputFromIdx(femmodel->inputs,femmodel->elements,femmodel->parameters,modelid-1,SolidearthExternalGeoidRateEnum, P1Enum);

}

/*run sea level change core geometry only once, after the Model Processor is done:*/
sealevelchange_initialgeometry(femmodel);

Expand Down
1 change: 1 addition & 0 deletions src/c/analyses/analyses.h
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@
#include "./HydrologyPismAnalysis.h"
#include "./LevelsetAnalysis.h"
#include "./MasstransportAnalysis.h"
#include "./MmemasstransportAnalysis.h"
#include "./OceantransportAnalysis.h"
#include "./SamplingAnalysis.h"
#include "./SmbAnalysis.h"
Expand Down
3 changes: 3 additions & 0 deletions src/c/cores/cores.h
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ void controlm1qn3_core(FemModel* femmodel);
void controladm1qn3_core(FemModel* femmodel);
void controlvalidation_core(FemModel* femmodel);
void masstransport_core(FemModel* femmodel);
void mmemasstransport_core(FemModel* femmodel);
void oceantransport_core(FemModel* femmodel);
void depthaverage_core(FemModel* femmodel);
void extrudefrombase_core(FemModel* femmodel);
Expand All @@ -45,6 +46,7 @@ void slopecompute_core(FemModel* femmodel);
void steadystate_core(FemModel* femmodel);
void transient_core(FemModel* femmodel);
void transient_precore(FemModel* femmodel);
void transient_postcore(FemModel* femmodel);
void dakota_core(FemModel* femmodel);
void ad_core(FemModel* femmodel);
void adgradient_core(FemModel* femmodel);
Expand All @@ -62,6 +64,7 @@ void debris_core(FemModel* femmodel);
#ifdef _HAVE_SEALEVELCHANGE_
void sealevelchange_core(FemModel* femmodel);
void sealevelchange_initialgeometry(FemModel* femmodel);
void sealevelchange_finalize(FemModel* femmodel);
SealevelGeometry* sealevelchange_geometry(FemModel* femmodel);
#endif
void grd_core(FemModel* femmodel,SealevelGeometry* slgeom);
Expand Down
Loading

0 comments on commit 2cbbce0

Please sign in to comment.