From 78b1b9cbad536b8231a9c1862450208aaca65c96 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 24 May 2023 10:49:53 -0600 Subject: [PATCH] Add a test. --- tests/sundials/ida_03.cc | 138 +++++++++++++++++++++++++++++++++++ tests/sundials/ida_03.output | 49 +++++++++++++ tests/sundials/ida_03_in.prm | 19 +++++ 3 files changed, 206 insertions(+) create mode 100644 tests/sundials/ida_03.cc create mode 100644 tests/sundials/ida_03.output create mode 100644 tests/sundials/ida_03_in.prm diff --git a/tests/sundials/ida_03.cc b/tests/sundials/ida_03.cc new file mode 100644 index 0000000000..12735172e9 --- /dev/null +++ b/tests/sundials/ida_03.cc @@ -0,0 +1,138 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2017 - 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//----------------------------------------------------------- + +#include + +#include +#include + +#include + +#include "../tests.h" + + +/** + * Solve the exponential decay problem: + * + * y' = -k y + * y (0) = 1 + * + * with k=2. That is + * + * F(y', y, t) = y' + k y = 0 + * y_0 = 1 + * + * The exact solution is + * + * y(t) = exp(-k*t) + * + * The Jacobian we then need to assemble is the following: + * + * J = dF/dy + alpha dF/dy' + * = k + alpha 1 + * = (k+alpha) + */ +class HarmonicOscillator +{ +public: + HarmonicOscillator( + double _kappa, + const typename SUNDIALS::IDA>::AdditionalData &data) + : time_stepper(data) + , y(1) + , y_dot(1) + , kappa(_kappa) + { + using VectorType = Vector; + + time_stepper.reinit_vector = [&](VectorType &v) { v.reinit(1); }; + + + time_stepper.residual = [&](const double t, + const VectorType &y, + const VectorType &y_dot, + VectorType & res) { + res[0] = y_dot[0] + kappa * y[0]; + }; + + time_stepper.setup_jacobian = [&](const double, + const VectorType &, + const VectorType &, + const double alpha) { + J = kappa + alpha; + }; + + // Used only in ver < 4.0.0 + time_stepper.solve_jacobian_system = + [&](const VectorType &src, VectorType &dst) { dst[0] = src[0] / J; }; + + // Used in ver >= 4.0.0 + time_stepper.solve_with_jacobian = + [&](const VectorType &src, VectorType &dst, const double) { + dst[0] = src[0] / J; + }; + + time_stepper.output_step = [&](const double t, + const VectorType & sol, + const VectorType & sol_dot, + const unsigned int step_number) { + deallog << "Intermediate output:" << std::endl; + deallog << " t =" << t << std::endl; + deallog << " y =" << sol[0] << " (exact: " << std::exp(-kappa * t) + << ')' << std::endl; + deallog << " y'=" << sol_dot[0] + << " (exact: " << -kappa * std::exp(-kappa * t) << ')' + << std::endl; + }; + } + + void + run() + { + y[0] = 1; + y_dot[0] = -kappa; + time_stepper.solve_dae(y, y_dot); + } + SUNDIALS::IDA> time_stepper; + +private: + Vector y; + Vector y_dot; + double J; + double kappa; +}; + + +int +main() +{ + initlog(); + deallog << std::setprecision(10); + + SUNDIALS::IDA>::AdditionalData data; + ParameterHandler prm; + data.add_parameters(prm); + + // std::ofstream ofile(SOURCE_DIR "/ida_03.prm"); + // prm.print_parameters(ofile, ParameterHandler::ShortText); + // ofile.close(); + + std::ifstream ifile(SOURCE_DIR "/ida_03_in.prm"); + prm.parse_input(ifile); + + + HarmonicOscillator ode(2.0, data); + ode.run(); +} diff --git a/tests/sundials/ida_03.output b/tests/sundials/ida_03.output new file mode 100644 index 0000000000..2b3ff26500 --- /dev/null +++ b/tests/sundials/ida_03.output @@ -0,0 +1,49 @@ + +DEAL::Intermediate output: +DEAL:: t =0.000000000 +DEAL:: y =1.000000000 (exact: 1.000000000) +DEAL:: y'=-2.000000000 (exact: -2.000000000) +DEAL::Intermediate output: +DEAL:: t =0.2000000000 +DEAL:: y =0.6703200457 (exact: 0.6703200460) +DEAL:: y'=-1.340640088 (exact: -1.340640092) +DEAL::Intermediate output: +DEAL:: t =0.4000000000 +DEAL:: y =0.4493289638 (exact: 0.4493289641) +DEAL:: y'=-0.8986579256 (exact: -0.8986579282) +DEAL::Intermediate output: +DEAL:: t =0.6000000000 +DEAL:: y =0.3011942117 (exact: 0.3011942119) +DEAL:: y'=-0.6023884233 (exact: -0.6023884238) +DEAL::Intermediate output: +DEAL:: t =0.8000000000 +DEAL:: y =0.2018965179 (exact: 0.2018965180) +DEAL:: y'=-0.4037930380 (exact: -0.4037930360) +DEAL::Intermediate output: +DEAL:: t =1.000000000 +DEAL:: y =0.1353352831 (exact: 0.1353352832) +DEAL:: y'=-0.2706705663 (exact: -0.2706705665) +DEAL::Intermediate output: +DEAL:: t =1.200000000 +DEAL:: y =0.09071795314 (exact: 0.09071795329) +DEAL:: y'=-0.1814359026 (exact: -0.1814359066) +DEAL::Intermediate output: +DEAL:: t =1.400000000 +DEAL:: y =0.06081006262 (exact: 0.06081006263) +DEAL:: y'=-0.1216201252 (exact: -0.1216201253) +DEAL::Intermediate output: +DEAL:: t =1.600000000 +DEAL:: y =0.04076220398 (exact: 0.04076220398) +DEAL:: y'=-0.08152441186 (exact: -0.08152440796) +DEAL::Intermediate output: +DEAL:: t =1.800000000 +DEAL:: y =0.02732372240 (exact: 0.02732372245) +DEAL:: y'=-0.05464744495 (exact: -0.05464744489) +DEAL::Intermediate output: +DEAL:: t =2.000000000 +DEAL:: y =0.01831563874 (exact: 0.01831563889) +DEAL:: y'=-0.03663127757 (exact: -0.03663127778) +DEAL::Intermediate output: +DEAL:: t =2.000000000 +DEAL:: y =0.01831563874 (exact: 0.01831563889) +DEAL:: y'=-0.03663127757 (exact: -0.03663127778) diff --git a/tests/sundials/ida_03_in.prm b/tests/sundials/ida_03_in.prm new file mode 100644 index 0000000000..e6fa7952de --- /dev/null +++ b/tests/sundials/ida_03_in.prm @@ -0,0 +1,19 @@ +set Final time = 2 +set Initial time = 0 +set Time interval between each output = 0.2 +subsection Error control + set Absolute error tolerance = 1e-10 + set Ignore algebraic terms for error computations = true + set Relative error tolerance = 1e-10 +end +subsection Initial condition correction parameters + set Correction type at initial time = none + set Correction type after restart = none + set Maximum number of nonlinear iterations = 10 +end +subsection Running parameters + set Initial step size = 1e-6 + set Maximum number of nonlinear iterations = 10 + set Maximum order of BDF = 5 + set Minimum step size = 1e-7 +end -- 2.39.5