Loading...
Searching...
No Matches
PointIntegralSolver.h
// Copyright (C) 2013 Johan Hake
//
// This file is part of DOLFIN.
//
// DOLFIN is free software: you can 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 3 of the License, or
// (at your option) any later version.
//
// DOLFIN 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 Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
// First added: 2013-02-15
// Last changed: 2016-04-27
#ifndef __POINTINTEGRALSOLVER_H
#define __POINTINTEGRALSOLVER_H
#include <memory>
#include <set>
#include <vector>
#include <dolfin/common/Variable.h>
#include <dolfin/fem/Assembler.h>
namespace dolfin
{
// Forward declarations
class MultiStageScheme;
class UFC;
class PointIntegralSolver : public Variable
{
public:
explicit PointIntegralSolver(std::shared_ptr<MultiStageScheme> scheme);
~PointIntegralSolver();
void step(double dt);
void step_interval(double t0, double t1, double dt);
std::shared_ptr<MultiStageScheme> scheme() const
{ return _scheme; }
static Parameters default_parameters()
{
Parameters p("point_integral_solver");
p.add("reset_stage_solutions", true);
// Set parameters for NewtonSolver
Parameters pn("newton_solver");
pn.add("maximum_iterations", 40);
pn.add("always_recompute_jacobian", false);
pn.add("recompute_jacobian_each_solve", true);
pn.add("relaxation_parameter", 1., 0., 1.);
pn.add("relative_tolerance", 1e-10, 1e-20, 2.);
pn.add("absolute_tolerance", 1e-15, 1e-20, 2.);
pn.add("kappa", 0.1, 0.05, 1.0);
pn.add("eta_0", 1., 1e-15, 1.0);
pn.add("max_relative_previous_residual", 1e-1, 1e-5, 1.);
pn.add("reset_each_step", true);
pn.add("report", false);
pn.add("report_vertex", 0, 0, 32767);
pn.add("verbose_report", false);
p.add(pn);
return p;
}
void reset_newton_solver();
void reset_stage_solutions();
std::size_t num_jacobian_computations() const
{ return _num_jacobian_computations; }
private:
// In-place LU factorization of jacobian matrix
void _lu_factorize(std::vector<double>& A);
// Forward backward substitution, assume that mat is already
// in place LU factorized
void _forward_backward_subst(const std::vector<double>& A,
const std::vector<double>& b,
std::vector<double>& x) const;
// Compute jacobian using passed UFC form
void _compute_jacobian(std::vector<double>& jac,
const std::vector<double>& u,
unsigned int local_vert, UFC& loc_ufc,
const Cell& cell, const ufc::cell& ufc_cell,
int coefficient_index,
const std::vector<double>& coordinate_dofs);
// Compute the norm of a vector
double _norm(const std::vector<double>& vec) const;
// Check the forms making sure they only include piecewise linear
// test functions
void _check_forms();
// Build map between vertices, cells and the corresponding local
// vertex and initialize UFC data for each form
void _init();
// Solve an explicit stage
void _solve_explicit_stage(std::size_t vert_ind, unsigned int stage,
const ufc::cell& ufc_cell,
const std::vector<double>& coordinate_dofs);
// Solve an implicit stage
void _solve_implicit_stage(std::size_t vert_ind, unsigned int stage,
const Cell& cell, const ufc::cell& ufc_cell,
const std::vector<double>& coordinate_dofs);
void
_simplified_newton_solve(std::size_t vert_ind, unsigned int stage,
const Cell& cell, const ufc::cell& ufc_cell,
const std::vector<double>& coordinate_dofs);
// The MultiStageScheme
std::shared_ptr<MultiStageScheme> _scheme;
// Reference to mesh
std::shared_ptr<const Mesh> _mesh;
// The dofmap (Same for all stages and forms)
const GenericDofMap& _dofmap;
// Size of ODE system
const std::size_t _system_size;
// Offset into local dofmap
// FIXME: Consider put in local loop
const unsigned int _dof_offset;
// Number of stages
const unsigned int _num_stages;
// Local to local dofs to be used in tabulate entity dofs
std::vector<std::size_t> _local_to_local_dofs;
// Vertex map between vertices, cells and corresponding local
// vertex
std::vector<std::pair<std::size_t, unsigned int>> _vertex_map;
// Local to global dofs used when solution is fanned out to global
// vector
std::vector<dolfin::la_index> _local_to_global_dofs;
// Local stage solutions
std::vector<std::vector<double>> _local_stage_solutions;
// Local solutions
std::vector<double> _u0;
std::vector<double> _residual;
std::vector<double> _y;
std::vector<double> _dx;
// UFC objects, one for each form
std::vector<std::vector<std::shared_ptr<UFC>>> _ufcs;
// UFC objects for the last form
std::shared_ptr<UFC> _last_stage_ufc;
// Solution coefficient index in form
std::vector<std::vector<int>> _coefficient_index;
// Flag which is set to false once the jacobian has been computed
std::vector<bool> _recompute_jacobian;
// Jacobians/LU factorized jacobians matrices
std::vector<std::vector<double>> _jacobians;
// Variable used in the estimation of the error of the newton
// iteration for the first iteration (important for linear
// problems!)
double _eta;
// Number of computations of Jacobian
std::size_t _num_jacobian_computations;
};
}
#endif
dolfin::Cell
A Cell is a MeshEntity of topological codimension 0.
Definition Cell.h:43
dolfin::GenericDofMap
This class provides a generic interface for dof maps.
Definition GenericDofMap.h:50
dolfin::MultiStageScheme
Place-holder for forms and solutions for a multi-stage Butcher tableau based method.
Definition MultiStageScheme.h:42
dolfin::Parameters
Definition Parameters.h:95
dolfin::Parameters::add
void add(std::string key)
Definition Parameters.h:128
dolfin::PointIntegralSolver::num_jacobian_computations
std::size_t num_jacobian_computations() const
Return number of computations of jacobian.
Definition PointIntegralSolver.h:101
dolfin::PointIntegralSolver::~PointIntegralSolver
~PointIntegralSolver()
Destructor.
Definition PointIntegralSolver.cpp:75
dolfin::PointIntegralSolver::PointIntegralSolver
PointIntegralSolver(std::shared_ptr< MultiStageScheme > scheme)
Definition PointIntegralSolver.cpp:53
dolfin::PointIntegralSolver::scheme
std::shared_ptr< MultiStageScheme > scheme() const
Return the MultiStageScheme.
Definition PointIntegralSolver.h:62
dolfin::PointIntegralSolver::reset_stage_solutions
void reset_stage_solutions()
Reset stage solutions.
Definition PointIntegralSolver.cpp:88
dolfin::PointIntegralSolver::step
void step(double dt)
Step solver with time step dt.
Definition PointIntegralSolver.cpp:102
dolfin::PointIntegralSolver::reset_newton_solver
void reset_newton_solver()
Reset newton solver.
Definition PointIntegralSolver.cpp:80
dolfin::PointIntegralSolver::step_interval
void step_interval(double t0, double t1, double dt)
Step solver an interval using dt as time step.
Definition PointIntegralSolver.cpp:309
dolfin::PointIntegralSolver::default_parameters
static Parameters default_parameters()
Default parameter values.
Definition PointIntegralSolver.h:66
dolfin::UFC
Definition UFC.h:47
dolfin::Variable::Variable
Variable()
Create unnamed variable.
Definition Variable.cpp:31
dolfin
Definition adapt.h:30