diff --git a/hermes2d/src/CMakeLists.txt b/hermes2d/src/CMakeLists.txt index 1d9f1f7e3..5eef9955c 100644 --- a/hermes2d/src/CMakeLists.txt +++ b/hermes2d/src/CMakeLists.txt @@ -78,8 +78,7 @@ set(SRC ${HERMES_COMMON_DIR}/utils.cpp ${HERMES_COMMON_DIR}/matrix.cpp ${HERMES_COMMON_DIR}/tables.cpp - ${HERMES_COMMON_DIR}/Teuchos_stacktrace.cpp - ${HERMES_COMMON_DIR}/butcher_tables.cpp + ${HERMES_COMMON_DIR}/Teuchos_stacktrace.cpp ${HERMES_COMMON_DIR}/solver/nox.cpp ${HERMES_COMMON_DIR}/solver/epetra.cpp ${HERMES_COMMON_DIR}/solver/aztecoo.cpp diff --git a/hermes2d/src/hermes2d.h b/hermes2d/src/hermes2d.h index 51089e376..4a4b42fc2 100644 --- a/hermes2d/src/hermes2d.h +++ b/hermes2d/src/hermes2d.h @@ -37,9 +37,6 @@ // boundary conditions #include "../hermes_common/bctypes.h" -// butcher tables -#include "../../hermes_common/butcher_tables.h" - #include "h2d_common.h" #include "hermes_logging.h" @@ -104,6 +101,8 @@ #include "ogprojection.h" #include "numerical_flux.h" +#include "runge_kutta.h" +#include "tables.h" /** \mainpage diff --git a/hermes2d/src/runge_kutta.cpp b/hermes2d/src/runge_kutta.cpp index 2ac41fd65..665e34568 100644 --- a/hermes2d/src/runge_kutta.cpp +++ b/hermes2d/src/runge_kutta.cpp @@ -178,8 +178,8 @@ void create_stage_wf(double current_time, double time_step, ButcherTable* bt, bool HERMES_RESIDUAL_AS_VECTOR_RK = false; bool rk_time_step(double current_time, double time_step, ButcherTable* const bt, scalar* coeff_vec, DiscreteProblem* dp, MatrixSolverType matrix_solver, - double newton_tol, int newton_max_iter, bool verbose = false, - double newton_damping_coeff = 1.0, double newton_max_allowed_residual_norm = 1e6) + double newton_tol, int newton_max_iter, bool verbose, + double newton_damping_coeff, double newton_max_allowed_residual_norm) { // Matrix and vector for the time derivative part of // the equation (left-hand side). diff --git a/hermes2d/src/runge_kutta.h b/hermes2d/src/runge_kutta.h index 3657c37d4..8a815705b 100644 --- a/hermes2d/src/runge_kutta.h +++ b/hermes2d/src/runge_kutta.h @@ -21,10 +21,10 @@ /// matrix, Y the coefficient vector, and F the (nonlinear) stationary residual. /// Below, "stage_wf_left" and "stage_wf_right" refer to the left-hand side /// and right-hand side of the equation, respectively. -void create_stage_wf(double current_time, double time_step, ButcherTable* bt, +void HERMES_API create_stage_wf(double current_time, double time_step, ButcherTable* bt, DiscreteProblem* dp, WeakForm* stage_wf_right, WeakForm* stage_wf_left); -bool rk_time_step(double current_time, double time_step, ButcherTable* const bt, +bool HERMES_API rk_time_step(double current_time, double time_step, ButcherTable* const bt, scalar* coeff_vec, DiscreteProblem* dp, MatrixSolverType matrix_solver, double newton_tol, int newton_max_iter, bool verbose = false, double newton_damping_coeff = 1.0, double newton_max_allowed_residual_norm = 1e6); diff --git a/hermes2d/tutorial/09-timedep/main.cpp b/hermes2d/tutorial/09-timedep/main.cpp index ff9cd1e8b..dc00beecb 100644 --- a/hermes2d/tutorial/09-timedep/main.cpp +++ b/hermes2d/tutorial/09-timedep/main.cpp @@ -1,7 +1,6 @@ #define HERMES_REPORT_ALL #define HERMES_REPORT_FILE "application.log" #include "hermes2d.h" -#include "runge_kutta.h" // This example shows how to solve a time-dependent PDE discretized // in time via the implicit Euler method. The St. Vitus Cathedral @@ -36,7 +35,7 @@ MatrixSolverType matrix_solver = SOLVER_UMFPACK; // Possibilities: SOLVER_AMESO // Explicit_RK_1, Implicit_RK_1, Explicit_RK_2, Implicit_Crank_Nicolson_2, Implicit_SDIRK_2, // Implicit_Lobatto_IIIA_2, Implicit_Lobatto_IIIB_2, Implicit_Lobatto_IIIC_2, Explicit_RK_3, Explicit_RK_4, // Implicit_Lobatto_IIIA_4, Implicit_Lobatto_IIIB_4, Implicit_Lobatto_IIIC_4. -ButcherTableType butcher_table = Implicit_RK_1; +ButcherTableType butcher_table_type = Implicit_RK_1; // Boundary markers. const std::string BDY_GROUND = "Boundary ground"; @@ -62,7 +61,7 @@ Real temp_ext(Real t) { int main(int argc, char* argv[]) { // Choose a Butcher's table. - ButcherTable *bt = init_butcher_table(butcher_table); + ButcherTable bt(butcher_table_type); // Load the mesh. Mesh mesh; @@ -80,7 +79,7 @@ int main(int argc, char* argv[]) // Enter Dirichlet boundary values. BCValues bc_values; - bc_values.add_const(BDY_GROUND, TEMP_INIT); + bc_values.add_const(BDY_GROUND, TEMP_BND); // Initialize an H1 space with default shapeset. H1Space space(&mesh, &bc_types, &bc_values, P_INIT); @@ -90,6 +89,7 @@ int main(int argc, char* argv[]) // Initialize the solution. Solution u_prev_time; + // Set the initial condition. u_prev_time.set_const(&mesh, TEMP_INIT); @@ -124,9 +124,9 @@ int main(int argc, char* argv[]) { // Perform one Runge-Kutta time step according to the selected Butcher's table. info("Runge-Kutta time step (t = %g, tau = %g, stages: %d).", - current_time, time_step, bt->get_size()); + current_time, time_step, bt.get_size()); bool verbose = true; - if (!rk_time_step(current_time, time_step, bt, coeff_vec, &dp, matrix_solver, + if (!rk_time_step(current_time, time_step, &bt, coeff_vec, &dp, matrix_solver, NEWTON_TOL, NEWTON_MAX_ITER, verbose)) { error("Runge-Kutta time step failed, try to decrease time step size."); } @@ -149,7 +149,6 @@ int main(int argc, char* argv[]) // Cleanup. delete [] coeff_vec; - delete bt; // Wait for the view to be closed. View::wait(); diff --git a/hermes_common/butcher_tables.cpp b/hermes_common/butcher_tables.cpp deleted file mode 100644 index 720a678f1..000000000 --- a/hermes_common/butcher_tables.cpp +++ /dev/null @@ -1,176 +0,0 @@ -// This file is part of Hermes2D. -// -// Hermes2D is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 2 of the License, or -// (at your option) any later version. -// -// Hermes2D is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with Hermes2D. If not, see . - -#include "tables.h" - -ButcherTable* init_butcher_table(ButcherTableType butcher_table) -{ - ButcherTable *bt = new ButcherTable(); - - double gamma = 1./sqrt(2.); - - switch (butcher_table) { - case Explicit_RK_1: // Explicit Euler. - bt->alloc(1); - bt->set_B(0, 1.); - break; - - case Implicit_RK_1: // Implicit Euler. - bt->alloc(1); - bt->set_A(0, 0, 1.); - bt->set_B(0, 1.); - bt->set_C(0, 1.); - break; - - case Explicit_RK_2: // Explicit RK-2. - bt->alloc(2); - bt->set_A(1, 0, 2./3.); - bt->set_A(1, 1, 0.); - bt->set_B(0, 1./4.); - bt->set_B(1, 3./4.); - bt->set_C(0, 2./3.); - break; - - case Implicit_Crank_Nicolson_2: // Implicit Crank Nicolson. - bt->alloc(2); - bt->set_A(0, 0, 1./2.); - bt->set_A(0, 1, 1./2.); - bt->set_B(0, 1./2.); - bt->set_B(1, 1./2.); - bt->set_C(0, 1.); - break; - - case Implicit_SDIRK_2: // Implicit SDIRK-2 (second-order). - bt->alloc(2); - bt->set_A(0, 0, 1. - gamma); - bt->set_A(0, 1, 0.); - bt->set_A(1, 0, gamma); - bt->set_A(1, 1, 1. - gamma); - bt->set_B(0, gamma); - bt->set_B(1, 1. - gamma); - bt->set_C(0, 1. - gamma); - bt->set_C(1, 1.); - break; - - case Implicit_Lobatto_IIIA_2: // Implicit Lobatto IIIA (second-order). - bt->alloc(2); - bt->set_A(1, 0, 1./2.); - bt->set_A(1, 1, 1./2.); - bt->set_B(0, 1./2.); - bt->set_B(1, 1./2.); - bt->set_C(1, 1.); - break; - - case Implicit_Lobatto_IIIB_2: // Implicit Lobatto IIIB (second-order). - bt->alloc(2); - bt->set_A(0, 0, 1./2.); - bt->set_A(0, 1, 1./2.); - bt->set_B(0, 1./2.); - bt->set_B(1, 1./2.); - bt->set_C(0, 1./2.); - bt->set_C(1, 1./2.); - break; - - case Implicit_Lobatto_IIIC_2: // Implicit Lobatto IIIC (second-order). - bt->alloc(2); - bt->set_A(0, 0, 1./2.); - bt->set_A(0, 1, -1./2.); - bt->set_A(1, 0, 1./2.); - bt->set_A(1, 1, 1./2.); - bt->set_B(0, 1./2.); - bt->set_B(1, 1./2.); - bt->set_C(0, 0.); - bt->set_C(1, 1.); - break; - - case Explicit_RK_3: // Explicit RK-3. - bt->alloc(3); - bt->set_A(1, 0, 1./2.); - bt->set_A(2, 0, -1.); - bt->set_A(2, 1, 2.); - bt->set_B(0, 1./6.); - bt->set_B(1, 2./3.); - bt->set_B(2, 1./6.); - bt->set_C(1, 1./2.); - bt->set_C(2, 1.); - break; - - case Explicit_RK_4: // Explicit RK-4. - bt->alloc(4); - bt->set_A(1, 0, 1./2.); - bt->set_A(2, 1, 1./2.); - bt->set_A(3, 2, 1.); - bt->set_B(0, 1./6.); - bt->set_B(1, 1./3.); - bt->set_B(2, 1./3.); - bt->set_B(3, 1./6.); - bt->set_C(1, 1./2.); - bt->set_C(2, 1./2.); - bt->set_C(3, 1.); - break; - - case Implicit_Lobatto_IIIA_4: // Implicit Lobatto IIIA (fourth-order). - bt->alloc(3); - bt->set_A(1, 0, 5./24.); - bt->set_A(2, 0, 1./6.); - bt->set_A(1, 1, 1./3.); - bt->set_A(2, 1, 2./3.); - bt->set_A(1, 2, -1./24.); - bt->set_A(2, 2, 1./6.); - bt->set_B(0, 1./6.); - bt->set_B(1, 2./3.); - bt->set_B(2, 1./6.); - bt->set_C(1, 1./2.); - bt->set_C(2, 1.); - break; - - case Implicit_Lobatto_IIIB_4: // Implicit Lobatto IIIB (fourth-order). - bt->alloc(3); - bt->set_A(0, 0, 1./6.); - bt->set_A(1, 0, 1./6.); - bt->set_A(2, 0, 1./6.); - bt->set_A(0, 1, -1./6.); - bt->set_A(1, 1, 1./3.); - bt->set_A(2, 1, 5./6.); - bt->set_B(0, 1./6.); - bt->set_B(1, 2./3.); - bt->set_B(2, 1./6.); - bt->set_C(1, 1./2.); - bt->set_C(2, 1.); - break; - - case Implicit_Lobatto_IIIC_4: // Implicit Lobatto IIIC (fourth-order). - bt->alloc(3); - bt->set_A(0, 0, 1./6.); - bt->set_A(1, 0, 1./6.); - bt->set_A(2, 0, 1./6.); - bt->set_A(0, 1, -1./3.); - bt->set_A(1, 1, 5./12.); - bt->set_A(2, 1, 2./3.); - bt->set_A(0, 2, 1./6.); - bt->set_A(1, 2, -1./12.); - bt->set_A(2, 2, 1./6.); - bt->set_B(0, 1./6.); - bt->set_B(1, 2./3.); - bt->set_B(2, 1./6.); - bt->set_C(1, 1./2.); - bt->set_C(2, 1.); - break; - - default: error("Unknown Butcher's table."); - } - - return bt; -} diff --git a/hermes_common/butcher_tables.h b/hermes_common/butcher_tables.h deleted file mode 100644 index 2c485e887..000000000 --- a/hermes_common/butcher_tables.h +++ /dev/null @@ -1,25 +0,0 @@ -// This file is part of Hermes2D. -// -// Hermes2D is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 2 of the License, or -// (at your option) any later version. -// -// Hermes2D is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with Hermes2D. If not, see . - -#ifndef __BUTCHER_TABLES_H -#define __BUTCHER_TABLES_H - -#include "common.h" - -class ButcherTable; - -ButcherTable *init_butcher_table(ButcherTableType butcher_table); - -#endif diff --git a/hermes_common/common.h b/hermes_common/common.h index b01531cf1..70174817b 100644 --- a/hermes_common/common.h +++ b/hermes_common/common.h @@ -91,24 +91,6 @@ enum MatrixSolverType SOLVER_AZTECOO }; -// Butcher's tables. -enum ButcherTableType -{ - Explicit_RK_1, // Explicit Runge-Kutta RK-1, or explicit Euler method. - Implicit_RK_1, // Implicit Runge-Kutta RK-1, or implicit Euler method. - Explicit_RK_2, // Explicit Runge-Kutta RK-2 method. - Implicit_Crank_Nicolson_2, // Implicit Crank_Nicolson method. - Implicit_SDIRK_2, // Implicit SDIRK-2 method. - Implicit_Lobatto_IIIA_2, // Implicit Lobatto IIIA-2 method. - Implicit_Lobatto_IIIB_2, // Implicit Lobatto IIIB-2 method. - Implicit_Lobatto_IIIC_2, // Implicit Lobatto IIIB-2 method. - Explicit_RK_3, // Explicit Runge-Kutta RK-3 method. - Explicit_RK_4, // Explicit Runge-Kutta RK-4 method. - Implicit_Lobatto_IIIA_4, // Implicit Lobatto IIIA-4 method. - Implicit_Lobatto_IIIB_4, // Implicit Lobatto IIIB-4 method. - Implicit_Lobatto_IIIC_4 // Implicit Lobatto IIIB-4 method. -}; - // Should be in the same order as MatrixSolverTypes above, so that the // names may be accessed by the same enumeration variable. const std::string MatrixSolverNames[7] = { diff --git a/hermes_common/tables.cpp b/hermes_common/tables.cpp index 0c07209ed..cb2fb8b5d 100644 --- a/hermes_common/tables.cpp +++ b/hermes_common/tables.cpp @@ -77,6 +77,164 @@ ButcherTable::ButcherTable(int size) : Table(size) for (int j=0; jC[j] = 0; } + + +ButcherTable::ButcherTable(ButcherTableType butcher_table) +{ + double gamma = 1./sqrt(2.); + + switch (butcher_table) { + case Explicit_RK_1: // Explicit Euler. + this->alloc(1); + this->set_B(0, 1.); + break; + + case Implicit_RK_1: // Implicit Euler. + this->alloc(1); + this->set_A(0, 0, 1.); + this->set_B(0, 1.); + this->set_C(0, 1.); + break; + + case Explicit_RK_2: // Explicit RK-2. + this->alloc(2); + this->set_A(1, 0, 2./3.); + this->set_A(1, 1, 0.); + this->set_B(0, 1./4.); + this->set_B(1, 3./4.); + this->set_C(0, 2./3.); + break; + + case Implicit_Crank_Nicolson_2: // Implicit Crank Nicolson. + this->alloc(2); + this->set_A(0, 0, 1./2.); + this->set_A(0, 1, 1./2.); + this->set_B(0, 1./2.); + this->set_B(1, 1./2.); + this->set_C(0, 1.); + break; + + case Implicit_SDIRK_2: // Implicit SDIRK-2 (second-order). + this->alloc(2); + this->set_A(0, 0, 1. - gamma); + this->set_A(0, 1, 0.); + this->set_A(1, 0, gamma); + this->set_A(1, 1, 1. - gamma); + this->set_B(0, gamma); + this->set_B(1, 1. - gamma); + this->set_C(0, 1. - gamma); + this->set_C(1, 1.); + break; + + case Implicit_Lobatto_IIIA_2: // Implicit Lobatto IIIA (second-order). + this->alloc(2); + this->set_A(1, 0, 1./2.); + this->set_A(1, 1, 1./2.); + this->set_B(0, 1./2.); + this->set_B(1, 1./2.); + this->set_C(1, 1.); + break; + + case Implicit_Lobatto_IIIB_2: // Implicit Lobatto IIIB (second-order). + this->alloc(2); + this->set_A(0, 0, 1./2.); + this->set_A(0, 1, 1./2.); + this->set_B(0, 1./2.); + this->set_B(1, 1./2.); + this->set_C(0, 1./2.); + this->set_C(1, 1./2.); + break; + + case Implicit_Lobatto_IIIC_2: // Implicit Lobatto IIIC (second-order). + this->alloc(2); + this->set_A(0, 0, 1./2.); + this->set_A(0, 1, -1./2.); + this->set_A(1, 0, 1./2.); + this->set_A(1, 1, 1./2.); + this->set_B(0, 1./2.); + this->set_B(1, 1./2.); + this->set_C(0, 0.); + this->set_C(1, 1.); + break; + + case Explicit_RK_3: // Explicit RK-3. + this->alloc(3); + this->set_A(1, 0, 1./2.); + this->set_A(2, 0, -1.); + this->set_A(2, 1, 2.); + this->set_B(0, 1./6.); + this->set_B(1, 2./3.); + this->set_B(2, 1./6.); + this->set_C(1, 1./2.); + this->set_C(2, 1.); + break; + + case Explicit_RK_4: // Explicit RK-4. + this->alloc(4); + this->set_A(1, 0, 1./2.); + this->set_A(2, 1, 1./2.); + this->set_A(3, 2, 1.); + this->set_B(0, 1./6.); + this->set_B(1, 1./3.); + this->set_B(2, 1./3.); + this->set_B(3, 1./6.); + this->set_C(1, 1./2.); + this->set_C(2, 1./2.); + this->set_C(3, 1.); + break; + + case Implicit_Lobatto_IIIA_4: // Implicit Lobatto IIIA (fourth-order). + this->alloc(3); + this->set_A(1, 0, 5./24.); + this->set_A(2, 0, 1./6.); + this->set_A(1, 1, 1./3.); + this->set_A(2, 1, 2./3.); + this->set_A(1, 2, -1./24.); + this->set_A(2, 2, 1./6.); + this->set_B(0, 1./6.); + this->set_B(1, 2./3.); + this->set_B(2, 1./6.); + this->set_C(1, 1./2.); + this->set_C(2, 1.); + break; + + case Implicit_Lobatto_IIIB_4: // Implicit Lobatto IIIB (fourth-order). + this->alloc(3); + this->set_A(0, 0, 1./6.); + this->set_A(1, 0, 1./6.); + this->set_A(2, 0, 1./6.); + this->set_A(0, 1, -1./6.); + this->set_A(1, 1, 1./3.); + this->set_A(2, 1, 5./6.); + this->set_B(0, 1./6.); + this->set_B(1, 2./3.); + this->set_B(2, 1./6.); + this->set_C(1, 1./2.); + this->set_C(2, 1.); + break; + + case Implicit_Lobatto_IIIC_4: // Implicit Lobatto IIIC (fourth-order). + this->alloc(3); + this->set_A(0, 0, 1./6.); + this->set_A(1, 0, 1./6.); + this->set_A(2, 0, 1./6.); + this->set_A(0, 1, -1./3.); + this->set_A(1, 1, 5./12.); + this->set_A(2, 1, 2./3.); + this->set_A(0, 2, 1./6.); + this->set_A(1, 2, -1./12.); + this->set_A(2, 2, 1./6.); + this->set_B(0, 1./6.); + this->set_B(1, 2./3.); + this->set_B(2, 1./6.); + this->set_C(1, 1./2.); + this->set_C(2, 1.); + break; + + default: error("Unknown Butcher's table."); + } +} + void ButcherTable::alloc(int size) { // Size. diff --git a/hermes_common/tables.h b/hermes_common/tables.h index 1032a9892..3b46036e1 100644 --- a/hermes_common/tables.h +++ b/hermes_common/tables.h @@ -18,10 +18,28 @@ #include "common.h" +// Butcher's tables type. +enum ButcherTableType +{ + Explicit_RK_1, // Explicit Runge-Kutta RK-1, or explicit Euler method. + Implicit_RK_1, // Implicit Runge-Kutta RK-1, or implicit Euler method. + Explicit_RK_2, // Explicit Runge-Kutta RK-2 method. + Implicit_Crank_Nicolson_2, // Implicit Crank_Nicolson method. + Implicit_SDIRK_2, // Implicit SDIRK-2 method. + Implicit_Lobatto_IIIA_2, // Implicit Lobatto IIIA-2 method. + Implicit_Lobatto_IIIB_2, // Implicit Lobatto IIIB-2 method. + Implicit_Lobatto_IIIC_2, // Implicit Lobatto IIIB-2 method. + Explicit_RK_3, // Explicit Runge-Kutta RK-3 method. + Explicit_RK_4, // Explicit Runge-Kutta RK-4 method. + Implicit_Lobatto_IIIA_4, // Implicit Lobatto IIIA-4 method. + Implicit_Lobatto_IIIB_4, // Implicit Lobatto IIIB-4 method. + Implicit_Lobatto_IIIC_4 // Implicit Lobatto IIIB-4 method. +}; + // General square table of real numbers. class HERMES_API Table { - public: +public: Table(); Table(int size); virtual void alloc(int size); @@ -29,7 +47,7 @@ class HERMES_API Table double get_A(int i, int j); void set_A(int i, int j, double val); - protected: +protected: int size; double** A; }; @@ -37,16 +55,17 @@ class HERMES_API Table /// Butcher's tables for Runge-Kutta methods. class HERMES_API ButcherTable: public Table { - public: - ButcherTable(); - ButcherTable(int size); +public: + ButcherTable(); + ButcherTable(int size); + ButcherTable(ButcherTableType butcher_table); virtual void alloc(int size); double get_B(int i); double get_C(int i); void set_B(int i, double val); void set_C(int i, double val); - protected: +protected: double* B; double* C; };