Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Functions
01_Utilize_Thyra.cpp File Reference
#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:
 

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Description:

Example 1: Utilize Thyra

This problem takes 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 Basic Problem tutorial are replaced with snippets using Thyra to create 01_Utilize_Thyra.cpp for the 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
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);
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;
}
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;
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
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;
}
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);
We can also take advantage other Thyra features, Thyra::norm, to check if the solution is passing.
// Test if solution is passing.
if ( std::isnan(x_n[0]) || std::isnan(x_n[1]) ) {
// Test if solution is passing.
if ( std::isnan(Thyra::norm(*x_np1) ) {
The solution update is achieved via Thyra::V_V.
// Promote to next step (n <- n+1).
n++;
x_n[0] = x_np1[0];
x_n[1] = x_np1[1];
// Promote to next step (n <- n+1).
n++;
Thyra::V_V(x_n.ptr(), *x_np1);

The remainder of 01_Utilize_Thyra.cpp (e.g., regression testing) has similar changes to those from above.

Links

Definition at line 205 of file 01_Utilize_Thyra.cpp.