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