Skip to content

Tempus: Infrastructure for the initial conditions.

James Willenbring requested to merge tempus-add-stepper-initial-conditions-3715 into develop

Created by: ccober6

@trilinos/tempus

Description

In IntegratorBasic::initializeSolutionHistory() and IntegratorBasic::setSolutionHistory(), setInitialConditions() calls the Stepper to get the first SolutionState set to the initial conditions and make sure the initial conditions are "consistent", i.e., x, xDot and xDotDot satisfy the governing equations.

  • Change name of IntegratorBasic::setInitialState() to IntegratorBasic::initializeSolutionHistory() and also in all derived classes and tests.

  • Added consistency parameter "Initial Condition Consistency" which indicates which type of consistency should be applied to the initial conditions (ICs): 'None' - Do nothing to the ICs provided in the SolutionHistory. 'Zero' - Set the derivative of the SolutionState to zero in the SolutionHistory provided, e.g., xDot^0 = 0, or xDotDot^0 = 0. 'App' - Use the application's ICs, e.g., getNominalValues(). 'Consistent' - Make the initial conditions for x and xDot consistent with the governing equations, e.g., xDot = f(x,t), and f(x, xDot, t) = 0. For implicit ODEs, this requires a solve of f(x, xDot, t) = 0 for xDot, and another Jacobian and residual may be needed, e.g., boundary conditions on xDot may need to replace boundary conditions on x. For explicit ODEs, the default is 'Consistent', because it is one additional evaluation. For implicit ODEs, the default is 'None', because the application very often knows its IC and is easier and cheaper than the solve for 'Consistent'.

  • Added consistency check parameter "Initial Condition Consistency Check" which will check if the initial condition, x and xDot, is consistent with the governing equations, xDot = f(x,t), or f(x, xDot, t) = 0. The default value is true.

  • Started the migration to new "ModelEvaluator" interface by changing the solve and evaluate interface. For implicit steppers, the solve call is solveImplicitODE(x, xDot, time, p), which closely matches the mathematics, f(x, xDot, time, p). "p" are the parameters needed by the models, and currently includes :

    • TimeDerivative object which defines the time derivative for the stepper
    • Time step size
    • Stage number for Runge-Kutta methods
    • Alpha, the derviative of d(xDot)/d(x)
    • Beta, the derviative of d(x)/d(x) (often = 1)
    • EvaluationType, an enum to indicate how to evaluate the model . EVALUATE_RESIDUAL - Evaluate residual for the implicit ODE . SOLVE_FOR_X - Solve for x and determine xDot from x. . SOLVE_FOR_XDOT_CONST_X - Solve for xDot keeping x constant (for ICs). Note that this makes Alpha = d(xDot)/d(xDot) = 1 and Beta = d(x)/d(xDot) = 0. The implicit interface also has a residual evaluation call, evaluateImplicitODE(f, x, xDot, time, p), to evaluate f(x, xDot, t, p) = r for the initial condition consistency check. For explicit steppers, the evaluate call is evaluateExplicitODE(xDot, x, time), which closely matches the mathematics, xDot = f(x, t).

  • Added new methods, getUseFSAL() and setUseFSAL(), to indicate when the Stepper can use the First-Step-As-Last (FSAL) principle, where the last function evaluation, f(x^{n-1},t^{n-1}) [a.k.a. xDot^{n-1}], can be used for the first function evaluation, f(x^n,t^n). Often the Stepper can reuse xDot^{n-1}, but an exception is with operator splitting when xDot^{n-1} is invalid due to other Steppers changing the solution.

  • With "Use FSAL"=true, the time levels have consistent x, xDot and xDotDot, so isSynced=true for the state. For "Use FSAL"=false, x, xDot and xDotDot are not necessarily at the same time level, so isSynced=false.

  • Added new methods to get/set x, xDot, and/or xDotDot to the Stepper, (getStepperX(), getStepperXDot() and getStepperXDotDot()) which will return an RCP to either SolutionState memory or Stepper allocated memory. This reduces the memory needs for the solution in SolutionState while provides the Stepper with minimum required memory.

Miscellaneous Mods:

  • Create StepperExplicit base class. Base class many of the functions from explicit steppers.
  • Added a default constructor for SolutionState, which does not set the solution vectors, x, xdot and xdotdot, which should be set via setX(), setXDot(), and/or setXDotDot() prior to being added to SolutionHistory.
  • Removed unused constructors for SolutionState.
  • Removed experimental and unused SolutionState::swapSolutonData().
  • Added set methods for x, xDot and xDotDot in SolutionState, which can be used with the default constructor.
  • SolutionState::copySolutionData() now clones x, xDot and xDotDot when those vectors are null.
  • Added new method to set the solution status based on the solver return status, e.g., workingState->setSolutionStatus(sStatus).

Motivation and Context

This change adds the ability to set and make consistent the initial conditions (x, xDot and xDotDot) and use the FSAL principle, which can reduce the computational costs.

Related Issues

How Has This Been Tested?

Passed all tests.

  • Explicit Steppers

    • General defaults unless noted . "Initial Condition Consistency" = "Consistent" . "Initial Condition Consistency Check" = true
    • ForwardEuler: . Default ("Use FSAL"=true) reproduces previous solutions, and "Use FSAL"=false works and is tested with StepperOperatorSplitting.
    • ExplicitRK: . Default ("Use FSAL"=false) reproduces previous solutions, and "Use FSAL"=true works and is tested with Bogacki-Shampine 3(2) method.
    • Leapfrog: . Default ("Use FSAL"=false) reproduces previous solutions, and "Use FSAL"=true will also work (i.e., no-op), but issue a warning that it will have no affect.
    • NewmarkExplicitAForm: . Default ("Use FSAL"=true) reproduces previous solutions, and "Use FSAL"=false works and is tested with SinCos, BallParabolic and HarmonicOscillar problems.

  • Implicit Stepper

    • General defaults unless noted . "Use FSAL"=false . "Initial Condition Consistency" = "None" . "Initial Condition Consistency Check" = false
    • BDF2:
    • BackwardEuler: . "Initial Condition Consistency Check" = true for CDR and van der Pol test.
    • DIRK:
    • HHTAlpha: . Not fully updated because it does not test beyond Newmark. Likely will be replaced by Generalize-Alpha method stepper.
    • IMEX_RK_Partition:
    • IMEX_RK:
    • NewmarkImplicitAForm: . Default ("Use FSAL"=true) reproduces previous solutions.
    • NewmarkImplicitDForm: . Not fully updated because it does not have any tests. Likely will be replaced by Generalize-Alpha method stepper.
    • Trapezoidal: . Required ("Use FSAL"=true) reproduces previous solutions.

  • OperatorSplit Stepper

    • In most cases, subSteppers cannot use xDotOld, thus the default is "Use FSAL"=false. In some cases, the xDotOld can be used and save compute cycles. The user can set this for each subStepper.

  • StaggeredForwardSensitivity

    • Not fully updated because constructor needs reworking.

  • CDR test needed modification to handle solving for xDot. Solving for consistent initial conditions (i.e., solving for xDot) can change the Jacobian, especially wrt boundary conditions (BCs). Often BCs are in terms of the solution vector, x. If the boundary condition are included in the Jacobian, solving for xDot could create a singular matrix, as the boundary value, x, is not a function of xDot.

  • Fix Piro Test due to Tempus changes. Piro_TempusSolver_UnitTests.cpp

Checklist

  • My commit messages mention the appropriate GitHub issue numbers.
  • My code follows the code style of the affected package(s).
  • My change requires a change to the documentation.
  • I have updated the documentation accordingly.
  • I have read the code contribution guidelines for this project.
  • I have added tests to cover my changes.
  • All new and existing tests passed.
  • No new compiler warnings were introduced.
  • These changes break backwards compatibility.

Merge request reports