Line data Source code
1 : /*
2 : * Copyright (c) 2010 - 2017: G-CSC, Goethe University Frankfurt
3 : * Author: Arne Naegel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NI_SOLVER__NI_IMPL__
34 : #define __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NI_SOLVER__NI_IMPL__
35 :
36 : #include <iostream>
37 : #include <sstream>
38 :
39 : #include "nested_iteration.h"
40 : #include "lib_disc/function_spaces/grid_function_util.h"
41 : #include "common/util/string_util.h"
42 :
43 : #define PROFILE_NESTED_ITER
44 : #ifdef PROFILE_NESTED_ITER
45 : #define NESTED_ITER_PROFILE_FUNC() PROFILE_FUNC_GROUP("NestedIteration")
46 : #define NESTED_ITER_PROFILE_BEGIN(name) PROFILE_BEGIN_GROUP(name, "NestedIteration")
47 : #define NESTED_ITER_PROFILE_END() PROFILE_END()
48 : #else
49 : #define NESTED_ITER_PROFILE_FUNC()
50 : #define NESTED_ITER_PROFILE_BEGIN(name)
51 : #define NESTED_ITER_PROFILE_END()
52 : #endif
53 :
54 : namespace ug{
55 :
56 : template <typename TDomain, typename TAlgebra>
57 : NestedIterationSolver<TDomain,TAlgebra>::
58 : NestedIterationSolver(SmartPtr<ILinearOperatorInverse<vector_type> > LinearSolver) :
59 : m_spLinearSolver(LinearSolver),
60 : m_N(NULL),
61 : m_J(NULL),
62 : m_spAss(NULL),
63 : m_dgbCall(0),
64 : m_lastNumSteps(0),
65 : m_bUseAdaptiveRefinement(true),
66 : m_maxSteps(100), m_TOL(1e-3), m_absTOL(1e-12)
67 : {};
68 :
69 : template <typename TDomain, typename TAlgebra>
70 0 : NestedIterationSolver<TDomain,TAlgebra>::
71 : NestedIterationSolver() :
72 : m_spLinearSolver(NULL),
73 : m_N(NULL),
74 : m_J(NULL),
75 : m_spAss(NULL),
76 0 : m_dgbCall(0),
77 0 : m_lastNumSteps(0),
78 0 : m_bUseAdaptiveRefinement(true),
79 0 : m_maxSteps(100), m_TOL(1e-3), m_absTOL(1e-12)
80 : {};
81 :
82 : template <typename TDomain, typename TAlgebra>
83 0 : NestedIterationSolver<TDomain,TAlgebra>::
84 : NestedIterationSolver(SmartPtr<IOperator<vector_type> > N) :
85 : m_spLinearSolver(NULL),
86 : m_N(NULL),
87 : m_J(NULL),
88 : m_spAss(NULL),
89 0 : m_dgbCall(0),
90 0 : m_lastNumSteps(0),
91 0 : m_bUseAdaptiveRefinement(true),
92 0 : m_maxSteps(100), m_TOL(1e-3), m_absTOL(1e-12)
93 : {
94 0 : init(N);
95 0 : };
96 :
97 : template <typename TDomain, typename TAlgebra>
98 0 : NestedIterationSolver<TDomain,TAlgebra>::
99 : NestedIterationSolver(SmartPtr<IAssemble<TAlgebra> > spAss) :
100 : m_spLinearSolver(NULL),
101 : m_N(NULL),
102 : m_J(NULL),
103 : m_spAss(spAss),
104 : m_spDomErr(spAss),
105 0 : m_dgbCall(0),
106 0 : m_lastNumSteps(0),
107 0 : m_baseLevel(0),
108 0 : m_bUseAdaptiveRefinement(true),
109 0 : m_maxSteps(100), m_TOL(1e-3), m_absTOL(1e-12)
110 : {
111 0 : m_N = SmartPtr<AssembledOperator<TAlgebra> >(new AssembledOperator<TAlgebra>(m_spAss));
112 0 : };
113 :
114 :
115 : template <typename TDomain, typename TAlgebra>
116 0 : NestedIterationSolver<TDomain,TAlgebra>::
117 : NestedIterationSolver(SmartPtr<IAssemble<TAlgebra> > spAss, SmartPtr<IAssemble<TAlgebra> > spDomErr) :
118 : m_spLinearSolver(NULL),
119 : m_N(NULL),
120 : m_J(NULL),
121 : m_spAss(spAss),
122 : m_spDomErr(spDomErr),
123 0 : m_dgbCall(0),
124 0 : m_lastNumSteps(0),
125 0 : m_baseLevel(0),
126 0 : m_bUseAdaptiveRefinement(true)
127 : {
128 0 : m_N = SmartPtr<AssembledOperator<TAlgebra> >(new AssembledOperator<TAlgebra>(m_spAss));
129 0 : };
130 :
131 :
132 : /*! Requires an AssembledOperator */
133 : template <typename TDomain, typename TAlgebra>
134 0 : bool NestedIterationSolver<TDomain,TAlgebra>::init(SmartPtr<IOperator<vector_type> > N)
135 : {
136 : NESTED_ITER_PROFILE_BEGIN(NestedIterationSolver_init);
137 0 : m_N = N.template cast_dynamic<AssembledOperator<TAlgebra> >();
138 0 : if(m_N.invalid())
139 0 : UG_THROW("NestedIterationSolver: currently only works for AssembledDiscreteOperator.");
140 :
141 0 : m_spAss = m_N->discretization();
142 0 : return true;
143 : }
144 :
145 : //! todo: remove from interface
146 : template <typename TDomain, typename TAlgebra>
147 0 : bool NestedIterationSolver<TDomain,TAlgebra>::prepare(vector_type& u)
148 : {
149 :
150 0 : return true;
151 : }
152 :
153 :
154 : //! Refines domain and provides error estimate.
155 : /*! Values depend on m_spDomErr */
156 : template <typename TDomain, typename TAlgebra>
157 0 : void NestedIterationSolver<TDomain,TAlgebra>::estimate_and_mark_domain(const grid_function_type& u, SmartPtr<IElementMarkingStrategy<TDomain> > spMarking, bool bClearMarks)
158 : {
159 : UG_LOG("Computing error... "<< std::endl);
160 :
161 : typedef DomainDiscretization<TDomain, TAlgebra> domain_disc_type;
162 : // typedef IDomainErrorIndicator<TAlgebra> domain_indicator_type;
163 :
164 0 : SmartPtr<domain_disc_type> spDomainEstimator = m_spDomErr.template cast_dynamic<domain_disc_type>();
165 :
166 : // some checks
167 : UG_ASSERT(spDomainEstimator.valid(), "Huhh: No estimatable object???")
168 : UG_ASSERT(m_spRefiner.valid(), "Huhh: Invalid refiner???");
169 :
170 : // compute error
171 0 : if (m_spElemError.valid()) {
172 : // debug version
173 0 : spDomainEstimator->calc_error(u, u.dd(), m_spElemError.get());
174 : } else {
175 : // standard version
176 0 : spDomainEstimator->calc_error(u, u.dd());
177 : }
178 :
179 : // set (new) marks
180 0 : if (bClearMarks) m_spRefiner->clear_marks();
181 0 : spDomainEstimator->mark_with_strategy(*m_spRefiner, spMarking);
182 : UG_LOG("Estimated error: " << spMarking->global_estimated_error());
183 :
184 :
185 :
186 :
187 0 : }
188 :
189 :
190 : //! Coarsen all elements
191 : template <typename TDomain, typename TAlgebra>
192 : number NestedIterationSolver<TDomain,TAlgebra>::coarsen_domain(const grid_function_type& u)
193 : {
194 : /*typedef typename domain_traits<TDomain::dim>::element_type TElem;
195 : SmartPtr<DoFDistribution> spDD=u.dof_distribution();
196 : m_spRefiner->mark(spDD->begin<TElem>(), spDD->end<TElem>(), RM_COARSEN);
197 : m_spRefiner->coarsen();*/
198 :
199 : return 0; // dummy value
200 : }
201 :
202 :
203 : //! Apply solver for top level grid function
204 : /*! returns: approximation by nested iteration. Input values are ignored.
205 : *
206 : * */
207 : template <typename TDomain, typename TAlgebra>
208 0 : bool NestedIterationSolver<TDomain,TAlgebra>::apply(vector_type& u)
209 : {
210 : #ifdef UG_PARALLEL
211 : // UG_ASSERT(0, "Not implemented!")
212 : #endif
213 : NESTED_ITER_PROFILE_BEGIN(NestedIterationSolver_apply);
214 :
215 : // increase call count
216 0 : m_dgbCall++;
217 0 : grid_function_type &usol = dynamic_cast<grid_function_type&>(u);
218 :
219 : UG_ASSERT (usol.grid_level().is_surface(), "Must supply surface grid function");
220 : UG_ASSERT (usol.grid_level().top(), "Must supply top level grid function");
221 : const GridLevel& surfGridLevel = usol.grid_level();
222 :
223 : // Check for linear solver.
224 0 : if(m_spLinearSolver.invalid())
225 0 : UG_THROW("NestedIterationSolver::apply: Linear Solver not set.");
226 :
227 : // Check for Jacobian.
228 0 : if(m_J.invalid() || m_J->discretization() != m_spAss) {
229 0 : m_J = make_sp(new AssembledLinearOperator<TAlgebra>(m_spAss));
230 : }
231 : m_J->set_level(m_N->level());
232 :
233 0 : UG_LOG("N level: "<<m_N->level()<< std::endl);
234 0 : UG_LOG("u level: "<< usol.grid_level() << std::endl);
235 :
236 : // Create tmp vectors
237 0 : SmartPtr<vector_type> spD = u.clone_without_values();
238 :
239 : std::string ext; //char ext[50];
240 : int loopCnt = 0;
241 0 : m_lastNumSteps = 0;
242 : {
243 : // write start defect for debug
244 : //snprintf(ext, sizeof(ext), "_call%03d", m_dgbCall);
245 0 : std::string name("NESTED_ITER_StartSolution");
246 0 : ext = GetStringPrintf("_call%03d", m_dgbCall);
247 : name.append(ext);
248 0 : write_debug(u, name.c_str());
249 : }
250 :
251 : // Assembly loop (from some coarse level to TOP)
252 : int surfLevel = usol.grid_level().level(); // todo: should start on coarse(r) levels
253 0 : m_topLevel = surfLevel+m_maxSteps;
254 0 : for (int lev=surfLevel; lev<m_topLevel; ++lev)
255 : {
256 : /* OLD LUA CODE:
257 : * globalDisc:assemble_linear(A, b)
258 : * globalDisc:adjust_solution(u)
259 : * solver:init(A)
260 : * solver:apply_return_defect(u,b)
261 : * // globalDisc:adjust_solution(u)
262 : */
263 :
264 : // Assemble.
265 : NESTED_ITER_PROFILE_BEGIN(NestedIterationAssemble);
266 0 : m_spAss->assemble_linear(*m_J, *spD, surfGridLevel); // todo: replace for non-linear problems
267 0 : m_spAss->adjust_solution(u, surfGridLevel);
268 : NESTED_ITER_PROFILE_END();
269 :
270 : //snprintf(ext, sizeof(ext),"_call%03d_iter%03d", m_dgbCall, loopCnt);
271 0 : ext = GetStringPrintf("_call%03d_iter%03d", m_dgbCall, loopCnt);
272 : {
273 : // write initial data for debug
274 0 : std::string name0("NESTED_ITER_InitialSolution"); name0.append(ext);
275 0 : write_debug(u, name0.c_str());
276 0 : std::string name("NESTED_ITER_InitialDefect"); name.append(ext);
277 0 : write_debug(*spD, name.c_str());
278 :
279 : // Write Jacobian for debug
280 0 : std::string matname("NESTED_ITER_Jacobian");
281 : matname.append(ext);
282 0 : write_debug(m_J->get_matrix(), matname.c_str());
283 : }
284 :
285 : // Initialize inverse of Jacobian.
286 : try{
287 : NESTED_ITER_PROFILE_BEGIN(NestedIterationPrepareLinSolver);
288 0 : if(!m_spLinearSolver->init(m_J, u))
289 : {
290 : UG_LOG("ERROR in 'NestedIterationSolver::apply': Cannot init Inverse Linear "
291 : "Operator for Jacobi-Operator.\n");
292 : return false;
293 : }
294 : NESTED_ITER_PROFILE_END();
295 0 : }UG_CATCH_THROW("NestedIterationSolver::apply: Initialization of Linear Solver failed.");
296 :
297 : // Solve linearized system.
298 : try{
299 : NESTED_ITER_PROFILE_BEGIN(NewtonApplyLinSolver);
300 0 : if(!m_spLinearSolver->apply(u, *spD))
301 : {
302 : UG_LOG("ERROR in 'NestedIterationSolver::apply': Cannot apply Inverse Linear "
303 : "Operator for Jacobi-Operator.\n");
304 : return false;
305 : }
306 : NESTED_ITER_PROFILE_END();
307 0 : }UG_CATCH_THROW("NestedIterationSolver::apply: Application of Linear Solver failed.");
308 :
309 : // Adjust solution.
310 0 : m_spAss->adjust_solution(u, surfGridLevel);
311 : {
312 : // write defect for debug
313 0 : std::string name("NESTED_ITER_FinalDefect"); name.append(ext);
314 0 : write_debug(*spD, name.c_str());
315 0 : std::string name3("NESTED_ITER_FinalSolution"); name3.append(ext);
316 0 : write_debug(u, name3.c_str());
317 : }
318 :
319 : // Update counter.
320 0 : loopCnt++;
321 :
322 : // Quit?
323 0 : if (use_adaptive_refinement() == false) {
324 : UG_LOG("SUCCESS: Non-adaptive version!" << std::endl);
325 : break;
326 : }
327 :
328 : // Estimate and mark domain.
329 0 : this->estimate_and_mark_domain(usol, m_spRefinementMarking);
330 :
331 : // OPTIONAL: write defect for debug
332 0 : if (m_spElemError.valid())
333 : {
334 0 : std::string name("NESTED_ITER_Error");
335 : name.append(ext);
336 0 : VTKOutput<TDomain::dim> outError;
337 : outError.print(name.c_str(), *m_spElemError);
338 0 : }
339 :
340 : // Check error.
341 : const number err = m_spRefinementMarking->global_estimated_error();
342 0 : double desiredTOL = (m_spAssociatedSpace.valid()) ? m_TOL*m_spAssociatedSpace->norm(usol) + m_absTOL : m_TOL;
343 :
344 0 : UG_LOG("NI *** Desired tolerance: " << m_TOL << " * "<< m_spAssociatedSpace->norm(usol) << " + "<< m_absTOL <<std::endl);
345 :
346 : // Quit or continue?
347 0 : if(err < desiredTOL)
348 : {
349 : // Quit.
350 : UG_LOG("SUCCESS: Error smaller than tolerance: " << err << " < "<< desiredTOL << std::endl);
351 : break;
352 : } else {
353 : // Refine and continue.
354 : UG_LOG("FAILED: Error (still) greater than tolerance: " << err << " > "<< desiredTOL << std::endl);
355 0 : m_spRefiner->refine();
356 : }
357 : }
358 :
359 :
360 : return true;
361 : }
362 :
363 :
364 : template <typename TDomain, typename TAlgebra>
365 0 : std::string NestedIterationSolver<TDomain,TAlgebra>::config_string() const
366 : {
367 0 : std::stringstream ss;
368 0 : ss << "NestedIterationSolver\n";
369 0 : ss << " LinearSolver: ";
370 0 : if(m_spLinearSolver.valid()) ss << ConfigShift(m_spLinearSolver->config_string()) << "\n";
371 0 : else ss << " NOT SET!\n";
372 0 : return ss.str();
373 0 : }
374 :
375 :
376 : }
377 :
378 : #endif /* __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NESTED_ITER_SOLVER__NESTED_ITER_IMPL__ */
379 :
|
