Skip to content

Commit 19cd0e7

Browse files
mperegomperego
committed
Initial work on evolution of ice thickness.
Modified StokesFOThickness problem to enable thickness evolution using Tempus implicit schemes. Added a few fixes for Albany to run when MESH_DEPENDS_ON_SOLUTION is on. Added two tests computing thickness evolution and sensitivities (sensitivity test is enabled when MESH_DEPENDS_ON_SOLUTION is on).
1 parent 6ea155c commit 19cd0e7

26 files changed

+1046
-138
lines changed

‎CMakeLists.txt

+14-1
Original file line numberDiff line numberDiff line change
@@ -469,11 +469,24 @@ ELSE()
469469
SET(ALBANY_PYTHON FALSE)
470470
ENDIF()
471471

472+
# set optional dependency of mesh on solution, defaults to Disabled
473+
OPTION(ENABLE_MESH_DEPENDS_ON_SOLUTION "Flag to turn on dependency of mesh on parameters, e.g for shape optimization" OFF)
474+
IF (ENABLE_MESH_DEPENDS_ON_SOLUTION)
475+
MESSAGE("-- MESH_DEPENDS_ON_SOLUTION Enabled.")
476+
SET(ALBANY_MESH_DEPENDS_ON_SOLUTION TRUE)
477+
ELSE()
478+
MESSAGE("-- MESH_DEPENDS_ON_SOLUTION NOT Enabled.")
479+
SET(ALBANY_MESH_DEPENDS_ON_SOLUTION FALSE)
480+
ENDIF()
481+
472482
# set optional dependency of mesh on parameters, defaults to Disabled
473-
OPTION(ENABLE_MESH_DEPENDS_ON_PARAMETERS "Flag to turn on dependency of mesh on parameters, e.g for shape optimization" OFF)
483+
OPTION(ENABLE_MESH_DEPENDS_ON_PARAMETERS "Flag to turn on dependency of mesh on parameters, e.g for shape optimization" ${ALBANY_MESH_DEPENDS_ON_SOLUTION})
474484
IF (ENABLE_MESH_DEPENDS_ON_PARAMETERS)
475485
MESSAGE("-- MESH_DEPENDS_ON_PARAMETERS Enabled.")
476486
SET(ALBANY_MESH_DEPENDS_ON_PARAMETERS TRUE)
487+
ELSEIF(ALBANY_MESH_DEPENDS_ON_SOLUTION)
488+
MESSAGE(FATAL_ERROR
489+
"\nMESH_DEPENDS_ON_PARAMETERS cannot be OFF if MESH_DEPENDS_ON_SOLUTION is ON: ")
477490
ELSE()
478491
MESSAGE("-- MESH_DEPENDS_ON_PARAMETERS NOT Enabled.")
479492
SET(ALBANY_MESH_DEPENDS_ON_PARAMETERS FALSE)

‎src/Albany_config.h.in

+4-1
Original file line numberDiff line numberDiff line change
@@ -21,9 +21,12 @@
2121
#cmakedefine CMAKE_CXX_COMPILER_ID "${CMAKE_CXX_COMPILER_ID}"
2222
#cmakedefine CMAKE_CXX_COMPILER_VERSION "${CMAKE_CXX_COMPILER_VERSION}"
2323

24-
// Whether mesh and parameters depend on each others and/or solution
24+
// Whether mesh coordinates depend on parameters
2525
#cmakedefine ALBANY_MESH_DEPENDS_ON_PARAMETERS
2626

27+
// Whether mesh coordinates depend on solution
28+
#cmakedefine ALBANY_MESH_DEPENDS_ON_SOLUTION
29+
2730
// Whether the Omega_h interface is enabled
2831
#cmakedefine ALBANY_OMEGAH
2932

‎src/disc/stk/Albany_OrdinarySTKFieldContainer.cpp

+14-15
Original file line numberDiff line numberDiff line change
@@ -114,7 +114,7 @@ OrdinarySTKFieldContainer::OrdinarySTKFieldContainer(
114114
#endif
115115
//IKT FIXME? - currently won't write dxdp to output file if problem is steady,
116116
//as this output doesn't work in same way. May want to change in the future.
117-
const auto& sens_method = params_->get<std::string>("Sensitivity Method","NONE");
117+
const auto& sens_method = params_->get<std::string>("Sensitivity Method","None");
118118

119119
output_sens_field = this->num_params > 0 && num_time_deriv > 0 && sens_method != "None";
120120

@@ -179,24 +179,23 @@ OrdinarySTKFieldContainer::OrdinarySTKFieldContainer(
179179
}
180180

181181
//Transient sensitivities output to Exodus
182-
const int num_sens = (sens_method == "Forward") ? this->num_params : 1;
183-
for (int np = 0; np < num_sens; np++) {
184-
if (output_sens_field == true) {
185-
solution_field_dxdp[np] = &metaData_->declare_field<double>(
186-
stk::topology::NODE_RANK,
187-
params_->get<std::string>(
188-
sol_sens_tag_name_vec[np], sol_sens_id_name_vec[np]));
189-
stk::mesh::put_field_on_mesh(
190-
*solution_field_dxdp[np],
191-
metaData_->universal_part(),
192-
neq_,
193-
nullptr);
194-
}
182+
if(output_sens_field) {
183+
const int num_sens = (sens_method == "Forward") ? this->num_params : 1;
184+
for (int np = 0; np < num_sens; np++) {
185+
solution_field_dxdp[np] = &metaData_->declare_field<double>(
186+
stk::topology::NODE_RANK,
187+
params_->get<std::string>(
188+
sol_sens_tag_name_vec[np], sol_sens_id_name_vec[np]));
189+
stk::mesh::put_field_on_mesh(
190+
*solution_field_dxdp[np],
191+
metaData_->universal_part(),
192+
neq_,
193+
nullptr);
195194
#ifdef ALBANY_SEACAS
196-
if (output_sens_field == true)
197195
stk::io::set_field_role(
198196
*solution_field_dxdp[np], Ioss::Field::TRANSIENT);
199197
#endif
198+
}
200199
}
201200
}
202201

‎src/evaluators/interpolation/PHAL_DOFGradInterpolation.cpp

+3-1
Original file line numberDiff line numberDiff line change
@@ -10,4 +10,6 @@
1010
#include "PHAL_DOFGradInterpolation_Def.hpp"
1111

1212
PHAL_INSTANTIATE_TEMPLATE_CLASS_WITH_ONE_SCALAR_TYPE(PHAL::DOFGradInterpolationBase)
13-
PHAL_INSTANTIATE_TEMPLATE_CLASS_WITH_ONE_SCALAR_TYPE(PHAL::FastSolutionGradInterpolationBase)
13+
#if !defined(ALBANY_MESH_DEPENDS_ON_SOLUTION)
14+
PHAL_INSTANTIATE_TEMPLATE_CLASS_WITH_ONE_SCALAR_TYPE(PHAL::FastSolutionGradInterpolationBase)
15+
#endif

‎src/evaluators/interpolation/PHAL_DOFTensorGradInterpolation.cpp

+3-1
Original file line numberDiff line numberDiff line change
@@ -10,4 +10,6 @@
1010
#include "PHAL_DOFTensorGradInterpolation_Def.hpp"
1111

1212
PHAL_INSTANTIATE_TEMPLATE_CLASS_WITH_ONE_SCALAR_TYPE(PHAL::DOFTensorGradInterpolationBase)
13-
PHAL_INSTANTIATE_TEMPLATE_CLASS_WITH_ONE_SCALAR_TYPE(PHAL::FastSolutionTensorGradInterpolationBase)
13+
#if !defined(ALBANY_MESH_DEPENDS_ON_SOLUTION)
14+
PHAL_INSTANTIATE_TEMPLATE_CLASS_WITH_ONE_SCALAR_TYPE(PHAL::FastSolutionTensorGradInterpolationBase)
15+
#endif

‎src/evaluators/interpolation/PHAL_DOFVecGradInterpolation.cpp

+3
Original file line numberDiff line numberDiff line change
@@ -10,4 +10,7 @@
1010
#include "PHAL_DOFVecGradInterpolation_Def.hpp"
1111

1212
PHAL_INSTANTIATE_TEMPLATE_CLASS_WITH_ONE_SCALAR_TYPE(PHAL::DOFVecGradInterpolationBase)
13+
14+
#if !defined(ALBANY_MESH_DEPENDS_ON_SOLUTION)
1315
PHAL_INSTANTIATE_TEMPLATE_CLASS_WITH_ONE_SCALAR_TYPE(PHAL::FastSolutionVecGradInterpolationBase)
16+
#endif

‎src/landIce/evaluators/LandIce_Gather2DField.hpp

+2
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,7 @@ class Gather2DField : public PHX::EvaluatorWithBaseImpl<Traits>,
5252

5353
// Output:
5454
PHX::MDField<ScalarT,Cell,Node> field2D;
55+
PHX::MDField<ScalarT,Cell,Node> field2D_dot;
5556

5657
int offset; // Offset of first DOF being gathered
5758
int numNodes;
@@ -60,6 +61,7 @@ class Gather2DField : public PHX::EvaluatorWithBaseImpl<Traits>,
6061
std::string meshPart;
6162

6263
bool extruded;
64+
bool enableTransient;
6365
};
6466

6567
} // namespace LandIce

‎src/landIce/evaluators/LandIce_Gather2DField_Def.hpp

+174-1
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,13 @@ Gather2DField(const Teuchos::ParameterList& p,
3434
{
3535
this->addEvaluatedField(field2D);
3636

37+
const std::string transientFieldName = "Time Dependent 2D Field Name";
38+
enableTransient = p.isParameter(transientFieldName);
39+
if (enableTransient) {
40+
field2D_dot = PHX::MDField<ScalarT,Cell,Node>(p.get<std::string>(transientFieldName), dl->node_scalar);
41+
this->addEvaluatedField(field2D_dot);
42+
}
43+
3744
fieldLevel = p.get<int>("Field Level");
3845
TEUCHOS_TEST_FOR_EXCEPTION (fieldLevel<0, Teuchos::Exceptions::InvalidParameter,
3946
"[Gather2DField] Error! Field level must be non-negative.\n");
@@ -66,6 +73,8 @@ postRegistrationSetup(typename Traits::SetupData /* d */,
6673
PHX::FieldManager<Traits>& fm)
6774
{
6875
this->utils.setFieldData(field2D,fm);
76+
if (enableTransient)
77+
this->utils.setFieldData(field2D_dot,fm);
6978
}
7079

7180
template<typename EvalT, typename Traits>
@@ -91,6 +100,13 @@ evaluateFields (typename PHALTraits::EvalData workset)
91100
const auto& elem_lids = workset.disc->getElementLIDs_host(workset.wsIndex);
92101

93102
const auto x_data = Albany::getLocalData(workset.x);
103+
104+
Teuchos::ArrayRCP<const RealType> xdot_data;
105+
const bool gather_xdot = workset.transientTerms && enableTransient;
106+
if(gather_xdot)
107+
xdot_data = Albany::getLocalData(workset.xdot);
108+
109+
94110
const auto& dof_mgr = workset.disc->getDOFManager();
95111
const auto& elem_dof_lids = dof_mgr->elem_dof_lids().host();
96112
const auto& node_dof_mgr = workset.disc->getNodeDOFManager();
@@ -115,6 +131,8 @@ evaluateFields (typename PHALTraits::EvalData workset)
115131
for (int inode=0; inode<num_nodes; ++inode) {
116132
const LO ldof = elem_dof_lids(field_elem_LID,field_nodes[inode]);
117133
field2D(cell,nodes[inode]) = x_data[ldof];
134+
if(gather_xdot)
135+
field2D_dot(cell,nodes[inode]) = xdot_data[ldof];
118136
}
119137
};
120138

@@ -143,6 +161,8 @@ evaluateFields (typename PHALTraits::EvalData workset)
143161
for (int i=0; i<numSideNodes; ++i){
144162
const LO ldof = elem_dof_lids(elem_LID,offsets[i]);
145163
field2D(cell,nodes[i]) = x_data[ldof];
164+
if(gather_xdot)
165+
field2D_dot(cell,nodes[i]) = xdot_data[ldof];
146166
}
147167
}
148168
}
@@ -165,6 +185,11 @@ evaluateFields (typename PHALTraits::EvalData workset)
165185
" - num layers : " + std::to_string(layers_data->numLayers) + "\n");
166186

167187
const auto x_data = Albany::getLocalData(workset.x);
188+
Teuchos::ArrayRCP<const RealType> xdot_data;
189+
const bool gather_xdot = workset.transientTerms && enableTransient;
190+
if(gather_xdot)
191+
xdot_data = Albany::getLocalData(workset.xdot);
192+
168193
const auto& dof_mgr = workset.disc->getDOFManager();
169194
const auto& elem_dof_lids = dof_mgr->elem_dof_lids().host();
170195
const auto& node_dof_mgr = workset.disc->getNodeDOFManager();
@@ -193,6 +218,11 @@ evaluateFields (typename PHALTraits::EvalData workset)
193218
val = ScalarT(val.size(),x_data[ldof]);
194219
val.setUpdateValue(!workset.ignore_residual);
195220
val.fastAccessDx(firstunk) = workset.j_coeff;
221+
if(gather_xdot) {
222+
ref_t valdot = field2D_dot(cell,nodes[inode]);
223+
valdot = ScalarT(valdot.size(),xdot_data[ldof]);
224+
valdot.fastAccessDx(firstunk) = workset.m_coeff;
225+
}
196226
}
197227
};
198228

@@ -225,6 +255,137 @@ evaluateFields (typename PHALTraits::EvalData workset)
225255
ref_t val = field2D(cell,nodes[i]);
226256
val = ScalarT(val.size(), x_data[ldof]);
227257
val.fastAccessDx(nodes[i]*neq + offset) = workset.j_coeff;
258+
if(gather_xdot) {
259+
ref_t valdot = field2D_dot(cell,nodes[i]);
260+
valdot = ScalarT(valdot.size(),xdot_data[ldof]);
261+
valdot.fastAccessDx(nodes[i]*neq + offset) = workset.m_coeff;
262+
}
263+
}
264+
}
265+
}
266+
}
267+
268+
269+
//**********************************************************************
270+
template<>
271+
void Gather2DField<PHALTraits::Tangent, PHALTraits>::
272+
evaluateFields (typename PHALTraits::EvalData workset)
273+
{
274+
// Mesh data
275+
const auto& layers_data = workset.disc->getLayeredMeshNumberingLO();
276+
const auto& elem_lids = workset.disc->getElementLIDs_host(workset.wsIndex);
277+
const int bot = layers_data->bot_side_pos;
278+
const int top = layers_data->top_side_pos;
279+
280+
ALBANY_EXPECT (fieldLevel==0 || fieldLevel==layers_data->numLayers,
281+
"Field level must be 0 or match the number of layers in the mesh.\n"
282+
" - field level: " + std::to_string(fieldLevel) + "\n"
283+
" - num layers : " + std::to_string(layers_data->numLayers) + "\n");
284+
285+
const auto x_data = Albany::getLocalData(workset.x);
286+
Teuchos::ArrayRCP<const RealType> xdot_data;
287+
Teuchos::ArrayRCP<Teuchos::ArrayRCP<const ST>> Vx_data, Vxdot_data;
288+
const bool gather_Vx = Teuchos::nonnull(workset.Vx);
289+
if(gather_Vx)
290+
Vx_data = Albany::getLocalData(workset.Vx);
291+
const bool gather_xdot = Teuchos::nonnull(workset.xdot) && enableTransient;
292+
const bool gather_Vxdot = Teuchos::nonnull(workset.Vxdot) && enableTransient;
293+
if(gather_xdot)
294+
xdot_data = Albany::getLocalData(workset.xdot);
295+
if(gather_Vxdot)
296+
Vxdot_data = Albany::getLocalData(workset.Vxdot);
297+
298+
const auto& dof_mgr = workset.disc->getDOFManager();
299+
const auto& elem_dof_lids = dof_mgr->elem_dof_lids().host();
300+
const auto& node_dof_mgr = workset.disc->getNodeDOFManager();
301+
302+
const int neq = dof_mgr->getNumFields();
303+
304+
if (extruded) {
305+
#ifdef ALBANY_DEBUG
306+
check_topology(dof_mgr->get_topology());
307+
#endif
308+
const int field_layer = fieldLevel==0 ? 0 : fieldLevel-1;
309+
const int field_side_pos = field_layer==0 ? bot : top;
310+
const auto& field_nodes = dof_mgr->getGIDFieldOffsetsSide(offset,field_side_pos);
311+
for (std::size_t cell=0; cell<workset.numCells; ++cell ) {
312+
const int elem_LID = elem_lids(cell);
313+
const int basal_elem_LID = layers_data->getColumnId(elem_LID);
314+
const int field_elem_LID = layers_data->getId(basal_elem_LID,field_layer);
315+
316+
const auto f = [&] (const std::vector<int>& nodes) {
317+
const int num_nodes = nodes.size();
318+
for (int inode=0; inode<num_nodes; ++inode) {
319+
LO ldof = elem_dof_lids(field_elem_LID,field_nodes[inode]);
320+
321+
ref_t val = field2D(cell,nodes[inode]);
322+
if (gather_Vx) {
323+
val = TanFadType(val.size(),x_data[ldof]);
324+
for (int k=0; k<workset.num_cols_x; ++k)
325+
val.fastAccessDx(k) = workset.j_coeff*Vx_data[ldof][k];
326+
} else {
327+
val = TanFadType(x_data[ldof]);
328+
}
329+
330+
if(gather_xdot) {
331+
ref_t valdot = field2D_dot(cell,nodes[inode]);
332+
if (gather_Vxdot) {
333+
valdot = TanFadType(valdot.size(),xdot_data[ldof]);
334+
for (int k=0; k<workset.num_cols_x; ++k)
335+
valdot.fastAccessDx(k) = workset.m_coeff*Vxdot_data[ldof][k];
336+
} else {
337+
valdot = TanFadType(xdot_data[ldof]);
338+
}
339+
}
340+
}
341+
};
342+
343+
// Run lambda on both top and bottom nodes in the cell, but make sure
344+
// we order them in whatever way the nodes on $field_side_pos would be ordered
345+
f(node_dof_mgr->getGIDFieldOffsetsSide(0,bot,field_side_pos));
346+
f(node_dof_mgr->getGIDFieldOffsetsSide(0,top,field_side_pos));
347+
}
348+
} else {
349+
TEUCHOS_TEST_FOR_EXCEPTION (workset.sideSets.is_null(), std::logic_error,
350+
"Side sets defined in input file but not properly specified on the mesh.\n");
351+
352+
// Check for early return
353+
if (workset.sideSets->count(meshPart)==0) {
354+
return;
355+
}
356+
357+
for (const auto& side : workset.sideSets->at(meshPart)) {
358+
// Get the data that corresponds to the side
359+
const int cell = side.ws_elem_idx;
360+
const int elem_LID = elem_lids(cell);
361+
const int side_pos = side.side_pos;
362+
363+
// Note: explicitly as for the offsets to be ordered as the dofs on this side pos
364+
const auto& offsets = dof_mgr->getGIDFieldOffsetsSide(offset,side_pos);
365+
const auto& nodes = node_dof_mgr->getGIDFieldOffsetsSide(0, side_pos);
366+
const int numSideNodes = nodes.size();
367+
for (int i=0; i<numSideNodes; ++i){
368+
const LO ldof = elem_dof_lids(elem_LID,offsets[i]);
369+
ref_t val = field2D(cell,nodes[i]);
370+
371+
if (gather_Vx) {
372+
val = TanFadType(val.size(),x_data[ldof]);
373+
for (int k=0; k<workset.num_cols_x; ++k)
374+
val.fastAccessDx(k) = workset.j_coeff*Vx_data[ldof][k];
375+
} else {
376+
val = TanFadType(x_data[ldof]);
377+
}
378+
379+
if(gather_xdot) {
380+
ref_t valdot = field2D_dot(cell,nodes[i]);
381+
if (gather_Vxdot) {
382+
valdot = TanFadType(valdot.size(),xdot_data[ldof]);
383+
for (int k=0; k<workset.num_cols_x; ++k)
384+
valdot.fastAccessDx(k) = workset.m_coeff*Vxdot_data[ldof][k];
385+
} else {
386+
valdot = TanFadType(xdot_data[ldof]);
387+
}
388+
}
228389
}
229390
}
230391
}
@@ -279,6 +440,11 @@ evaluateFields (typename PHALTraits::EvalData workset)
279440
const auto& elem_lids = workset.disc->getElementLIDs_host(workset.wsIndex);
280441

281442
const auto x_data = Albany::getLocalData(workset.x);
443+
Teuchos::ArrayRCP<const RealType> xdot_data;
444+
const bool gather_xdot = workset.transientTerms && enableTransient;
445+
if(gather_xdot)
446+
xdot_data = Albany::getLocalData(workset.xdot);
447+
282448
const auto& dof_mgr = workset.disc->getDOFManager();
283449
const auto& elem_dof_lids = dof_mgr->elem_dof_lids().host();
284450
const auto& node_dof_mgr = workset.disc->getNodeDOFManager();
@@ -306,6 +472,10 @@ evaluateFields (typename PHALTraits::EvalData workset)
306472
const LO ldof = elem_dof_lids(field_elem_LID,field_nodes[inode]);
307473
ref_t val = field2D(cell,nodes[inode]);
308474
val = HessianVecFad(val.size(), x_data[ldof]);
475+
476+
if(gather_xdot)
477+
field2D_dot(cell,nodes[inode]) = xdot_data[ldof];
478+
309479
if (is_x_active) {
310480
int firstunk = neq*nodes[inode] + offset;
311481
val.fastAccessDx(firstunk).val() = workset.j_coeff;
@@ -348,7 +518,10 @@ evaluateFields (typename PHALTraits::EvalData workset)
348518
const LO ldof = elem_dof_lids(elem_LID,offsets[i]);
349519
ref_t val = field2D(cell,nodes[i]);
350520
val = HessianVecFad(val.size(), x_data[ldof]);
351-
521+
522+
if(gather_xdot)
523+
field2D_dot(cell,nodes[i]) = xdot_data[ldof];
524+
352525
// If we differentiate w.r.t. the solution, we have to set the first
353526
// derivative to workset.j_coeff
354527
if (is_x_active)

‎src/landIce/evaluators/LandIce_LaplacianRegularizationResidual.hpp

+2-2
Original file line numberDiff line numberDiff line change
@@ -55,8 +55,8 @@ namespace LandIce {
5555
PHX::MDField<const ScalarT> gradField;
5656
PHX::MDField<const ParamScalarT> forcing;
5757
PHX::MDField<const MeshScalarT> gradBF;
58-
PHX::MDField<const MeshScalarT> BF;
59-
PHX::MDField<const MeshScalarT> side_BF;
58+
PHX::MDField<const RealType> BF;
59+
PHX::MDField<const RealType> side_BF;
6060
PHX::MDField<const MeshScalarT> w_measure;
6161
PHX::MDField<const MeshScalarT> w_side_measure;
6262

0 commit comments

Comments
 (0)