Tempus
Version of the Day
Time Integration
|
#include <iomanip>
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include "Teuchos_StandardCatchMacros.hpp"
#include "Thyra_VectorStdOps.hpp"
#include "Thyra_DefaultSpmdVectorSpace.hpp"
#include "Thyra_DetachedVectorView.hpp"
Go to the source code of this file.
Functions | |
int | main (int argc, char *argv[]) |
Description: More... | |
int main | ( | int | argc, |
char * | argv[] | ||
) |
Description:
This problem takes 00_Basic_Problem.cpp "Basic Problem" and utilizes Thyra vectors, replacing the C++ double arrays, 01_Utilize_Thyra.cpp. Tempus uses Thyra for its Abstract Numerical Algorithms (ANAs), which are the mathematical concepts of vectors, vector spaces, and linear operators. All other ANA interfaces and support software are built on these fundamental operator/vector interfaces (including ModelEvaluators).
In the following table, code snippets from the 00_Basic_Problem.cpp "Basic Problem" tutorial are replaced with snippets using Thyra to create 01_Utilize_Thyra.cpp for the 01_Utilize_Thyra.cpp "Utilize Thyra" tutorial. This is similar to performing a diff between 00_Basic_Problem.cpp and 01_Utilize_Thyra.cpp, but the first column provides comments related to the changes. You may want to to do a diff (e.g., vimdiff or tkdiff) to see these changes within main (e.g., vimdiff 00_Basic_Problem/00_Basic_Problem.cpp 01_Utilize_Thyra/01_Utilize_Thyra.cpp).
Comments | Original "Basic Problem" Code Snippet | Replacement "Utilize Thyra" Code Snippet |
---|---|---|
Setup Thyra vectors. We first need to replace the C++ double arrays with a vector space to construct the Thyra::Vector. We additionally have introduced the use of Teuchos Reference-Counted Pointers (Teuchos:RCP), which are Trilinos's smart pointers. Details on RCP can be found at https://www.osti.gov/servlets/purl/919177. | // Solution and its time-derivative.
double x_n[2]; // at time index n
double xDot_n[2]; // at time index n
| // Solution and its time-derivative.
int vectorLength = 2; // number state unknowns
RCP<const Thyra::VectorSpaceBase<double> > xSpace =
Thyra::defaultSpmdVectorSpace<double>(vectorLength);
RCP<Thyra::VectorBase<double> > x_n = Thyra::createMember(xSpace);
RCP<Thyra::VectorBase<double> > xDot_n = Thyra::createMember(xSpace);
|
Initialize Thyra vectors. The initialization can be achieved with the Thyra::DetachedVectorView's. The scoping ensures they are deleted after use. | // Initial Conditions
double time = 0.0;
double epsilon = 1.0e-1;
x_n [0] = 2.0;
x_n [1] = 0.0;
xDot_n[0] = 0.0;
xDot_n[1] = -2.0/epsilon;
| // Initial Conditions
double time = 0.0;
double epsilon = 1.0e-1;
{ // Scope to delete DetachedVectorViews
Thyra::DetachedVectorView<double> x_n_view(*x_n);
x_n_view[0] = 2.0;
x_n_view[1] = 0.0;
Thyra::DetachedVectorView<double> xDot_n_view(*xDot_n);
xDot_n_view[0] = 0.0;
xDot_n_view[1] = -2.0/epsilon;
}
|
Access Thyra vectors. Elements of the Thyra::Vector can be quickly accessed with get_ele. | cout << n << " " << time << " " << x_n[0] << " " << x_n[1] << endl;
| cout << n << " " << time << " " << get_ele(*(x_n), 0)
<< " " << get_ele(*(x_n), 1) << endl;
|
Clone Thyra vectors. To initialize the solution at the next time step, we can simply clone the current timestep. | // Initialize next time step
double x_np1[2]; // at time index n+1
x_np1[0] = x_n[0];
x_np1[1] = x_n[1];
| // Initialize next time step
RCP<Thyra::VectorBase<double> > x_np1 = x_n->clone_v(); // at time index n+1
|
Accessing elements of Thyra vectors. The model evaluation is achieved through Thyra::DetachedVectorView. | // Righthand side evaluation and time-derivative at n.
xDot_n[0] = x_n[1];
xDot_n[1] = ((1.0 - x_n[0]*x_n[0])*x_n[1] - x_n[0])/epsilon;
| // Righthand side evaluation and time-derivative at n.
{
Thyra::ConstDetachedVectorView<double> x_n_view(*x_n);
Thyra::DetachedVectorView<double> xDot_n_view(*xDot_n);
xDot_n_view[0] = x_n_view[1];
xDot_n_view[1] =
((1.0-x_n_view[0]*x_n_view[0])*x_n_view[1]-x_n_view[0])/epsilon;
}
|
Compute with Thyra vectors. The Forward Euler time stepping is achieved by using Thyra::V_VpStV, which performs an axpy. | // Take the timestep - Forward Euler
x_np1[0] = x_n[0] + dt*xDot_n[0];
x_np1[1] = x_n[1] + dt*xDot_n[1];
| // Take the timestep - Forward Euler
Thyra::V_VpStV(x_np1.ptr(), *x_n, dt, *xDot_n);
|
Other Thyra vector utilities. We can also take advantage other Thyra features, Thyra::norm, to check if the solution has passed. | // Test if solution has passed.
if ( std::isnan(x_n[0]) || std::isnan(x_n[1]) ) {
| // Test if solution has passed.
if ( std::isnan(Thyra::norm(*x_np1) ) {
|
Use Thyra vector assignment. The solution update is achieved via Thyra::V_V. | // Promote to next step (n <- n+1).
x_n[0] = x_np1[0];
x_n[1] = x_np1[1];
n++;
| // Promote to next step (n <- n+1).
Thyra::V_V(x_n.ptr(), *x_np1);
n++;
|
The remainder of 01_Utilize_Thyra.cpp (e.g., regression testing) has similar changes to those from above.
Definition at line 214 of file 01_Utilize_Thyra.cpp.