Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
nested_iteration_impl.h
Go to the documentation of this file.
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"
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
54namespace ug{
55
56template <typename TDomain, typename TAlgebra>
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
69template <typename TDomain, typename TAlgebra>
72 m_spLinearSolver(NULL),
73 m_N(NULL),
74 m_J(NULL),
75 m_spAss(NULL),
76 m_dgbCall(0),
77 m_lastNumSteps(0),
78 m_bUseAdaptiveRefinement(true),
79 m_maxSteps(100), m_TOL(1e-3), m_absTOL(1e-12)
80{};
81
82template <typename TDomain, typename TAlgebra>
85 m_spLinearSolver(NULL),
86 m_N(NULL),
87 m_J(NULL),
88 m_spAss(NULL),
89 m_dgbCall(0),
90 m_lastNumSteps(0),
91 m_bUseAdaptiveRefinement(true),
92 m_maxSteps(100), m_TOL(1e-3), m_absTOL(1e-12)
93{
94 init(N);
95};
96
97template <typename TDomain, typename TAlgebra>
100 m_spLinearSolver(NULL),
101 m_N(NULL),
102 m_J(NULL),
103 m_spAss(spAss),
104 m_spDomErr(spAss),
105 m_dgbCall(0),
106 m_lastNumSteps(0),
107 m_baseLevel(0),
108 m_bUseAdaptiveRefinement(true),
109 m_maxSteps(100), m_TOL(1e-3), m_absTOL(1e-12)
110{
112};
113
114
115template <typename TDomain, typename TAlgebra>
118 m_spLinearSolver(NULL),
119 m_N(NULL),
120 m_J(NULL),
121 m_spAss(spAss),
122 m_spDomErr(spDomErr),
123 m_dgbCall(0),
124 m_lastNumSteps(0),
125 m_baseLevel(0),
126 m_bUseAdaptiveRefinement(true)
127{
129};
130
131
133template <typename TDomain, typename TAlgebra>
135{
136 NESTED_ITER_PROFILE_BEGIN(NestedIterationSolver_init);
137 m_N = N.template cast_dynamic<AssembledOperator<TAlgebra> >();
138 if(m_N.invalid())
139 UG_THROW("NestedIterationSolver: currently only works for AssembledDiscreteOperator.");
140
141 m_spAss = m_N->discretization();
142 return true;
143}
144
146template <typename TDomain, typename TAlgebra>
148{
149
150 return true;
151}
152
153
155
156template <typename TDomain, typename TAlgebra>
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 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 if (m_spElemError.valid()) {
172 // debug version
173 spDomainEstimator->calc_error(u, u.dd(), m_spElemError.get());
174 } else {
175 // standard version
176 spDomainEstimator->calc_error(u, u.dd());
177 }
178
179 // set (new) marks
180 if (bClearMarks) m_spRefiner->clear_marks();
181 spDomainEstimator->mark_with_strategy(*m_spRefiner, spMarking);
182 UG_LOG("Estimated error: " << spMarking->global_estimated_error());
183
184
185
186
187}
188
189
191template <typename TDomain, typename TAlgebra>
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
204
207template <typename TDomain, typename TAlgebra>
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 m_dgbCall++;
217 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 if(m_spLinearSolver.invalid())
225 UG_THROW("NestedIterationSolver::apply: Linear Solver not set.");
226
227 // Check for Jacobian.
228 if(m_J.invalid() || m_J->discretization() != m_spAss) {
229 m_J = make_sp(new AssembledLinearOperator<TAlgebra>(m_spAss));
230 }
231 m_J->set_level(m_N->level());
232
233 UG_LOG("N level: "<<m_N->level()<< std::endl);
234 UG_LOG("u level: "<< usol.grid_level() << std::endl);
235
236 // Create tmp vectors
237 SmartPtr<vector_type> spD = u.clone_without_values();
238
239 std::string ext; //char ext[50];
240 int loopCnt = 0;
241 m_lastNumSteps = 0;
242 {
243 // write start defect for debug
244 //snprintf(ext, sizeof(ext), "_call%03d", m_dgbCall);
245 std::string name("NESTED_ITER_StartSolution");
246 ext = GetStringPrintf("_call%03d", m_dgbCall);
247 name.append(ext);
248 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 m_topLevel = surfLevel+m_maxSteps;
254 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 m_spAss->assemble_linear(*m_J, *spD, surfGridLevel); // todo: replace for non-linear problems
267 m_spAss->adjust_solution(u, surfGridLevel);
269
270 //snprintf(ext, sizeof(ext),"_call%03d_iter%03d", m_dgbCall, loopCnt);
271 ext = GetStringPrintf("_call%03d_iter%03d", m_dgbCall, loopCnt);
272 {
273 // write initial data for debug
274 std::string name0("NESTED_ITER_InitialSolution"); name0.append(ext);
275 write_debug(u, name0.c_str());
276 std::string name("NESTED_ITER_InitialDefect"); name.append(ext);
277 write_debug(*spD, name.c_str());
278
279 // Write Jacobian for debug
280 std::string matname("NESTED_ITER_Jacobian");
281 matname.append(ext);
282 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 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 }
295 }UG_CATCH_THROW("NestedIterationSolver::apply: Initialization of Linear Solver failed.");
296
297 // Solve linearized system.
298 try{
299 NESTED_ITER_PROFILE_BEGIN(NewtonApplyLinSolver);
300 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 }
307 }UG_CATCH_THROW("NestedIterationSolver::apply: Application of Linear Solver failed.");
308
309 // Adjust solution.
310 m_spAss->adjust_solution(u, surfGridLevel);
311 {
312 // write defect for debug
313 std::string name("NESTED_ITER_FinalDefect"); name.append(ext);
314 write_debug(*spD, name.c_str());
315 std::string name3("NESTED_ITER_FinalSolution"); name3.append(ext);
316 write_debug(u, name3.c_str());
317 }
318
319 // Update counter.
320 loopCnt++;
321
322 // Quit?
323 if (use_adaptive_refinement() == false) {
324 UG_LOG("SUCCESS: Non-adaptive version!" << std::endl);
325 break;
326 }
327
328 // Estimate and mark domain.
329 this->estimate_and_mark_domain(usol, m_spRefinementMarking);
330
331 // OPTIONAL: write defect for debug
332 if (m_spElemError.valid())
333 {
334 std::string name("NESTED_ITER_Error");
335 name.append(ext);
337 outError.print(name.c_str(), *m_spElemError);
338 }
339
340 // Check error.
341 const number err = m_spRefinementMarking->global_estimated_error();
342 double desiredTOL = (m_spAssociatedSpace.valid()) ? m_TOL*m_spAssociatedSpace->norm(usol) + m_absTOL : m_TOL;
343
344 UG_LOG("NI *** Desired tolerance: " << m_TOL << " * "<< m_spAssociatedSpace->norm(usol) << " + "<< m_absTOL <<std::endl);
345
346 // Quit or continue?
347 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 m_spRefiner->refine();
356 }
357 }
358
359
360 return true;
361}
362
363
364template <typename TDomain, typename TAlgebra>
366{
367 std::stringstream ss;
368 ss << "NestedIterationSolver\n";
369 ss << " LinearSolver: ";
370 if(m_spLinearSolver.valid()) ss << ConfigShift(m_spLinearSolver->config_string()) << "\n";
371 else ss << " NOT SET!\n";
372 return ss.str();
373}
374
375
376}
377
378#endif /* __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NESTED_ITER_SOLVER__NESTED_ITER_IMPL__ */
379
location name
Definition checkpoint_util.lua:128
Definition smart_pointer.h:108
bool valid() const
returns true if the pointer is valid, false if not.
Definition smart_pointer.h:206
matrix operator based on the assembling of a problem
Definition assembled_linear_operator.h:60
Definition assembled_non_linear_operator.h:43
domain discretization implementing the interface
Definition domain_disc.h:651
represents numerical solutions on a grid using an algebraic vector
Definition grid_function.h:121
const GridLevel & grid_level() const
returns the grid level
Definition grid_function.h:247
SmartPtr< DoFDistribution > dd()
Definition grid_function.h:237
Definition grid_level.h:42
bool top() const
returns if top level
Definition grid_level.h:83
bool is_surface() const
returns if type is surface
Definition grid_level.h:89
int level() const
returns the level
Definition grid_level.h:74
Interface providing Jacobian and Defect of a discretization.
Definition assemble_interface.h:110
Abstract base class for element marking (in adaptive refinement)
Definition error_elem_marking_strategy.h:125
describes an inverse linear mapping X->Y
Definition linear_operator_inverse.h:80
describes a mapping X->Y
Definition operator.h:86
linear solver using abstract preconditioner interface
Definition linear_solver.h:57
SmartPtr< AssembledOperator< algebra_type > > m_N
assembling routine
Definition nested_iteration.h:193
virtual std::string config_string() const
returns information about configuration parameters
Definition nested_iteration_impl.h:365
virtual bool apply(vector_type &u)
apply Operator, i.e. N^{-1}(0) = u
Definition nested_iteration_impl.h:208
SmartPtr< IAssemble< TAlgebra > > m_spAss
assembling
Definition nested_iteration.h:197
number coarsen_domain(const grid_function_type &u)
Coarsen all elements.
Definition nested_iteration_impl.h:192
virtual bool init(SmartPtr< IOperator< vector_type > > N)
Definition nested_iteration_impl.h:134
void estimate_and_mark_domain(const grid_function_type &u, SmartPtr< IElementMarkingStrategy< TDomain > > spMarking, bool bClearMarks=true)
Refines domain and provides error estimate.
Definition nested_iteration_impl.h:157
TAlgebra::vector_type vector_type
Vector type.
Definition nested_iteration.h:70
NestedIterationSolver()
default constructor
Definition nested_iteration_impl.h:71
virtual bool prepare(vector_type &u)
prepare Operator
Definition nested_iteration_impl.h:147
output writer to the VTK file format
Definition vtkoutput.h:150
void print(const char *filename, TFunction &u, int step, number time, bool makeConsistent)
Definition vtkoutput_impl.h:184
virtual void init()
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#define UG_CATCH_THROW(msg)
Definition error.h:64
#define UG_THROW(msg)
Definition error.h:57
#define UG_LOG(msg)
Definition log.h:367
double number
Definition types.h:124
the ug namespace
string ConfigShift(string s)
Definition string_util.cpp:457
std::string GetStringPrintf(const char *format, Args... args)
Definition string_util.h:449
#define NESTED_ITER_PROFILE_END()
Definition nested_iteration_impl.h:47
#define NESTED_ITER_PROFILE_BEGIN(name)
Definition nested_iteration_impl.h:46
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836