← Back to branch summary

~chaffra/+junk/trilinos

« back to all changes in this revision

Viewing changes to packages/rythmos/test/UnitTest/Rythmos_ImplicitBDF_UnitTest.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme, Johannes Ring
  • Date: 2009-12-13 12:53:22 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20091213125322-in0nrdjc55deqsw9
Tags: 10.0.3.dfsg-1
[Christophe Prud'homme]
* New upstream release

[Johannes Ring]
* debian/patches/libname.patch: Add prefix 'libtrilinos_' to all
  libraries. 
* debian/patches/soname.patch: Add soversion to libraries.
* debian/watch: Update download URL.
* debian/control:
  - Remove python-numeric from Build-Depends (virtual package).
  - Remove automake and autotools from Build-Depends and add cmake to
    reflect switch to CMake.
  - Add python-support to Build-Depends.
* debian/rules: 
  - Cleanup and updates for switch to CMake.

Show diffs side-by-side

added added

removed removed

Lines of Context:
packages/rythmos/test/UnitTest/Rythmos_ImplicitBDF_UnitTest.cpp
 
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 
23
// USA
 
24
// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
 
25
//
 
26
// ***********************************************************************
 
27
//@HEADER
 
28
 
 
29
#include "Teuchos_UnitTestHarness.hpp"
 
30
 
 
31
#include "Rythmos_Types.hpp"
 
32
#include "Rythmos_UnitTestHelpers.hpp"
 
33
#include "Rythmos_ImplicitBDFStepper.hpp"
 
34
#include "Rythmos_TimeStepNonlinearSolver.hpp"
 
35
#include "../SinCos/SinCosModel.hpp"
 
36
#include "../VanderPol/VanderPolModel.hpp"
 
37
#include "Rythmos_UnitTestModels.hpp"
 
38
#include "Thyra_DetachedVectorView.hpp"
 
39
 
 
40
namespace Rythmos {
 
41
 
 
42
TEUCHOS_UNIT_TEST( Rythmos_ImplicitBDFStepper, minOrder ) {
 
43
  // Model
 
44
  RCP<ParameterList> modelParamList = Teuchos::parameterList();
 
45
  sublist(modelParamList,Stratimikos_name);
 
46
  sublist(modelParamList,DiagonalTransientModel_name);
 
47
  RCP<Thyra::ModelEvaluator<double> > model = getDiagonalModel<double>(modelParamList);
 
48
  Thyra::ModelEvaluatorBase::InArgs<double> model_ic = model->getNominalValues();
 
49
  // Solver
 
50
  RCP<TimeStepNonlinearSolver<double> > nlSolver = timeStepNonlinearSolver<double>();
 
51
  std::vector<StepSizeType> stepTypeVec(2);
 
52
  stepTypeVec[0] = STEP_TYPE_VARIABLE;
 
53
  stepTypeVec[1] = STEP_TYPE_FIXED;
 
54
  double step = -1.0;
 
55
  int maxOrder = 5;
 
56
  for (int minOrder=1 ; minOrder <= 5 ; ++minOrder ) {
 
57
    // parameter list
 
58
    RCP<ParameterList> stepperParamList = Teuchos::parameterList();
 
59
    ParameterList& pl = stepperParamList->sublist("Step Control Settings");
 
60
    pl.set("minOrder",minOrder);
 
61
    pl.set("maxOrder",maxOrder);
 
62
    ParameterList& vopl = pl.sublist("VerboseObject");
 
63
    vopl.set("Verbosity Level","none");
 
64
    for (int i=0 ; i<Teuchos::as<int>(stepTypeVec.size()) ; ++i) {
 
65
      StepSizeType stepType = stepTypeVec[i];
 
66
 
 
67
      // Stepper
 
68
      RCP<ImplicitBDFStepper<double> > stepper = Teuchos::rcp(new ImplicitBDFStepper<double>(model,nlSolver,stepperParamList));
 
69
      TEST_EQUALITY_CONST( Teuchos::is_null(stepper), false );
 
70
      stepper->setInitialCondition(model_ic);
 
71
 
 
72
      for (int order=1 ; order<=minOrder ; ++order) {
 
73
        step = stepper->takeStep(1.0,stepType);
 
74
        const StepStatus<double> status = stepper->getStepStatus();
 
75
        TEST_EQUALITY_CONST( status.order, order );
 
76
      }
 
77
      for (int steps=0 ; steps<4 ; ++steps) {
 
78
        step = stepper->takeStep(1.0,stepType);
 
79
        const StepStatus<double> status = stepper->getStepStatus();
 
80
        TEST_COMPARE( status.order, >=, minOrder );
 
81
        TEST_COMPARE( status.order, <=, maxOrder );
 
82
      }
 
83
    }
 
84
  }
 
85
}
 
86
 
 
87
TEUCHOS_UNIT_TEST( Rythmos_ImplicitBDFStepper, exactNumericalAnswer_BE ) {
 
88
  double a = 1.5;
 
89
  double f = 1.6;
 
90
  double L = 1.7;
 
91
  RCP<ParameterList> modelPL = Teuchos::parameterList();
 
92
  {
 
93
    modelPL->set("Implicit model formulation",true);
 
94
    modelPL->set("Coeff a",a);
 
95
    modelPL->set("Coeff f",f);
 
96
    modelPL->set("Coeff L",L);
 
97
  }
 
98
  RCP<SinCosModel> model = sinCosModel();
 
99
  model->setParameterList(modelPL);
 
100
  Thyra::ModelEvaluatorBase::InArgs<double> model_ic = model->getNominalValues();
 
101
  RCP<TimeStepNonlinearSolver<double> > nlSolver = timeStepNonlinearSolver<double>();
 
102
  RCP<ParameterList> stepperPL = Teuchos::parameterList();
 
103
  {
 
104
    ParameterList& pl = stepperPL->sublist("Step Control Settings");
 
105
    pl.set("minOrder",1);
 
106
    pl.set("maxOrder",1);
 
107
    ParameterList& vopl = pl.sublist("VerboseObject");
 
108
    vopl.set("Verbosity Level","none");
 
109
  }
 
110
  RCP<ImplicitBDFStepper<double> > stepper = implicitBDFStepper<double>(model,nlSolver,stepperPL);
 
111
  stepper->setInitialCondition(model_ic);
 
112
  // Take a few steps and compare the output.
 
113
  int N = 3;
 
114
  double dt = 0.1;
 
115
  double x_exact_0 = 0.0;
 
116
  double x_exact_1 = 1.0;
 
117
  for (int i=1 ; i<=N ; ++i) {
 
118
    double t = i*dt;
 
119
    double dt_taken = stepper->takeStep(dt,STEP_TYPE_FIXED);
 
120
    TEST_ASSERT( dt_taken == dt );
 
121
    RCP<const VectorBase<double> > x;
 
122
    RCP<const VectorBase<double> > xdot;
 
123
    {
 
124
      // Get x out of stepper.
 
125
      Array<double> t_vec;
 
126
      Array<RCP<const VectorBase<double> > > x_vec;
 
127
      Array<RCP<const VectorBase<double> > > xdot_vec;
 
128
      t_vec.resize(1); t_vec[0] = t;
 
129
      stepper->getPoints(t_vec,&x_vec,&xdot_vec,NULL);
 
130
      x = x_vec[0];
 
131
      xdot = xdot_vec[0];
 
132
    }
 
133
    // Compute exact solution:
 
134
    double c = dt/(1+f*f*dt*dt/(L*L));
 
135
    double x_e_0 = c*(x_exact_0/dt + x_exact_1 + dt*f*f/(L*L)*a);
 
136
    double x_e_1 = c*(-f*f/(L*L)*x_exact_0 + x_exact_1/dt + f*f/(L*L)*a);
 
137
    double xd_e_0 = (x_e_0-x_exact_0)/dt;
 
138
    double xd_e_1 = (x_e_1-x_exact_1)/dt;
 
139
    double tol = 1.0e-12;
 
140
    {
 
141
      Thyra::ConstDetachedVectorView<double> x_view( *x );
 
142
      TEST_FLOATING_EQUALITY( x_view[0], x_e_0, tol );
 
143
      TEST_FLOATING_EQUALITY( x_view[1], x_e_1, tol );
 
144
 
 
145
      Thyra::ConstDetachedVectorView<double> xdot_view( *xdot );
 
146
      TEST_FLOATING_EQUALITY( xdot_view[0], xd_e_0, tol );
 
147
      TEST_FLOATING_EQUALITY( xdot_view[1], xd_e_1, tol );
 
148
    }
 
149
    x_exact_0 = x_e_0;
 
150
    x_exact_1 = x_e_1;
 
151
  }
 
152
}
 
153
 
 
154
TEUCHOS_UNIT_TEST( Rythmos_ImplicitBDFStepper, exactNumericalAnswer_BE_nonlinear ) {
 
155
  double epsilon = 0.5;
 
156
  RCP<ParameterList> modelPL = Teuchos::parameterList();
 
157
  {
 
158
    modelPL->set("Implicit model formulation",true);
 
159
    modelPL->set("Coeff epsilon",epsilon);
 
160
  }
 
161
  RCP<VanderPolModel> model = vanderPolModel();
 
162
  model->setParameterList(modelPL);
 
163
  Thyra::ModelEvaluatorBase::InArgs<double> model_ic = model->getNominalValues();
 
164
  RCP<TimeStepNonlinearSolver<double> > nlSolver = timeStepNonlinearSolver<double>();
 
165
  {
 
166
    RCP<ParameterList> nlPL = Teuchos::parameterList();
 
167
    nlPL->set("Default Tol",1.0e-10);
 
168
    nlPL->set("Default Max Iters",20);
 
169
    nlSolver->setParameterList(nlPL);
 
170
  }
 
171
  RCP<ParameterList> stepperPL = Teuchos::parameterList();
 
172
  {
 
173
    ParameterList& pl = stepperPL->sublist("Step Control Settings");
 
174
    pl.set("minOrder",1);
 
175
    pl.set("maxOrder",1);
 
176
    ParameterList& vopl = pl.sublist("VerboseObject");
 
177
    vopl.set("Verbosity Level","none");
 
178
  }
 
179
  RCP<ImplicitBDFStepper<double> > stepper = implicitBDFStepper<double>(model,nlSolver,stepperPL);
 
180
  stepper->setInitialCondition(model_ic);
 
181
  double h = 0.1;
 
182
  std::vector<double> x_0_exact;
 
183
  std::vector<double> x_1_exact;
 
184
  std::vector<double> x_0_dot_exact;
 
185
  std::vector<double> x_1_dot_exact;
 
186
  {
 
187
    x_0_exact.push_back(2.0); // IC
 
188
    x_1_exact.push_back(0.0);
 
189
 
 
190
    x_0_exact.push_back(1.982896621392518e+00); // matlab 
 
191
    x_1_exact.push_back(-1.710337860748234e-01); 
 
192
 
 
193
    x_0_exact.push_back(1.951487185706842e+00); // matlab 
 
194
    x_1_exact.push_back(-3.140943568567556e-01); 
 
195
    
 
196
    x_0_exact.push_back(1.908249109758246e+00); // matlab 
 
197
    x_1_exact.push_back(-4.323807594859574e-01); 
 
198
    
 
199
    x_0_dot_exact.push_back(0.0);
 
200
    x_1_dot_exact.push_back(0.0);
 
201
 
 
202
    for ( int i=1 ; i< Teuchos::as<int>(x_0_exact.size()) ; ++i ) {
 
203
      x_0_dot_exact.push_back( (x_0_exact[i]-x_0_exact[i-1])/h );
 
204
      x_1_dot_exact.push_back( (x_1_exact[i]-x_1_exact[i-1])/h );
 
205
      //std::cout << "x_0_dot_exact["<<i<<"] = "<<x_0_dot_exact[i] << std::endl;
 
206
      //std::cout << "x_1_dot_exact["<<i<<"] = "<<x_1_dot_exact[i] << std::endl;
 
207
    }
 
208
  }
 
209
  double tol_discrete = 1.0e-12;
 
210
  double tol_continuous = 1.0e-2;
 
211
  {
 
212
    // Get IC out
 
213
    double t = 0.0;
 
214
    RCP<const VectorBase<double> > x;
 
215
    RCP<const VectorBase<double> > xdot;
 
216
    {
 
217
      // Get x out of stepper.
 
218
      Array<double> t_vec;
 
219
      Array<RCP<const VectorBase<double> > > x_vec;
 
220
      Array<RCP<const VectorBase<double> > > xdot_vec;
 
221
      t_vec.resize(1); t_vec[0] = t;
 
222
      stepper->getPoints(t_vec,&x_vec,&xdot_vec,NULL);
 
223
      x = x_vec[0];
 
224
      xdot = xdot_vec[0];
 
225
    }
 
226
    {
 
227
      Thyra::ConstDetachedVectorView<double> x_view( *x );
 
228
      TEST_FLOATING_EQUALITY( x_view[0], x_0_exact[0], tol_discrete );
 
229
      TEST_FLOATING_EQUALITY( x_view[1], x_1_exact[0], tol_discrete );
 
230
 
 
231
      Thyra::ConstDetachedVectorView<double> xdot_view( *xdot );
 
232
      TEST_FLOATING_EQUALITY( xdot_view[0], x_0_dot_exact[0], tol_discrete );
 
233
      TEST_FLOATING_EQUALITY( xdot_view[1], x_1_dot_exact[0], tol_discrete );
 
234
    }
 
235
  }
 
236
  for (int i=1 ; i < Teuchos::as<int>(x_0_exact.size()); ++i) {
 
237
    // Each time step
 
238
    double t = 0.0+i*h;
 
239
    double h_taken = stepper->takeStep(h,STEP_TYPE_FIXED);
 
240
    TEST_ASSERT( h_taken == h );
 
241
    RCP<const VectorBase<double> > x;
 
242
    RCP<const VectorBase<double> > xdot;
 
243
    {
 
244
      // Get x out of stepper.
 
245
      Array<double> t_vec;
 
246
      Array<RCP<const VectorBase<double> > > x_vec;
 
247
      Array<RCP<const VectorBase<double> > > xdot_vec;
 
248
      t_vec.resize(1); t_vec[0] = t;
 
249
      stepper->getPoints(t_vec,&x_vec,&xdot_vec,NULL);
 
250
      x = x_vec[0];
 
251
      xdot = xdot_vec[0];
 
252
    }
 
253
    {
 
254
      Thyra::ConstDetachedVectorView<double> x_view( *x );
 
255
      TEST_FLOATING_EQUALITY( x_view[0], x_0_exact[i], tol_discrete );
 
256
      TEST_FLOATING_EQUALITY( x_view[1], x_1_exact[i], tol_discrete );
 
257
 
 
258
      Thyra::ConstDetachedVectorView<double> xdot_view( *xdot );
 
259
      TEST_FLOATING_EQUALITY( xdot_view[0], x_0_dot_exact[i], tol_discrete );
 
260
      TEST_FLOATING_EQUALITY( xdot_view[1], x_1_dot_exact[i], tol_discrete );
 
261
    }
 
262
    // Now compare this to the continuous exact solution:
 
263
    {
 
264
      Thyra::ModelEvaluatorBase::InArgs<double> inArgs = model->getExactSolution(t);
 
265
      RCP<const VectorBase<double> > x_continuous_exact = inArgs.get_x();
 
266
      RCP<const VectorBase<double> > xdot_continuous_exact = inArgs.get_x_dot();
 
267
      {
 
268
        Thyra::ConstDetachedVectorView<double> x_view( *x );
 
269
        Thyra::ConstDetachedVectorView<double> xce_view( *x_continuous_exact );
 
270
        TEST_FLOATING_EQUALITY( x_view[0], xce_view[0], tol_continuous );
 
271
        TEST_FLOATING_EQUALITY( x_view[1], xce_view[1], tol_continuous*10 );
 
272
 
 
273
        Thyra::ConstDetachedVectorView<double> xdot_view( *xdot );
 
274
        Thyra::ConstDetachedVectorView<double> xdotce_view( *xdot_continuous_exact );
 
275
        TEST_FLOATING_EQUALITY( xdot_view[0], xdotce_view[0], tol_continuous*10 );
 
276
        TEST_FLOATING_EQUALITY( xdot_view[1], xdotce_view[1], tol_continuous*10 );
 
277
      }
 
278
    }
 
279
  }
 
280
}
 
281
 
 
282
} // namespace Rythmos
 
283