Line data Source code
1 : /*
2 : * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Arne Nägel
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__LINEAR_OPERATOR__DEBUG_ITERATOR__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__DEBUG_ITERATOR__
35 :
36 : #include "lib_algebra/operator/interface/preconditioner.h"
37 : #include "lib_algebra/operator/interface/preconditioned_linear_operator_inverse.h"
38 : #include "lib_algebra/algebra_common/sparsematrix_util.h"
39 : #include "lib_algebra/operator/debug_writer.h"
40 :
41 :
42 : #ifdef UG_PARALLEL
43 : #include "lib_algebra/parallelization/parallelization.h"
44 : #include "lib_algebra/parallelization/parallel_matrix_overlap_impl.h"
45 : #endif
46 :
47 : namespace ug{
48 :
49 : /// Debugging iterator
50 : /**
51 : * This class implements an iterator that can be used for debugging:
52 : *
53 : * It wraps for an ILinearIterator<typename TAlgebra::vector_type> object, which can be used as usual,
54 : * i.e., all call are forwarded. Moreover, this class can determine the algebarically smooth error. In order to doso, one must specify a solver.
55 : * This solver is called once during initialization.
56 : *
57 : * This class is a VectorDebugWritingObject <typename TAlgebra::vector_type>. If a DebugWriter is supplied using use the set_debug routines, the smooth error is written to file.
58 : *
59 : */
60 : template <typename TAlgebra>
61 : class DebugIterator :
62 : public virtual ILinearIterator<typename TAlgebra::vector_type>,
63 : public VectorDebugWritingObject <typename TAlgebra::vector_type>
64 : {
65 : public:
66 : // Algebra type
67 : typedef TAlgebra algebra_type;
68 :
69 : // Vector type
70 : typedef typename TAlgebra::vector_type vector_type;
71 :
72 : // Matrix type
73 : typedef typename TAlgebra::matrix_type matrix_type;
74 :
75 : // Base type
76 : typedef ILinearIterator<typename TAlgebra::vector_type> base_type;
77 :
78 : typedef IPreconditionedLinearOperatorInverse<vector_type> solver_type;
79 : private:
80 : typedef VectorDebugWritingObject <typename TAlgebra::vector_type> vdwo_type;
81 :
82 : public:
83 : /// Constructor
84 0 : DebugIterator() : from(0.0), to (1.0)
85 : {
86 0 : MatrixOperator<matrix_type, vector_type> *ptr =new MatrixOperator<matrix_type, vector_type>();
87 0 : m_pOperator = SmartPtr<MatrixOperator<matrix_type, vector_type> >(ptr);
88 0 : };
89 :
90 0 : DebugIterator(SmartPtr<base_type> pprecond, SmartPtr<solver_type> psolver) : from(0.0), to (1.0)
91 : {
92 0 : set_preconditioner(pprecond);
93 0 : set_solver(psolver);
94 0 : MatrixOperator<matrix_type, vector_type> *ptr =new MatrixOperator<matrix_type, vector_type>();
95 0 : m_pOperator = SmartPtr<MatrixOperator<matrix_type, vector_type> >(ptr);
96 0 : };
97 :
98 : /// Clone
99 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
100 : {
101 0 : SmartPtr<DebugIterator<algebra_type> > newInst(new DebugIterator<algebra_type> ());
102 0 : newInst->set_damp(this->damping());
103 0 : newInst->set_preconditioner(this->get_preconditioner()->clone()); // clone preconditioner
104 0 : newInst->set_solver(this->get_solver()); // keep solver (initialized below)
105 0 : newInst->set_debug(this->vdwo_type::vector_debug_writer());
106 0 : newInst->from = this->from;
107 0 : newInst->to = this->to;
108 0 : return newInst;
109 : }
110 :
111 : /// returns if parallel solving is supported
112 0 : virtual bool supports_parallel() const
113 : {
114 0 : if(m_pprecond.valid())
115 0 : return m_pprecond->supports_parallel();
116 : return true;
117 : }
118 :
119 : /// specify the real preconditioner (used for iterating)
120 0 : void set_preconditioner(SmartPtr<base_type> pprecond) {m_pprecond=pprecond;}
121 : /// specify the solver that will be used for debugging (optional)
122 0 : void set_solver(SmartPtr<solver_type> psolver) {m_solver=psolver;}
123 : /// specify bounds for random initial guess (optional)
124 0 : void set_random_bounds(double a, double b){from =a; to=b;}
125 :
126 0 : void set_solution(SmartPtr<vector_type> sol)
127 0 : {m_spSolVector = sol;}
128 :
129 : protected:
130 : /// get reference to 'real' preconditioner
131 : SmartPtr<base_type> get_preconditioner() {return m_pprecond;}
132 : /// get reference to aux. solver
133 : SmartPtr<solver_type> get_solver() {return m_solver;}
134 :
135 :
136 : /// name of preconditioner
137 0 : virtual const char* name() const
138 0 : {return "DebugIterator";}
139 :
140 :
141 : /// init (expensive)
142 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type> > J,
143 : const vector_type& u)
144 : {
145 : // cast to matrix based operator
146 0 : SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp =
147 : J.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
148 :
149 : // Check that matrix if of correct type
150 0 : if(pOp.invalid())
151 0 : UG_THROW(name() << "::init': Passed Operator is "
152 : "not based on matrix. This Preconditioner can only "
153 : "handle matrix-based operators.");
154 :
155 : // forward request to matrix based implementation
156 0 : return init(pOp);
157 : }
158 :
159 : /// init (expensive)
160 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type> > L)
161 : {
162 : // cast to matrix based operator
163 0 : SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp =
164 : L.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
165 :
166 : // Check that matrix if of correct type
167 0 : if(pOp.invalid())
168 0 : UG_THROW(name() << "::init': Passed Operator is "
169 : "not based on matrix. This Preconditioner can only "
170 : "handle matrix-based operators.");
171 :
172 : // forward request to matrix based implementation
173 0 : return init(pOp);
174 : }
175 :
176 : /// init (expensive, since it calls \sa find_smooth_error)
177 0 : bool init(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
178 : {
179 0 : m_pOperator->get_matrix() = pOp->get_matrix();
180 : //m_pOperator->get_matrix().set_as_copy_of(pOp->get_matrix());
181 : #ifdef UG_PARALLEL
182 : //m_pOperator->get_matrix().set_storage_Type(pOp->get_matrix()->get_storage_type());
183 : #endif
184 0 : m_pprecond->init(pOp);
185 :
186 : //write_debug(pOp->get_matrix(), "DebugMatrix");
187 0 : find_smooth_error();
188 :
189 0 : return true;
190 : }
191 :
192 : /// Determines algebraically smooth error
193 : /** Solves Ae=0 starting from a random initial guess.
194 : * \returns true, iff error was determined*/
195 0 : bool find_smooth_error()
196 : {
197 0 : if (m_solver.invalid())
198 : {
199 : UG_LOG("WARNING: cannot find smooth error; no solver supplied!");
200 0 : return false;
201 : }
202 :
203 0 : if (m_spSolVector.invalid())
204 : {
205 : UG_LOG("WARNING: cannot find smooth error; no valid vector !");
206 0 : return false;
207 : }
208 :
209 : //vector_type myRhs(m_pOperator->get_matrix().num_rows());
210 : //vector_type myError(m_pOperator->get_matrix().num_rows());
211 :
212 : SmartPtr<vector_type> myError = m_spSolVector;
213 0 : SmartPtr<vector_type> myRhs = m_spSolVector->clone_without_values();
214 :
215 : myError->set(0.0);
216 0 : myError->set_random(from, to);
217 :
218 : myRhs->set(0.0);
219 :
220 0 : this->write_debug(*myError, "DebugIterError0");
221 0 : m_solver->init(m_pOperator);
222 0 : m_solver->apply(*myError, *myRhs);
223 :
224 0 : this->write_debug(*myError, "DebugIterErrorS");
225 :
226 : return true;
227 : }
228 :
229 : /// forwarding call to original preconditioner
230 0 : virtual bool apply(vector_type& c, const vector_type& d)
231 : {
232 0 : return m_pprecond->apply(c, d);
233 : }
234 :
235 : /// forwarding call to original preconditioner
236 0 : virtual bool apply_update_defect(vector_type& c, vector_type& d)
237 : {
238 0 : return (m_pprecond->apply_update_defect(c, d));
239 : }
240 :
241 : protected:
242 :
243 : SmartPtr<base_type> m_pprecond;
244 :
245 : SmartPtr<solver_type> m_solver;
246 : SmartPtr<MatrixOperator<matrix_type, vector_type> > m_pOperator;
247 : SmartPtr<vector_type> m_spSolVector;
248 :
249 : double from, to;
250 :
251 : };
252 :
253 :
254 :
255 : } // end namespace ug
256 :
257 : #endif // __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__TRANSFORMING__
|
