Skip to content

Albany Knowledge

Jump to bottom
Cameron Smith edited this page Aug 13, 2021 · 4 revisions

SCOREC use in Albany as of August 2021

Evaluator System

Albany has a flexible system of “evaluators” which is how it builds different kinds of solvers. An evaluator is just a function that takes some fields as input and outputs other fields. Example evaluators are:

  1. GatherSolution : ‘gathers’ (copies) a solution field from the Tpetra solution vector with one entry per unique node to a PHX::MDField (per-workset array) with one entry per node use by an element (for each element, stores field values at its nodes)

  2. ScatterResidual: ‘scatters’ (copies) values from a residual field (PHX::MDField / per-workset array) to the Tpetra solution vector. See here for the residual specialization evaluator for this field: https://github.com/gahansen/Albany/blob/master/src/evaluators/PHAL_ScatterResidual_Def.hpp#L151

  3. LCM’s ConstitutiveModelInterface is the central “physics” evaluator if you will. it is a general wrapper around every single possible physical model available in LCM, and chooses one based on a given name as seen here: https://github.com/gahansen/Albany/blob/master/src/LCM/models/ConstitutiveModelInterface_Def.hpp#L291. The material model is set by the materials.xml input file.

  4. TransportResidual handles, among other things, the solution of the heat equation, i.e. sources, diffusion, etc.

  5. ThermoMechanicalCoefficients is responsible for (among others) calculating the derivative of temperature with respect to time. It does this with the simple dT = (T - T_old) / dt formula. T_old is a key problem: it is taken from the state array in the Workset structure: https://github.com/gahansen/Albany/blob/master/src/LCM/evaluators/ThermoMechanicalCoefficients_Def.hpp#L97 .

  6. SaveStateField copies values from a PHX::MDField into a state array in the Workset structure: https://github.com/gahansen/Albany/blob/master/src/evaluators/PHAL_SaveStateField_Def.hpp#L74.

All the evaluator code is templated on whether we are computing the residual or the jacobian (EvalT). This means the C++ compiler is generating two versions of the code. In some cases, automatic differentiation means there is no additional work to do for the jacobian case. For other evaluators, special cases are hand-coded for the residual and jacobian. Case in point, for the Jacobian, GatherSolution is responsible for initializing the automatic derivatives:

https://github.com/gahansen/Albany/blob/master/src/evaluators/PHAL_GatherSolution_Def.hpp#L594

The whole evaluator system can be seen as a graph, the output of one evaluator is plugged in as the input to another one. Adding this line:

<Parameter name="Phalanx Graph Visualization Detail" type="int" value="2"/>

To the “Problem” ParameterList in the Albany input file will output the evaluator graph(s) in GraphViz format.

There are at least three systems involved in handling state for Albany: The Workset structure https://github.com/gahansen/Albany/blob/master/src/PHAL_Workset.hpp#L50. Among other things, this keeps the solution and its derivatives as Tpetra_Vector objects. Some details on this: first, a very important variable, the workste.stateArrayPtr, is taken from the StateManager’s element state.

  1. The StateManager. One of the things this class does is that its updateStates function copies “foo” data to “foo_old” in the discretization’s own state arrays: https://github.com/gahansen/Albany/blob/master/src/Albany_StateManager.cpp#L712. In fact, it seems the StateManager is always acting on the state arrays provided by the discretization. Another important function: it initializes those state arrays with its setStateArrays call: https://github.com/gahansen/Albany/blob/master/src/Albany_StateManager.cpp#L404.
  2. The Phalanx evaluators and fields
  3. Last but not least, the APF Discretization itself. In addition, we have functions for copying data in and out of both nodal and QP fields. These function operate between the APF field storage and the discretization state arrays, which are the basis for the StateManager.

Here is an example of a state variable getting registered and saved in the context of the mechanics problem (this may not be too useful for restarts, but is good to know nonetheless) A parameterlist for the state variable is created by the state manager function ‘registerStateVariable’: https://github.com/gahansen/Albany/blob/master/src/LCM/problems/MechanicsProblem.hpp#L1183 https://github.com/gahansen/Albany/blob/master/src/Albany_StateManager.cpp#L33 An evaluator ‘SaveStateField’ is created, with this parameterlist as its input. https://github.com/gahansen/Albany/blob/master/src/LCM/problems/MechanicsProblem.hpp#L1206 https://github.com/gahansen/Albany/blob/master/src/evaluators/PHAL_SaveStateField_Def.hpp#L69

The storage in the evaluator system is mainly in the form of PHX::MDField objects. Those objects get created by the evaluator that evaluates them.

Computing Basis Functions

Basis functions are computed in the file

https://github.com/gahansen/Albany/blob/master/src/evaluators/PHAL_ComputeBasisFunctions_Def.hpp

Determinants of Jacobians (of the mapping from the parametric element space to physical space) are assumed to always be positive, as can be seen in:

https://github.com/trilinos/trilinos/blob/master/packages/intrepid/core/src/Discretization/FunctionSpaceTools/Intrepid_FunctionSpaceToolsDef.hpp#L1350 publicTrilinos/packages/intrepid/core/src/Discretization/FunctionSpaceTools/Intrepid_FunctionSpaceToolsDef.hpp line 1336

We should look into whether we are communicating tet vertex orderings in the order that Intrepid expects them (this is for triangles, I didn’t check tets)

The Intrepid setJacobian method seems to be consistent with APF's Jacobian formula, modulo a transpose maybe:

https://github.com/trilinos/trilinos/blob/master/packages/intrepid/core/src/Cell/Intrepid_CellToolsDef.hpp#L1197 https://github.com/SCOREC/core/blob/master/apf/apfVectorElement.cc#L51

Their linear triangle basis also seems consistent with APF's triangle basis:

https://github.com/trilinos/trilinos/blob/master/packages/intrepid/core/src/Discretization/Basis/Intrepid_HGRAD_TRI_C1_FEMDef.hpp#L134 https://github.com/SCOREC/core/blob/master/apf/apfShape.cc#L111

The two things above should just amount to getting the v0-v1 and v0-v2 vectors of a triangle. Here is their determinant code (which is only copy pasted like five times):

https://github.com/trilinos/trilinos/blob/master/packages/intrepid/core/src/Shared/Intrepid_RealSpaceToolsDef.hpp#L1577

The 3D case has some pivoting and uses a sub-matrix determinant, while the 2D case is straightforward but can be negative. Contrast this with APF's determinant math, which is set up such that 2D determinants are always positive (via sqrt(d^2) in the simple case):

https://github.com/SCOREC/core/blob/master/apf/apfVectorElement.cc#L84

Clone this wiki locally