Skip to content

Commit 2cbbce0

Browse files
larourlarour
committed
NEW: adding continuous design QMU capabilities and Mmemasstransport
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.
1 parent 84b0849 commit 2cbbce0

30 files changed

+752
-85
lines changed

‎m4/analyses.m4‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -457,6 +457,20 @@ fi
457457
AM_CONDITIONAL([MASSTRANSPORT], [test x$HAVE_MASSTRANSPORT = xyes])
458458
AC_MSG_RESULT($HAVE_MASSTRANSPORT)
459459
dnl }}}
460+
dnl with-Mmemasstransport{{{
461+
AC_ARG_WITH([Mmemasstransport],
462+
AS_HELP_STRING([--with-Mmemasstransport = YES], [compile with Mmemasstransport capabilities (default is yes)]),
463+
[MMEMASSTRANSPORT=$withval],[MMEMASSTRANSPORT=yes])
464+
AC_MSG_CHECKING(for Mmemasstransport capability compilation)
465+
466+
HAVE_MMEMASSTRANSPORT=no
467+
if test "x$MMEMASSTRANSPORT" = "xyes"; then
468+
HAVE_MMEMASSTRANSPORT=yes
469+
AC_DEFINE([_HAVE_MMEMASSTRANSPORT_],[1],[with Mmemasstransport capability])
470+
fi
471+
AM_CONDITIONAL([MMEMASSTRANSPORT], [test x$HAVE_MMEMASSTRANSPORT = xyes])
472+
AC_MSG_RESULT($HAVE_MMEMASSTRANSPORT)
473+
dnl }}}
460474
dnl with-Melting{{{
461475
AC_ARG_WITH([Melting],
462476
AS_HELP_STRING([--with-Melting = YES], [compile with Melting capabilities (default is yes)]),

‎src/c/Makefile.am‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -197,6 +197,7 @@ issm_sources += \
197197
./toolkits/mpi/commops/GetOwnershipBoundariesFromRange.cpp \
198198
./toolkits/ToolkitOptions.cpp \
199199
./modules/MmeToInputFromIdx/MmeToInputFromIdx.cpp\
200+
./modules/MmeToInputx/MmeToInputx.cpp\
200201
./modules/ModelProcessorx/ModelProcessorx.cpp \
201202
./modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp \
202203
./modules/ModelProcessorx/NodesPartitioning.cpp \
@@ -248,6 +249,7 @@ issm_sources += \
248249
./modules/Solverx/Solverx.cpp \
249250
./modules/StochasticForcingx/StochasticForcingx.cpp \
250251
./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp \
252+
./modules/UpdateMmesx/UpdateMmesx.cpp \
251253
./cores/ProcessArguments.cpp \
252254
./cores/ResetBoundaryConditions.cpp \
253255
./cores/WrapperCorePointerFromSolutionEnum.cpp \
@@ -274,6 +276,7 @@ issm_sources += \
274276
./cores/transient_core.cpp \
275277
./cores/steadystate_core.cpp \
276278
./cores/masstransport_core.cpp \
279+
./cores/mmemasstransport_core.cpp \
277280
./cores/oceantransport_core.cpp \
278281
./cores/depthaverage_core.cpp \
279282
./cores/extrudefrombase_core.cpp \
@@ -497,7 +500,8 @@ if MELTING
497500
issm_sources += ./analyses/MeltingAnalysis.cpp
498501
endif
499502
if MASSTRANSPORT
500-
issm_sources += ./analyses/MasstransportAnalysis.cpp
503+
issm_sources += ./analyses/MasstransportAnalysis.cpp \
504+
./analyses/MmemasstransportAnalysis.cpp
501505
endif
502506
if OCEANTRANSPORT
503507
issm_sources += ./analyses/OceantransportAnalysis.cpp

‎src/c/analyses/EnumToAnalysis.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -109,6 +109,9 @@ Analysis* EnumToAnalysis(int analysis_enum){
109109
#ifdef _HAVE_MASSTRANSPORT_
110110
case MasstransportAnalysisEnum : return new MasstransportAnalysis();
111111
#endif
112+
#ifdef _HAVE_MMEMASSTRANSPORT_
113+
case MmemasstransportAnalysisEnum : return new MmemasstransportAnalysis();
114+
#endif
112115
#ifdef _HAVE_MELTING_
113116
case MeltingAnalysisEnum : return new MeltingAnalysis();
114117
#endif

‎src/c/analyses/MmemasstransportAnalysis.cpp‎

Lines changed: 143 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
1+
#include "./MmemasstransportAnalysis.h"
2+
#include <math.h>
3+
#include "../toolkits/toolkits.h"
4+
#include "../classes/classes.h"
5+
#include "../classes/Inputs/TransientInput.h"
6+
#include "../classes/Inputs/TriaInput.h"
7+
#include "../classes/gauss/Gauss.h"
8+
#include "../shared/shared.h"
9+
#include "../modules/modules.h"
10+
11+
/*Model processing*/
12+
void MmemasstransportAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
13+
/*No constraints*/
14+
}/*}}}*/
15+
void MmemasstransportAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
16+
/*No loads*/
17+
}/*}}}*/
18+
void MmemasstransportAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
19+
::CreateNodes(nodes,iomodel,MmemasstransportAnalysisEnum,P1Enum);
20+
}/*}}}*/
21+
int MmemasstransportAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
22+
return 3;
23+
}/*}}}*/
24+
void MmemasstransportAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
25+
26+
int nature=0;
27+
bool isdakota=0;
28+
29+
/*Update elements: */
30+
int counter=0;
31+
for(int i=0;i<iomodel->numberofelements;i++){
32+
if(iomodel->my_elements[i]){
33+
Element* element=(Element*)elements->GetObjectByOffset(counter);
34+
element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum);
35+
counter++;
36+
}
37+
}
38+
39+
/*Plug inputs into element:*/
40+
iomodel->FetchDataToInput(inputs,elements,"md.geometry.thickness",ThicknessEnum);
41+
iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
42+
iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
43+
iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0);
44+
iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MmemasstransportMaskIceLevelsetEnum);
45+
iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MmemasstransportMaskOceanLevelsetEnum);
46+
iomodel->FetchDataToInput(inputs,elements,"md.mmemasstransport.thickness", MmemasstransportThicknessEnum);
47+
48+
/*Initialize sea level cumulated sea level loads :*/
49+
iomodel->ConstantToInput(inputs,elements,0.,AccumulatedDeltaIceThicknessEnum,P0Enum);
50+
iomodel->ConstantToInput(inputs,elements,0.,OldAccumulatedDeltaIceThicknessEnum,P0Enum);
51+
iomodel->ConstantToInput(inputs,elements,0.,DeltaIceThicknessEnum,P0Enum);
52+
53+
}/*}}}*/
54+
void MmemasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
55+
56+
int numoutputs;
57+
char** requestedoutputs = NULL;
58+
59+
int nids,npart,nel;
60+
IssmDouble* ids=NULL;
61+
IssmDouble* partition = NULL;
62+
63+
iomodel->FetchData(&nel,"md.mesh.numberofelements");
64+
iomodel->FetchData(&ids,&nids,NULL,"md.mmemasstransport.ids");
65+
//_printf_("nids: " << nids << "\n"); for(int i=0;i<nids;i++)_printf_(ids[i] << "|"); _printf_("\n");
66+
parameters->AddObject(new DoubleMatParam(MmemasstransportModelidsEnum,ids,nids,1));
67+
iomodel->FetchData(&partition,&npart,NULL,"md.mmemasstransport.partition");
68+
if (npart!=nel)_error_("MmemasstransportAnalysis:UpdateParameters: partition vector should be distributed over elements, not vertices!");
69+
parameters->AddObject(new DoubleMatParam(MmemasstransportPartitionEnum,partition,nel,1));
70+
71+
xDelete<IssmDouble>(ids);
72+
xDelete<IssmDouble>(partition);
73+
74+
/*Requested outputs*/
75+
iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.mmemasstransport.requested_outputs");
76+
parameters->AddObject(new IntParam(MmemasstransportNumRequestedOutputsEnum,numoutputs));
77+
if(numoutputs)parameters->AddObject(new StringArrayParam(MmemasstransportRequestedOutputsEnum,requestedoutputs,numoutputs));
78+
iomodel->DeleteData(&requestedoutputs,numoutputs,"md.mmemasstransport.requested_outputs");
79+
80+
}/*}}}*/
81+
82+
/*Finite Element Analysis*/
83+
void MmemasstransportAnalysis::Core(FemModel* femmodel){/*{{{*/
84+
_error_("not implemented");
85+
}/*}}}*/
86+
void MmemasstransportAnalysis::PreCore(FemModel* femmodel){/*{{{*/
87+
_error_("not implemented");
88+
}/*}}}*/
89+
ElementVector* MmemasstransportAnalysis::CreateDVector(Element* element){/*{{{*/
90+
/*Default, return NULL*/
91+
return NULL;
92+
}/*}}}*/
93+
ElementMatrix* MmemasstransportAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
94+
_error_("Not implemented");
95+
}/*}}}*/
96+
ElementMatrix* MmemasstransportAnalysis::CreateKMatrix(Element* element){/*{{{*/
97+
_error_("not implemented yet");
98+
}/*}}}*/
99+
ElementVector* MmemasstransportAnalysis::CreatePVector(Element* element){/*{{{*/
100+
_error_("not implemented yet");
101+
}/*}}}*/
102+
void MmemasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
103+
104+
/*do nothing:*/
105+
return;
106+
}/*}}}*/
107+
void MmemasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
108+
_error_("Not implemented yet");
109+
}/*}}}*/
110+
void MmemasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
111+
112+
/*Fetch number of nodes and dof for this finite element*/
113+
IssmDouble time;
114+
IssmDouble ice[3];
115+
IssmDouble ocean[3];
116+
IssmDouble height;
117+
int numnodes = element->GetNumberOfNodes();
118+
119+
element->parameters->FindParam(&time,TimeEnum);
120+
121+
TriaInput* h_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportThicknessEnum,time)); _assert_(h_input);
122+
TriaInput* i_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportMaskIceLevelsetEnum,time)); _assert_(i_input);
123+
TriaInput* o_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportMaskOceanLevelsetEnum,time)); _assert_(o_input);
124+
125+
Gauss* gauss=element->NewGauss();
126+
for(int iv=0;iv<3;iv++){
127+
gauss->GaussVertex(iv);
128+
i_input->GetInputValue(&ice[iv],gauss);
129+
o_input->GetInputValue(&ocean[iv],gauss);
130+
}
131+
h_input->GetInputAverage(&height);
132+
delete gauss;
133+
134+
/*Add thickness ice and ocean levelsets as inputs to the tria element: */
135+
element->AddInput(ThicknessEnum,&height,P0Enum); //very important, do not change the type of ThicknessEnum to P1 when it should be P0!
136+
element->AddInput(MaskIceLevelsetEnum,ice,P1Enum);
137+
element->AddInput(MaskOceanLevelsetEnum,ocean,P1Enum);
138+
139+
}/*}}}*/
140+
void MmemasstransportAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
141+
/*Default, do nothing*/
142+
return;
143+
}/*}}}*/

‎src/c/analyses/MmemasstransportAnalysis.h‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
/*! \file MmemasstransportAnalysis.h
2+
* \brief: header file for generic external result object
3+
*/
4+
5+
#ifndef _MmemasstransportAnalysis_
6+
#define _MmemasstransportAnalysis_
7+
8+
/*Headers*/
9+
#include "./Analysis.h"
10+
11+
class MmemasstransportAnalysis: public Analysis{
12+
13+
public:
14+
/*Model processing*/
15+
void CreateConstraints(Constraints* constraints,IoModel* iomodel);
16+
void CreateLoads(Loads* loads, IoModel* iomodel);
17+
void CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr=false);
18+
int DofsPerNode(int** doflist,int domaintype,int approximation);
19+
void UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type);
20+
void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
21+
22+
/*Finite element Analysis*/
23+
void Core(FemModel* femmodel);
24+
void PreCore(FemModel* femmodel);
25+
ElementVector* CreateDVector(Element* element);
26+
ElementMatrix* CreateJacobianMatrix(Element* element);
27+
ElementMatrix* CreateKMatrix(Element* element);
28+
ElementVector* CreatePVector(Element* element);
29+
void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
30+
void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index);
31+
void InputUpdateFromSolution(IssmDouble* solution,Element* element);
32+
void UpdateConstraints(FemModel* femmodel);
33+
};
34+
#endif

‎src/c/analyses/OceantransportAnalysis.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -49,9 +49,9 @@ void OceantransportAnalysis::UpdateElements(Elements* elements,Inputs* inputs,Io
4949
iomodel->FetchData(&modelid,"md.dsl.modelid");
5050

5151
/*replace dataset of forcings with only one, the modelid'th:*/
52-
MmeToInputFromIdx(inputs,elements,modelid,OceantransportSpcbottompressureEnum, P1Enum);
53-
MmeToInputFromIdx(inputs,elements,modelid,OceantransportSpcdslEnum, P1Enum);
54-
MmeToInputFromIdx(inputs,elements,modelid,OceantransportSpcstrEnum, P0Enum);
52+
MmeToInputFromIdx(inputs,elements,NULL,modelid,OceantransportSpcbottompressureEnum, P1Enum);
53+
MmeToInputFromIdx(inputs,elements,NULL,modelid,OceantransportSpcdslEnum, P1Enum);
54+
MmeToInputFromIdx(inputs,elements,NULL,modelid,OceantransportSpcstrEnum, P0Enum);
5555
}
5656

5757
/*Initialize sea level cumulated sea level loads :*/

‎src/c/analyses/SealevelchangeAnalysis.cpp‎

Lines changed: 20 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -79,7 +79,7 @@ void SealevelchangeAnalysis::UpdateElements(Elements* elements,Inputs* inputs,Io
7979
iomodel->ConstantToInput(inputs,elements,0.,BedNorthEnum,P1Enum);
8080
iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum);
8181

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

254+
/*Indicate we have not yet run the Geometry Core module: */
255+
parameters->AddObject(new BoolParam(SealevelchangeGeometryDoneEnum,false));
256+
254257
parameters->FindParam(&grdmodel,GrdModelEnum);
255258
parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
256259
if(grdmodel==ElasticEnum && isgrd){
@@ -669,6 +672,22 @@ void SealevelchangeAnalysis::Core(FemModel* femmodel){/*{{{*/
669672
}/*}}}*/
670673
void SealevelchangeAnalysis::PreCore(FemModel* femmodel){/*{{{*/
671674

675+
int isuq=0;
676+
int modelid=0;
677+
678+
/*Resolve Mmes using the modelid, if necessary: meaning if we are running a transient model and that UQ computations have not been triggered:*/
679+
femmodel->parameters->FindParam(&isuq,QmuIsdakotaEnum);
680+
if (!isuq && femmodel->inputs->GetInputObjectEnum(SolidearthExternalDisplacementEastRateEnum)==DatasetInputEnum){
681+
femmodel->parameters->FindParam(&modelid,SolidearthExternalModelidEnum);
682+
683+
/*replace dataset of forcings with only one, the modelid'th:*/
684+
MmeToInputFromIdx(femmodel->inputs,femmodel->elements,femmodel->parameters,modelid-1,SolidearthExternalDisplacementNorthRateEnum, P1Enum);
685+
MmeToInputFromIdx(femmodel->inputs,femmodel->elements,femmodel->parameters,modelid-1,SolidearthExternalDisplacementEastRateEnum, P1Enum);
686+
MmeToInputFromIdx(femmodel->inputs,femmodel->elements,femmodel->parameters,modelid-1,SolidearthExternalDisplacementUpRateEnum, P1Enum);
687+
MmeToInputFromIdx(femmodel->inputs,femmodel->elements,femmodel->parameters,modelid-1,SolidearthExternalGeoidRateEnum, P1Enum);
688+
689+
}
690+
672691
/*run sea level change core geometry only once, after the Model Processor is done:*/
673692
sealevelchange_initialgeometry(femmodel);
674693

‎src/c/analyses/analyses.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@
3737
#include "./HydrologyPismAnalysis.h"
3838
#include "./LevelsetAnalysis.h"
3939
#include "./MasstransportAnalysis.h"
40+
#include "./MmemasstransportAnalysis.h"
4041
#include "./OceantransportAnalysis.h"
4142
#include "./SamplingAnalysis.h"
4243
#include "./SmbAnalysis.h"

‎src/c/cores/cores.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@ void controlm1qn3_core(FemModel* femmodel);
3434
void controladm1qn3_core(FemModel* femmodel);
3535
void controlvalidation_core(FemModel* femmodel);
3636
void masstransport_core(FemModel* femmodel);
37+
void mmemasstransport_core(FemModel* femmodel);
3738
void oceantransport_core(FemModel* femmodel);
3839
void depthaverage_core(FemModel* femmodel);
3940
void extrudefrombase_core(FemModel* femmodel);
@@ -45,6 +46,7 @@ void slopecompute_core(FemModel* femmodel);
4546
void steadystate_core(FemModel* femmodel);
4647
void transient_core(FemModel* femmodel);
4748
void transient_precore(FemModel* femmodel);
49+
void transient_postcore(FemModel* femmodel);
4850
void dakota_core(FemModel* femmodel);
4951
void ad_core(FemModel* femmodel);
5052
void adgradient_core(FemModel* femmodel);
@@ -62,6 +64,7 @@ void debris_core(FemModel* femmodel);
6264
#ifdef _HAVE_SEALEVELCHANGE_
6365
void sealevelchange_core(FemModel* femmodel);
6466
void sealevelchange_initialgeometry(FemModel* femmodel);
67+
void sealevelchange_finalize(FemModel* femmodel);
6568
SealevelGeometry* sealevelchange_geometry(FemModel* femmodel);
6669
#endif
6770
void grd_core(FemModel* femmodel,SealevelGeometry* slgeom);

0 commit comments

Comments
 (0)