Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_RKButcherTableauHelpers.hpp
1//@HEADER
2// ***********************************************************************
3//
4// Rythmos Package
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25//
26// ***********************************************************************
27//@HEADER
28
29
30#ifndef RYTHMOS_RK_BUTCHER_TABLEAU_HELPERS_HPP
31#define RYTHMOS_RK_BUTCHER_TABLEAU_HELPERS_HPP
32
33#include "Rythmos_Types.hpp"
34
35#include "Rythmos_RKButcherTableauBase.hpp"
36#include "Teuchos_Assert.hpp"
37#include "Thyra_ProductVectorBase.hpp"
38
39namespace Rythmos {
40
41/* \brief . */
42template<class Scalar>
43void assembleIRKState(
44 const int stageIndex,
45 const Teuchos::SerialDenseMatrix<int,Scalar> &A_in,
46 const Scalar dt,
47 const Thyra::VectorBase<Scalar> &x_base,
48 const Thyra::ProductVectorBase<Scalar> &x_stage_bar,
49 Teuchos::Ptr<Thyra::VectorBase<Scalar> > x_out_ptr
50 )
51{
52
53 // typedef ScalarTraits<Scalar> ST; // unused
54
55 const int numStages_in = A_in.numRows();
56 TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( stageIndex, 0, numStages_in );
57 TEUCHOS_ASSERT_EQUALITY( A_in.numRows(), numStages_in );
58 TEUCHOS_ASSERT_EQUALITY( A_in.numCols(), numStages_in );
59 TEUCHOS_ASSERT_EQUALITY( x_stage_bar.productSpace()->numBlocks(), numStages_in );
60 Thyra::VectorBase<Scalar>& x_out = *x_out_ptr;
61
62 V_V( outArg(x_out), x_base );
63 for ( int j = 0; j < numStages_in; ++j ) {
64 Vp_StV( outArg(x_out), dt * A_in(stageIndex,j), *x_stage_bar.getVectorBlock(j) );
65 }
66
67}
68
69
70/* \brief . */
71template<class Scalar>
72void assembleIRKSolution(
73 const Teuchos::SerialDenseVector<int,Scalar> &b_in,
74 const Scalar dt,
75 const Thyra::VectorBase<Scalar> &x_base,
76 const Thyra::ProductVectorBase<Scalar> &x_stage_bar,
77 Teuchos::Ptr<Thyra::VectorBase<Scalar> > x_out_ptr
78 )
79{
80
81 // typedef ScalarTraits<Scalar> ST; // unused
82
83 const int numStages_in = b_in.length();
84 TEUCHOS_ASSERT_EQUALITY( b_in.length(), numStages_in );
85 TEUCHOS_ASSERT_EQUALITY( x_stage_bar.productSpace()->numBlocks(), numStages_in );
86 Thyra::VectorBase<Scalar>& x_out = *x_out_ptr;
87
88 V_V( outArg(x_out), x_base );
89 for ( int j = 0; j < numStages_in; ++j ) {
90 Vp_StV( outArg(x_out), dt * b_in(j), *x_stage_bar.getVectorBlock(j) );
91 }
92
93}
94
95/* \brief . */
96template<class Scalar>
97void assembleERKState(
98 const int stageIndex,
99 const Teuchos::SerialDenseMatrix<int,Scalar> &A_in,
100 const Scalar dt,
101 const Thyra::VectorBase<Scalar> &x_base,
102 const Thyra::VectorBase<Scalar> &x_stage_bar,
103 Teuchos::Ptr<Thyra::VectorBase<Scalar> > x_out_ptr
104 )
105{
106 TEUCHOS_TEST_FOR_EXCEPT(true);
107}
108
109/* \brief . */
110template<class Scalar>
111void assembleERKSolution(
112 const Teuchos::SerialDenseVector<int,Scalar> &b_in,
113 const Scalar dt,
114 const Thyra::VectorBase<Scalar> &x_base,
115 const Thyra::VectorBase<Scalar> &x_stage_bar,
116 Teuchos::Ptr<Thyra::VectorBase<Scalar> > x_out_ptr
117 )
118{
119 TEUCHOS_TEST_FOR_EXCEPT(true);
120}
121
122template<class Scalar>
123bool isEmptyRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt ) {
124 typedef ScalarTraits<Scalar> ST;
125
126 // Check that numStages > 0
127 if (rkbt.numStages() == 0) {
128 return true;
129 }
130
131 // Check that the b vector has _some_ non-zero entry
132 int numNonZero = 0;
133 int numStages_local = rkbt.numStages();
134 const Teuchos::SerialDenseVector<int,Scalar> b_local = rkbt.b();
135 for (int i=0 ; i<numStages_local ; ++i) {
136 if (b_local(i) != ST::zero()) {
137 numNonZero++;
138 }
139 }
140 if (numNonZero == 0) {
141 return true;
142 }
143 // There is no reason to check A and c because they can be zero and you're
144 // producing an explicit method as long as b has something in it.
145 return false;
146}
147
148
149template<class Scalar>
150void assertNonEmptyRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
151{
152 TEUCHOS_TEST_FOR_EXCEPTION( isEmptyRKButcherTableau(rkbt), std::logic_error,
153 "Error, this RKButcherTableau is either empty or the b vector is all zeros!\n"
154 );
155}
156
157template<class Scalar>
158bool isDIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
159{
160 if (isEmptyRKButcherTableau(rkbt)) {
161 return false;
162 }
163 typedef ScalarTraits<Scalar> ST;
164 int numStages_local = rkbt.numStages();
165 const Teuchos::SerialDenseMatrix<int,Scalar> A_local = rkbt.A();
166 for (int i=0 ; i<numStages_local ; ++i) {
167 for (int j=0 ; j<numStages_local ; ++j) {
168 if ((j>i) && (A_local(i,j) != ST::zero())) {
169 return false;
170 }
171 }
172 }
173 return true;
174}
175
176template<class Scalar>
177bool isIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
178{
179 if (isEmptyRKButcherTableau(rkbt)) {
180 return false;
181 }
182 return true;
183}
184
185template<class Scalar>
186void validateIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
187{
188 TEUCHOS_TEST_FOR_EXCEPTION( !isIRKButcherTableau(rkbt), std::logic_error,
189 "Error! This implicit RK Butcher Tableau is empty!\n"
190 );
191}
192
193template<class Scalar>
194void validateDIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
195{
196 TEUCHOS_TEST_FOR_EXCEPTION( !isDIRKButcherTableau(rkbt), std::logic_error,
197 "Error! This Diagonal Implicit RK Butcher Tableau has non-zeros in the upper triangular part!\n"
198 );
199}
200
201template<class Scalar>
202bool isSDIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
203{
204 if (isEmptyRKButcherTableau(rkbt)) {
205 return false;
206 }
207 if (!isDIRKButcherTableau(rkbt)) {
208 return false;
209 }
210 // Verify the diagonal entries are all equal.
211 // typedef ScalarTraits<Scalar> ST; // unused
212 int numStages_local = rkbt.numStages();
213 const Teuchos::SerialDenseMatrix<int,Scalar> A_local = rkbt.A();
214 Scalar val = A_local(0,0);
215 for (int i=0 ; i<numStages_local ; ++i) {
216 if (A_local(i,i) != val) {
217 return false;
218 }
219 }
220 return true;
221}
222
223template<class Scalar>
224void validateSDIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
225{
226 TEUCHOS_TEST_FOR_EXCEPTION( !isSDIRKButcherTableau(rkbt), std::logic_error,
227 "Error! This Singly Diagonal Implicit RK Butcher Tableau does not have equal diagonal entries!\n"
228 );
229}
230
231template<class Scalar>
232bool isERKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt)
233{
234 if (isEmptyRKButcherTableau(rkbt)) {
235 return false;
236 }
237 // Verify the diagonal is zero and the upper triangular part is zero
238 typedef ScalarTraits<Scalar> ST;
239 int numStages_local = rkbt.numStages();
240 const Teuchos::SerialDenseMatrix<int,Scalar> A_local = rkbt.A();
241 for (int i=0 ; i<numStages_local ; ++i) {
242 for (int j=0 ; j<numStages_local ; ++j) {
243 if ((j>=i) && ((A_local(i,j) != ST::zero()))) {
244 return false;
245 }
246 }
247 }
248 const Teuchos::SerialDenseVector<int,Scalar> c_local = rkbt.c();
249 if( c_local(0) != ST::zero() ) {
250 return false;
251 }
252 // 08/13/08 tscoffe: I'm not sure what else I can check for b & c...
253 return true;
254}
255
256
257
258template<class Scalar>
259void validateERKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
260{
261 TEUCHOS_TEST_FOR_EXCEPTION( !isERKButcherTableau(rkbt), std::logic_error,
262 "Error! This ERK Butcher Tableau is not lower triangular or c(0) is not zero!\n"
263 );
264}
265
266/*
267template<class Scalar>
268void validateERKOrder( RKButcherTableauBase<Scalar> rkbt, int order_in )
269{
270 typedef ScalarTraits<Scalar> ST;
271 Teuchos::SerialDenseMatrix<int,Scalar> A_local = rkbt.A();
272 Teuchos::SerialDenseVector<int,Scalar> b_local = rkbt.b();
273 Teuchos::SerialDenseVector<int,Scalar> c_local = rkbt.c();
274 int N = rkbt.numStages();
275 TEUCHOS_TEST_FOR_EXCEPT(N == 0);
276
277 if (order_in == 3) {
278 Scalar sum1 = ST::zero();
279 Scalar sum2 = ST::zero();
280 Scalar sum3 = ST::zero();
281 Scalar sum4 = ST::zero();
282 for (int j=0 ; j<N ; ++j) {
283 sum1 += b_local(j);
284 for (int k=0 ; k<N ; ++k) {
285 sum2 += 2*b_local(j)*A_local(j,k);
286 for (int l=0 ; l<N ; ++l) {
287 sum3 += 3*b_local(j)*A_local(j,k)*A_local(j,l);
288 sum4 += 6*b_local(j)*A_local(j,k)*A_local(k,l);
289 }
290 }
291 }
292 TEUCHOS_TEST_FOR_EXCEPTION(
293 (
294 ( sum1 != ST::one() ) ||
295 ( sum2 != ST::one() ) ||
296 ( sum3 != ST::one() ) ||
297 ( sum4 != ST::one() )
298 ),
299 std::logic_error,
300 "Error!, this RK Butcher Tableau does not meet the order conditions for 3rd order\n"
301 );
302 } else {
303 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,
304 "Error! this function is only defined for order 3\n"
305 );
306 }
307}
308
309template<class Scalar>
310void validateIRKOrder( RKButcherTableauBase<Scalar> rkbt, int order_in )
311{
312 TEUCHOS_TEST_FOR_EXCEPT(true);
313}
314
315template<class Scalar>
316void validateDIRKOrder( RKButcherTableauBase<Scalar> rkbt, int order_in )
317{
318 TEUCHOS_TEST_FOR_EXCEPT(true);
319}
320
321template<class Scalar>
322void validateSDIRKOrder( RKButcherTableauBase<Scalar> rkbt, int order_in )
323{
324 TEUCHOS_TEST_FOR_EXCEPT(true);
325}
326*/
327
328enum E_RKButcherTableauTypes {
329 RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_INVALID,
330 RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_ERK,
331 RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_IRK,
332 RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK,
333 RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK
334};
335
336template<class Scalar>
337E_RKButcherTableauTypes determineRKBTType(const RKButcherTableauBase<Scalar>& rkbt) {
338 if (isEmptyRKButcherTableau(rkbt)) {
339 return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_INVALID;
340 }
341 if (isERKButcherTableau(rkbt)) {
342 return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_ERK;
343 }
344 if (rkbt.numStages() == 1) {
345 // In this case, its just an IRK method.
346 return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_IRK;
347 }
348 if (isSDIRKButcherTableau(rkbt)) {
349 return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK;
350 }
351 if (isDIRKButcherTableau(rkbt)) {
352 return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK;
353 }
354 if (isIRKButcherTableau(rkbt)) {
355 return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_IRK;
356 }
357 return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_INVALID;
358}
359
360
361
362} // namespace Rythmos
363
364
365#endif // RYTHMOS_RK_BUTCHER_TABLEAU_HELPERS_HPP