CppNoddy  0.92
All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Modules Pages
BVPKarmanJacobian.cpp
Go to the documentation of this file.
1/// \file BVPKarmanJacobian.cpp
2/// \ingroup Tests
3/// \ingroup BVP
4/// Solving the Karman rotating-disk equations for the flow above an
5/// infinite rotating disk by applying a user-provided Jacobian
6/// \f[ U''(y) = U^2(y) + V(y)U'(y) - W^2(y) \f]
7/// \f[ W''(y) = 2U(y)W(y) + V(y)W'(y) \f]
8/// \f[ 2U(y) + V'(y) = 0 \f]
9/// with boundary conditions \f$ U(0)=V(0)=0 \f$, \f$ W(0)=1 \f$
10/// and \f$ U(\infty ) \to 0 \f$, \f$ W(\infty ) \to 0 \f$.
11/// The class constructs and solves the
12/// global matrix problem using 2nd-order finite differences.
13
14#include <BVP_bundle.h>
15
16#include "../Utils_Fill.h"
17
18// enumerate the variables in the ODE system
19enum { U, Ud, V, W, Wd };
20namespace CppNoddy
21{
22 namespace Example
23 {
24 /// Define the Karman equations
25 class Karman_equations : public Equation_1matrix<double>
26 {
27 public:
28
29 /// The Karman system is a 5th order real system of ODEs
31
32 /// Define the Karman system
34 {
35 // The 5th order system for ( U, U', V, W, W' )
36 f[ U ] = z[ Ud ];
37 f[ Ud ] = z[ U ] * z[ U ] + z[ V ] * z[ Ud ] - z[ W ] * z[ W ];
38 f[ V ] = -2 * z[ U ];
39 f[ W ] = z[ Wd ];
40 f[ Wd ] = 2 * z[ U ] * z[ W ] + z[ V ] * z[ Wd ];
41 }
42
43 /// Provide the exact Jacobian of above RHS rather than using finite-differences
45 {
46 jac( 0, Ud ) = 1.0;
47 jac( 1, U ) = 2 * z[ U ];
48 jac( 1, Ud ) = z[ V ];
49 jac( 1, V ) = z[ Ud ];
50 jac( 1, W ) = -2 * z[ W ];
51 jac( 2, U ) = -2.0;
52 jac( 3, Wd ) = 1.0;
53 jac( 4, U ) = 2 * z[ W ];
54 jac( 4, V ) = z[ Wd ];
55 jac( 4, W ) = 2 * z[ U ];
56 jac( 4, Wd ) = z[ V ];
57 }
58
60 {
62 }
63 };
64
65 class Karman_left_BC : public Residual<double>
66 {
67 public:
68 // 3 BCs for 5 unknowns
69 Karman_left_BC() : Residual<double> ( 3, 5 ) {}
70
72 {
73 B[ 0 ] = z[ U ];
74 B[ 1 ] = z[ V ];
75 B[ 2 ] = z[ W ] - 1.0;
76 }
77 };
78
79 class Karman_right_BC : public Residual<double>
80 {
81 public:
82 // 2 BCs for 5 unknowns
83 Karman_right_BC() : Residual<double> ( 2, 5 ) {}
84
86 {
87 B[ 0 ] = z[ U ];
88 B[ 1 ] = z[ W ];
89 }
90 };
91
92
93 } // end Example namespace
94} // end CppNoddy namespace
95
96using namespace CppNoddy;
97using namespace std;
98
99int main()
100{
101 cout << "\n";
102 cout << "=== BVP: Karman equations with analytic Jacobian ====\n";
103 cout << "\n";
104
108
109 // domain for the ODE problem
110 double left = 0.0;
111 double right = 20.0;
112 // number of nodal points
113 int N = 801;
114 // mesh
115 DenseVector<double> nodes( Utility::power_node_vector( left, right, N, 1.2 ) );
116
117 // pass it to the ode
118 ODE_BVP<double> ode( &problem, nodes, &BC_left, &BC_right );
119
120 // set the solution to the initial guess
121 for ( int i = 0; i < N; ++i )
122 {
123 double y = ode.solution().coord( i );
124 ode.solution()( i, U ) = 0.0;
125 ode.solution()( i, Ud ) = 0.0;
126 ode.solution()( i, V ) = 0.0;
127 ode.solution()( i, W ) = exp( -y );
128 ode.solution()( i, Wd ) = -exp( -y );
129 }
130
131 try
132 {
133 ode.solve2();
134 }
135 catch ( const std::runtime_error &error )
136 {
137 cout << " \033[1;31;48m * FAILED THROUGH EXCEPTION BEING RAISED \033[0m\n";
138 return 1;
139 }
140
141 const double tol( 1.e-4 );
142 // check the BL transpiration vs the known solution
143 if ( abs( ode.solution()( N - 1, V ) + 0.88447 ) > tol )
144 {
145 cout << "\033[1;31;48m * FAILED \033[0m\n";
146 cout << " Difference = " << abs( ode.solution()( N - 1, V ) + 0.88447 ) << "\n";
147 return 1;
148 }
149 else
150 {
151 cout << "\033[1;32;48m * PASSED \033[0m\n";
152 return 0;
153 }
154
155}
@ f
Definition: BVPBerman.cpp:15
int main()
@ V
Definition: BVPKarman.cpp:20
@ Wd
Definition: BVPKarman.cpp:20
@ W
Definition: BVPKarman.cpp:20
@ U
Definition: BVPKarman.cpp:20
@ Ud
Definition: BVPKarman.cpp:20
A shorter bundled include file for ODE_BVP and PDE_IBVP codes.
A matrix class that constructs a DENSE matrix as a row major std::vector of DenseVectors.
Definition: DenseMatrix.h:25
An DenseVector class – a dense vector object.
Definition: DenseVector.h:34
An equation object base class used in the IBVP classes (and others).
Define the Karman equations.
Definition: BVPKarman.cpp:28
Karman_equations()
The Karman system is a 5th order real system of ODEs.
void jacobian(const DenseVector< double > &z, DenseMatrix< double > &jac) const
Provide the exact Jacobian of above RHS rather than using finite-differences.
void matrix0(const DenseVector< double > &x, DenseMatrix< double > &m) const
Define the matrix in terms of the current state vector.
void residual_fn(const DenseVector< double > &z, DenseVector< double > &f) const
Define the Karman system.
Define the boundary conditions.
Definition: BVPKarman.cpp:54
void residual_fn(const DenseVector< double > &z, DenseVector< double > &B) const
A blank virtual residual function method.
void residual_fn(const DenseVector< double > &z, DenseVector< double > &B) const
A blank virtual residual function method.
A templated object for real/complex vector system of first-order ordinary differential equations.
Definition: ODE_BVP.h:37
OneD_Node_Mesh< _Type, _Xtype > & solution()
Definition: ODE_BVP.h:187
void solve2()
Formulate and solve the ODE using Newton iteration and a second-order finite difference scheme.
Definition: ODE_BVP.cpp:83
A base class to be inherited by objects that define residuals.
Definition: Residual.h:15
DenseVector< double > power_node_vector(const double &lower, const double &upper, const std::size_t &N, const double &power)
Return a DENSE vector with the nodal points of a non-uniform mesh distributed between the upper/lower...
Definition: Utility.cpp:123
A collection of OO numerical routines aimed at simple (typical) applied problems in continuum mechani...
void fill_identity(CppNoddy::Sequential_Matrix_base< _Type > &A)
Fill diagonal with unit values.
Definition: Utils_Fill.h:22

© 2012

R.E. Hewitt