DOLFIN
DOLFIN C++ interface
Loading...
Searching...
No Matches
PointIntegralSolver.h
1// Copyright (C) 2013 Johan Hake
2//
3// This file is part of DOLFIN.
4//
5// DOLFIN is free software: you can redistribute it and/or modify
6// it under the terms of the GNU Lesser General Public License as published by
7// the Free Software Foundation, either version 3 of the License, or
8// (at your option) any later version.
9//
10// DOLFIN is distributed in the hope that it will be useful,
11// but WITHOUT ANY WARRANTY; without even the implied warranty of
12// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13// GNU Lesser General Public License for more details.
14//
15// You should have received a copy of the GNU Lesser General Public License
16// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
17//
18// First added: 2013-02-15
19// Last changed: 2016-04-27
20
21#ifndef __POINTINTEGRALSOLVER_H
22#define __POINTINTEGRALSOLVER_H
23
24#include <memory>
25#include <set>
26#include <vector>
27
28#include <dolfin/common/Variable.h>
29#include <dolfin/fem/Assembler.h>
30
31namespace dolfin
32{
33 // Forward declarations
34 class MultiStageScheme;
35 class UFC;
36
38
42
44 {
45 public:
46
50 explicit PointIntegralSolver(std::shared_ptr<MultiStageScheme> scheme);
51
54
56 void step(double dt);
57
59 void step_interval(double t0, double t1, double dt);
60
62 std::shared_ptr<MultiStageScheme> scheme() const
63 { return _scheme; }
64
67 {
68 Parameters p("point_integral_solver");
69
70 p.add("reset_stage_solutions", true);
71
72 // Set parameters for NewtonSolver
73 Parameters pn("newton_solver");
74 pn.add("maximum_iterations", 40);
75 pn.add("always_recompute_jacobian", false);
76 pn.add("recompute_jacobian_each_solve", true);
77 pn.add("relaxation_parameter", 1., 0., 1.);
78 pn.add("relative_tolerance", 1e-10, 1e-20, 2.);
79 pn.add("absolute_tolerance", 1e-15, 1e-20, 2.);
80
81 pn.add("kappa", 0.1, 0.05, 1.0);
82 pn.add("eta_0", 1., 1e-15, 1.0);
83 pn.add("max_relative_previous_residual", 1e-1, 1e-5, 1.);
84 pn.add("reset_each_step", true);
85 pn.add("report", false);
86 pn.add("report_vertex", 0, 0, 32767);
87 pn.add("verbose_report", false);
88
89 p.add(pn);
90
91 return p;
92 }
93
96
99
101 std::size_t num_jacobian_computations() const
102 { return _num_jacobian_computations; }
103
104 private:
105
106 // In-place LU factorization of jacobian matrix
107 void _lu_factorize(std::vector<double>& A);
108
109 // Forward backward substitution, assume that mat is already
110 // in place LU factorized
111 void _forward_backward_subst(const std::vector<double>& A,
112 const std::vector<double>& b,
113 std::vector<double>& x) const;
114
115 // Compute jacobian using passed UFC form
116 void _compute_jacobian(std::vector<double>& jac,
117 const std::vector<double>& u,
118 unsigned int local_vert, UFC& loc_ufc,
119 const Cell& cell, const ufc::cell& ufc_cell,
120 int coefficient_index,
121 const std::vector<double>& coordinate_dofs);
122
123 // Compute the norm of a vector
124 double _norm(const std::vector<double>& vec) const;
125
126 // Check the forms making sure they only include piecewise linear
127 // test functions
128 void _check_forms();
129
130 // Build map between vertices, cells and the corresponding local
131 // vertex and initialize UFC data for each form
132 void _init();
133
134 // Solve an explicit stage
135 void _solve_explicit_stage(std::size_t vert_ind, unsigned int stage,
136 const ufc::cell& ufc_cell,
137 const std::vector<double>& coordinate_dofs);
138
139 // Solve an implicit stage
140 void _solve_implicit_stage(std::size_t vert_ind, unsigned int stage,
141 const Cell& cell, const ufc::cell& ufc_cell,
142 const std::vector<double>& coordinate_dofs);
143
144 void
145 _simplified_newton_solve(std::size_t vert_ind, unsigned int stage,
146 const Cell& cell, const ufc::cell& ufc_cell,
147 const std::vector<double>& coordinate_dofs);
148
149 // The MultiStageScheme
150 std::shared_ptr<MultiStageScheme> _scheme;
151
152 // Reference to mesh
153 std::shared_ptr<const Mesh> _mesh;
154
155 // The dofmap (Same for all stages and forms)
156 const GenericDofMap& _dofmap;
157
158 // Size of ODE system
159 const std::size_t _system_size;
160
161 // Offset into local dofmap
162 // FIXME: Consider put in local loop
163 const unsigned int _dof_offset;
164
165 // Number of stages
166 const unsigned int _num_stages;
167
168 // Local to local dofs to be used in tabulate entity dofs
169 std::vector<std::size_t> _local_to_local_dofs;
170
171 // Vertex map between vertices, cells and corresponding local
172 // vertex
173 std::vector<std::pair<std::size_t, unsigned int>> _vertex_map;
174
175 // Local to global dofs used when solution is fanned out to global
176 // vector
177 std::vector<dolfin::la_index> _local_to_global_dofs;
178
179 // Local stage solutions
180 std::vector<std::vector<double>> _local_stage_solutions;
181
182 // Local solutions
183 std::vector<double> _u0;
184 std::vector<double> _residual;
185 std::vector<double> _y;
186 std::vector<double> _dx;
187
188 // UFC objects, one for each form
189 std::vector<std::vector<std::shared_ptr<UFC>>> _ufcs;
190
191 // UFC objects for the last form
192 std::shared_ptr<UFC> _last_stage_ufc;
193
194 // Solution coefficient index in form
195 std::vector<std::vector<int>> _coefficient_index;
196
197 // Flag which is set to false once the jacobian has been computed
198 std::vector<bool> _recompute_jacobian;
199
200 // Jacobians/LU factorized jacobians matrices
201 std::vector<std::vector<double>> _jacobians;
202
203 // Variable used in the estimation of the error of the newton
204 // iteration for the first iteration (important for linear
205 // problems!)
206 double _eta;
207
208 // Number of computations of Jacobian
209 std::size_t _num_jacobian_computations;
210
211 };
212
213}
214
215#endif
A Cell is a MeshEntity of topological codimension 0.
Definition Cell.h:43
This class provides a generic interface for dof maps.
Definition GenericDofMap.h:50
Definition Parameters.h:95
void add(std::string key)
Definition Parameters.h:128
This class is a time integrator for general Runge Kutta forms.
Definition PointIntegralSolver.h:44
std::size_t num_jacobian_computations() const
Return number of computations of jacobian.
Definition PointIntegralSolver.h:101
~PointIntegralSolver()
Destructor.
Definition PointIntegralSolver.cpp:75
PointIntegralSolver(std::shared_ptr< MultiStageScheme > scheme)
Definition PointIntegralSolver.cpp:53
std::shared_ptr< MultiStageScheme > scheme() const
Return the MultiStageScheme.
Definition PointIntegralSolver.h:62
void reset_stage_solutions()
Reset stage solutions.
Definition PointIntegralSolver.cpp:88
void step(double dt)
Step solver with time step dt.
Definition PointIntegralSolver.cpp:102
void reset_newton_solver()
Reset newton solver.
Definition PointIntegralSolver.cpp:80
void step_interval(double t0, double t1, double dt)
Step solver an interval using dt as time step.
Definition PointIntegralSolver.cpp:309
static Parameters default_parameters()
Default parameter values.
Definition PointIntegralSolver.h:66
Definition UFC.h:47
Common base class for DOLFIN variables.
Definition Variable.h:36
Definition adapt.h:30