Plugins
Loading...
Searching...
No Matches
newton_limex_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 * Author: Markus Breit,
4 * but largely copied from ugcore Newton implementation by Andreas Vogel
5 *
6 * This file is part of UG4.
7 *
8 * UG4 is free software: you can redistribute it and/or modify it under the
9 * terms of the GNU Lesser General Public License version 3 (as published by the
10 * Free Software Foundation) with the following additional attribution
11 * requirements (according to LGPL/GPL v3 §7):
12 *
13 * (1) The following notice must be displayed in the Appropriate Legal Notices
14 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
15 *
16 * (2) The following notice must be displayed at a prominent place in the
17 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
18 *
19 * (3) The following bibliography is recommended for citation and must be
20 * preserved in all covered files:
21 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
22 * parallel geometric multigrid solver on hierarchically distributed grids.
23 * Computing and visualization in science 16, 4 (2013), 151-164"
24 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
25 * flexible software system for simulating pde based models on high performance
26 * computers. Computing and visualization in science 16, 4 (2013), 165-179"
27 *
28 * This program is distributed in the hope that it will be useful,
29 * but WITHOUT ANY WARRANTY; without even the implied warranty of
30 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 * GNU Lesser General Public License for more details.
32 */
33
34#include <sstream>
35
37//#include "common/util/string_util.h"
38#include "newton_limex.h"
39
40
41#define PROFILE_NEWTON
42#ifdef PROFILE_NEWTON
43 #define NEWTON_PROFILE_FUNC() PROFILE_FUNC_GROUP("Newton")
44 #define NEWTON_PROFILE_BEGIN(name) PROFILE_BEGIN_GROUP(name, "Newton")
45 #define NEWTON_PROFILE_END() PROFILE_END()
46#else
47 #define NEWTON_PROFILE_FUNC()
48 #define NEWTON_PROFILE_BEGIN(name)
49 #define NEWTON_PROFILE_END()
50#endif
51
52
53namespace ug{
54
55template <typename TAlgebra>
58: m_spLinearSolver(NULL),
59 m_N(NULL),
60 m_J(NULL),
61 m_spAss(NULL),
62 m_linSolverSteps(0),
63 m_linSolverRate(0.0)
64{}
65
66
67template <typename TAlgebra>
70: m_spLinearSolver(NULL),
71 m_N(NULL),
72 m_J(NULL),
73 m_spAss(NULL),
74 m_linSolverSteps(0),
75 m_linSolverRate(0.0)
76{
77 init(N);
78}
79
80
81template <typename TAlgebra>
84: m_spLinearSolver(NULL),
85 m_N(new AssembledOperator<TAlgebra>(spAss)),
86 m_J(NULL),
87 m_spAss(spAss),
88 m_linSolverSteps(0),
89 m_linSolverRate(0.0)
90{}
91
92
93template <typename TAlgebra>
95{
96 NEWTON_PROFILE_BEGIN(NewtonLimexSolver_init);
97 m_N = N.template cast_dynamic<AssembledOperator<TAlgebra> >();
98 if (m_N.invalid())
99 UG_THROW("NewtonLimexSolver: currently only works for AssembledDiscreteOperator.");
100
101 m_spAss = m_N->discretization();
102 return true;
103}
104
105
106template <typename TAlgebra>
108{
109 return true;
110}
111
112
113template <typename TAlgebra>
115{
116 NEWTON_PROFILE_BEGIN(NewtonLimexSolver_apply);
117 // check for linear solver
118 if (m_spLinearSolver.invalid())
119 UG_THROW("NewtonLimexSolver::apply: Linear solver not set.");
120
121 // prepare Jacobian
122 if (m_J.invalid() || m_J->discretization() != m_spAss)
123 m_J = make_sp(new AssembledLinearOperator<TAlgebra>(m_spAss));
124
125 m_J->set_level(m_N->level());
126
127 // create tmp vectors
128 SmartPtr<vector_type> spD = u.clone_without_values();
129 SmartPtr<vector_type> spC = u.clone_without_values();
130
131 // set Dirichlet values
132 try {m_N->prepare(u);}
133 UG_CATCH_THROW("NewtonLimexSolver::prepare: Operator preparation failed.");
134
135 // compute defect
136 try
137 {
138 NEWTON_PROFILE_BEGIN(NewtonComputeDefect1);
139 m_N->apply(*spD, u);
141 }
142 UG_CATCH_THROW("NewtonLimexSolver::apply: Defect computation failed.");
143
144 // increase offset of output for linear solver
145 const int stdLinOffset = m_spLinearSolver->standard_offset();
146 m_spLinearSolver->convergence_check()->set_offset(stdLinOffset + 3);
147
148// perform exactly one Newton step
149 // set c = 0
150 NEWTON_PROFILE_BEGIN(NewtonSetCorretionZero);
151 spC->set(0.0);
153
154 // compute Jacobian
155 try
156 {
157 NEWTON_PROFILE_BEGIN(NewtonComputeJacobian);
158 m_J->init(u);
160 }
161 UG_CATCH_THROW("NewtonLimexSolver::apply: Jacobian initialization failed.");
162
163 // init Jacobian inverse
164 try
165 {
166 NEWTON_PROFILE_BEGIN(NewtonPrepareLinSolver);
167 if (!m_spLinearSolver->init(m_J, u))
168 {
169 UG_LOGN("ERROR in 'NewtonLimexSolver::apply': Cannot init inverse linear "
170 "operator for Jacobi operator.");
171 return false;
172 }
174 }
175 UG_CATCH_THROW("NewtonLimexSolver::apply: Initialization of Linear Solver failed.");
176
177 // solve linearized system
178 try
179 {
180 NEWTON_PROFILE_BEGIN(NewtonApplyLinSolver);
181 if (!m_spLinearSolver->apply(*spC, *spD))
182 {
183 UG_LOGN("ERROR in 'NewtonLimexSolver::apply': Cannot apply inverse linear "
184 "operator for Jacobi operator.");
185 return false;
186 }
188 }
189 UG_CATCH_THROW("NewtonLimexSolver::apply: Application of Linear Solver failed.");
190
191 // store convergence history
192 m_linSolverSteps = m_spLinearSolver->step();
193 m_linSolverRate = m_spLinearSolver->convergence_check()->avg_rate();
194
195 // update solution
196 u -= *spC;
197
198 // apply constraints
199 m_N->prepare(u);
200
201 // reset offset of output for linear solver to previous value
202 m_spLinearSolver->convergence_check()->set_offset(stdLinOffset);
203
204 return true;
205}
206
207
208template <typename TAlgebra>
210{
211 return m_linSolverRate;
212}
213
214template <typename TAlgebra>
216{
217 return m_linSolverSteps;
218}
219
220template <typename TAlgebra>
222{
223 std::stringstream ss;
224 ss << "NewtonLimexSolver\n";
225
226 ss << " LinearSolver: ";
227 if (m_spLinearSolver.valid()) ss << ConfigShift(m_spLinearSolver->config_string()) << "\n";
228 else ss << " NOT SET!\n";
229
230 return ss.str();
231}
232
233
234}
TAlgebra::vector_type vector_type
vector type
Definition newton_limex.h:59
LimexNewtonSolver()
default constructor
Definition newton_limex_impl.h:57
virtual bool prepare(vector_type &u)
prepare operator
Definition newton_limex_impl.h:107
number linear_solver_rate() const
prints average linear solver convergence
Definition newton_limex_impl.h:209
virtual bool init(SmartPtr< IOperator< vector_type > > N)
This operator inverts the operator N: Y -> X.
Definition newton_limex_impl.h:94
virtual std::string config_string() const
Returns information about configuration parameters. This should return necessary information about pa...
Definition newton_limex_impl.h:221
virtual bool apply(vector_type &u)
apply operator, i.e. N^{-1}(0) = u
Definition newton_limex_impl.h:114
int linear_solver_steps() const
information on linear solver convergence
Definition newton_limex_impl.h:215
#define UG_CATCH_THROW(msg)
#define UG_THROW(msg)
#define UG_LOGN(msg)
double number
UG_API std::string ConfigShift(std::string s)
#define NEWTON_PROFILE_END()
Definition newton_limex_impl.h:45
#define NEWTON_PROFILE_BEGIN(name)
Definition newton_limex_impl.h:44
SmartPtr< T, FreePolicy > make_sp(T *inst)
function ProblemDisc new(problemDesc, dom)