Line data Source code
1 : /*
2 : * SPDX-FileCopyrightText: Copyright (c) 2014-2020: G-CSC, Goethe University Frankfurt
3 : * SPDX-License-Identifier: LicenseRef-UG4-LGPL-3.0
4 : *
5 : * Author: Arne Naegel
6 : *
7 : * This file is part of UG4.
8 : *
9 : * UG4 is free software: you can redistribute it and/or modify it under the
10 : * terms of the GNU Lesser General Public License version 3 (as published by the
11 : * Free Software Foundation) with the following additional attribution
12 : * requirements (according to LGPL/GPL v3 §7):
13 : *
14 : * (1) The following notice must be displayed in the Appropriate Legal Notices
15 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
16 : *
17 : * (2) The following notice must be displayed at a prominent place in the
18 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
19 : *
20 : * (3) The following bibliography is recommended for citation and must be
21 : * preserved in all covered files:
22 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
23 : * parallel geometric multigrid solver on hierarchically distributed grids.
24 : * Computing and visualization in science 16, 4 (2013), 151-164"
25 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
26 : * flexible software system for simulating pde based models on high performance
27 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
28 : *
29 : * This program is distributed in the hope that it will be useful,
30 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
31 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 : * GNU Lesser General Public License for more details.
33 : */
34 :
35 : #ifndef TIME_EXTRAPOLATION_H_
36 : #define TIME_EXTRAPOLATION_H_
37 :
38 : #include <vector>
39 : #include <cmath>
40 :
41 : // ug libraries
42 : #include "common/common.h"
43 : #include "common/util/smart_pointer.h"
44 :
45 : #include "lib_algebra/lib_algebra.h"
46 : #include "lib_algebra/operator/debug_writer.h"
47 :
48 : #include "lib_disc/function_spaces/grid_function.h"
49 : #include "lib_disc/function_spaces/grid_function_user_data.h"
50 : #include "lib_disc/function_spaces/integrate.h"
51 : #include "lib_disc/function_spaces/metric_spaces.h"
52 :
53 : #include "lib_disc/time_disc/time_disc_interface.h"
54 : #include "lib_disc/spatial_disc/user_data/linker/scale_add_linker.h"
55 :
56 :
57 :
58 : namespace ug{
59 :
60 : /*
61 : std::vector<size_t> steps {1, 2, 4};
62 : timex = new AitkenNevilleTimex(steps);
63 :
64 : //
65 : timex.set_solution(sol0, 0); // single step
66 : timex.set_solution(sol1, 1); // double step
67 : timex.set_solution(sol2, 2); // four step
68 :
69 : timex.set_estimator(L2NormEstimator())
70 : timex.set_estimator(InfNormEstimator())
71 : timex.set_estimator(RelativeEstimator())
72 : //
73 : timex.apply() // updating sol1 and sol2
74 :
75 : timex.get_error_estimate()
76 :
77 : */
78 :
79 : namespace tools {
80 :
81 :
82 : //! calculates vUpdate = vUpdate + alpha2*vFine + alpha3*vCoarse. for doubles
83 : inline void VecScaleAddWithNormRel(double &vUpdate, double alpha2, const double &vFine, double alpha3, const double &vCoarse, double &norm)
84 : {
85 0 : const double update = alpha2*vFine + alpha3*vCoarse;
86 0 : vUpdate = vUpdate + update;
87 0 : norm = std::max(norm, 0.5*fabs(update)/(1.0+fabs(vFine)+fabs(vCoarse)));
88 : }
89 :
90 : //! calculates vUpdate = vUpdate + alpha2*vFine + alpha3*vCoarse
91 : template<typename vector_t, template <class T> class TE_VEC>
92 0 : inline void VecScaleAddWithNormRel(TE_VEC<vector_t> &vUpdate, double alpha2, const TE_VEC<vector_t> &vFine, double alpha3, const TE_VEC<vector_t> &vCoarse, double &norm)
93 : {
94 0 : for(size_t i=0; i<vUpdate.size(); i++)
95 0 : VecScaleAddWithNormRel(vUpdate[i], alpha2, vFine[i], alpha3, vCoarse[i], norm);
96 0 : }
97 :
98 : //! calculates vUpdate = vUpdate + alpha2*vFine + alpha3*vCoarse. for doubles
99 : inline void VecScaleAddWithNormInf(double &vUpdate, double alpha2, const double &vFine, double alpha3, const double &vCoarse, double &norm)
100 : {
101 0 : const double update = alpha2*vFine + alpha3*vCoarse;
102 0 : vUpdate = vUpdate + update;
103 0 : norm = std::max(norm, fabs(update));
104 : }
105 :
106 : //! calculates vUpdate = vUpdate + alpha2*vFine + alpha3*vCoarse
107 : template<typename vector_t, template <class T> class TE_VEC>
108 0 : inline void VecScaleAddWithNormInf(TE_VEC<vector_t> &vUpdate, double alpha2, const TE_VEC<vector_t> &vFine, double alpha3, const TE_VEC<vector_t> &vCoarse, double &norm, const int delta=1, const int offset=0)
109 : {
110 : // std::cerr << norm << " "<< delta << offset << std::endl;
111 0 : for(size_t i=offset; i<vUpdate.size(); i+=delta)
112 : VecScaleAddWithNormInf(vUpdate[i], alpha2, vFine[i], alpha3, vCoarse[i], norm);
113 :
114 : // std::cerr << norm << std::endl;
115 0 : }
116 :
117 :
118 : //! calculates vUpdate = vUpdate + alpha2*vFine + alpha3*vCoarse. for doubles
119 : inline void VecScaleAddWithNorm2(double &vUpdate, double alpha2, const double &vFine, double alpha3, const double &vCoarse, double &norm)
120 : {
121 0 : const double update = alpha2*vFine + alpha3*vCoarse;
122 0 : vUpdate = vUpdate+ update;
123 0 : norm += update*update;
124 : }
125 :
126 : //! calculates vUpdate = vUpdate + alpha2*vFine + alpha3*vCoarse
127 : template<typename vector_t, template <class T> class TE_VEC>
128 0 : inline void VecScaleAddWithNorm2(TE_VEC<vector_t> &vUpdate, double alpha2, const TE_VEC<vector_t> &vFine, double alpha3, const TE_VEC<vector_t> &vCoarse, double &norm, const int delta=1, const int offset=0)
129 : {
130 0 : for(size_t i=offset; i<vUpdate.size(); i+=delta)
131 : VecScaleAddWithNorm2(vUpdate[i], alpha2, vFine[i], alpha3, vCoarse[i], norm);
132 0 : }
133 :
134 :
135 :
136 : }
137 :
138 : //! Interface for sub-diagonal error estimator (w.r.t time in Aitken-Neville scheme)
139 : /*! Given u_fine and u_coarse,
140 : *
141 : * */
142 : template <class TVector>
143 : class ISubDiagErrorEst
144 : {
145 : public:
146 : // constructor
147 0 : ISubDiagErrorEst() : m_est(0.0) {};
148 :
149 : // destructor
150 : virtual ~ISubDiagErrorEst() {};
151 :
152 : //! compute update vUpdate = vFine + alpha * (vFine- vCoarse) AND estimate error \| alpha * (vFine- vCoarse)\|
153 : virtual bool update(SmartPtr<TVector> vUpdate, number alpha2, SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse) = 0;
154 :
155 : //! get estimate
156 0 : number get_current_estimate() {return m_est; };
157 : void reset_estimate() { m_est=0.0; };
158 :
159 0 : virtual std::string config_string() const {return "ISubDiagErrorEst";}
160 :
161 : protected:
162 : number m_est;
163 : };
164 :
165 :
166 : //! Evaluate using (algebraic) infinity norm
167 : /*! Optional: provide component offset o and stride s, such that vector is evaluated at vec[o+k*s]
168 : * */
169 : template <class TVector>
170 : class NormInfEstimator : public ISubDiagErrorEst<TVector>
171 : {
172 : protected:
173 : typedef ISubDiagErrorEst<TVector> base_type;
174 : int m_stride;
175 : int m_offset;
176 :
177 : public:
178 : // constructor
179 0 : NormInfEstimator() :
180 0 : ISubDiagErrorEst<TVector>(), m_stride(1), m_offset(0) {};
181 :
182 :
183 : // apply
184 0 : bool update(SmartPtr<TVector> vUpdate, number alpha2, SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
185 : {
186 0 : const int delta = m_stride;
187 0 : const int offset = m_offset;
188 0 : base_type::m_est=0.0;
189 0 : tools::VecScaleAddWithNormInf(*vUpdate, alpha2, *vFine, -alpha2, *vCoarse, base_type::m_est, delta, offset);
190 0 : return true;
191 : }
192 :
193 : // offset (e.g., component to work on for systems w/ CPU1)
194 0 : void set_offset(int offset) {m_offset=offset;}
195 :
196 : // delta (e.g., total number components for systems w/ CPU1)
197 0 : void set_stride(int delta) {m_stride=delta;}
198 :
199 : };
200 :
201 :
202 : //! Evaluate using (algebraic) L2 norm
203 : template <class TVector>
204 : class Norm2Estimator : public ISubDiagErrorEst<TVector>
205 : {
206 : protected:
207 : typedef ISubDiagErrorEst<TVector> base_type;
208 : int m_stride;
209 : int m_offset;
210 :
211 : public:
212 : // constructor
213 0 : Norm2Estimator() :
214 0 : ISubDiagErrorEst<TVector>(), m_stride(1), m_offset(0) {};
215 : Norm2Estimator(int stride) :
216 : ISubDiagErrorEst<TVector>(), m_stride(stride), m_offset(0) {};
217 : Norm2Estimator(int delta, int offset) :
218 : ISubDiagErrorEst<TVector>(), m_stride(delta), m_offset (offset) {};
219 :
220 : // apply
221 0 : bool update(SmartPtr<TVector> vUpdate, number alpha2, SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
222 : {
223 0 : const int delta = m_stride;
224 0 : const int offset = m_offset;
225 0 : base_type::m_est=0.0;
226 0 : tools::VecScaleAddWithNorm2(*vUpdate, alpha2, *vFine, -alpha2, *vCoarse, base_type::m_est, delta, offset);
227 : #ifdef UG_PARALLEL
228 : double locEst = base_type::m_est;
229 : double globEst =0.0;
230 : vUpdate->layouts()->proc_comm().allreduce(&locEst, &globEst, 1, PCL_DT_DOUBLE, PCL_RO_SUM);
231 : base_type::m_est = globEst;
232 : #endif
233 0 : base_type::m_est = sqrt(base_type::m_est);
234 0 : return true;
235 : }
236 :
237 : // offset (e.g., component to work on for systems w/ CPU1)
238 0 : void set_offset(int offset) { m_offset=offset;}
239 :
240 : // delta (e.g., total number components for systems w/ CPU1)
241 0 : void set_stride(int delta) { m_stride=delta;}
242 :
243 : };
244 :
245 :
246 : /// Evaluate using (algebraic) L2 norm
247 : template <class TVector>
248 : class NormRelEstimator : public ISubDiagErrorEst<TVector>
249 : {
250 : protected:
251 : typedef ISubDiagErrorEst<TVector> base_type;
252 :
253 : public:
254 : // constructor
255 0 : NormRelEstimator() : ISubDiagErrorEst<TVector>() {};
256 :
257 : // apply w/ rel norm
258 0 : bool update(SmartPtr<TVector> vUpdate, number alpha2, SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
259 : {
260 0 : base_type::m_est=0;
261 0 : tools::VecScaleAddWithNormRel(*vUpdate, alpha2, *vFine, -alpha2, *vCoarse, base_type::m_est);
262 0 : return true;
263 : }
264 : };
265 :
266 : /*
267 : template<class TVector>
268 : class VectorDebugWritingEstimator
269 : : public VectorDebugWritingObject<TVector>
270 : {
271 : public:
272 : /// type of vector
273 : typedef TVector vector_type;
274 :
275 : public:
276 : VectorDebugWritingEstimator()
277 : : VectorDebugWritingObject<vector_type>() {}
278 :
279 : VectorDebugWritingEstimator(SmartPtr<IVectorDebugWriter<vector_type> > spDebugWriter)
280 : : VectorDebugWritingObject<vector_type>(spDebugWriter) {}
281 :
282 : int get_call_id() { return m_dgbCall; }
283 : void inc_call_id() { m_dgbCall++; }
284 :
285 : protected:
286 : /// call counter
287 : int m_dgbCall;
288 : };
289 : */
290 :
291 :
292 :
293 :
294 : /** Evaluates difference between two grid functions in L_inf norm */
295 : template <typename TGridFunction>
296 : class SupErrorEvaluator
297 : : public IComponentSpace<TGridFunction>
298 : {
299 : public:
300 : typedef IComponentSpace<TGridFunction> base_type;
301 :
302 0 : SupErrorEvaluator(const char *fctNames) : base_type(fctNames) {};
303 : //SupErrorEvaluator(const char *fctNames, number scale) : base_type(fctNames, 1, scale) {};
304 : SupErrorEvaluator(const char *fctNames, const char* ssNames /*, number scale*/)
305 0 : : base_type(fctNames, ssNames, 1/*, scale*/) {};
306 0 : ~SupErrorEvaluator() {};
307 :
308 :
309 : using IComponentSpace<TGridFunction>::norm;
310 : using IComponentSpace<TGridFunction>::distance;
311 :
312 : double norm(SmartPtr<TGridFunction> uFine)
313 : { return norm(*uFine); }
314 :
315 0 : double norm(TGridFunction& uFine)
316 : {
317 : // gather subsets in group
318 0 : SubsetGroup ssGrp(uFine.domain()->subset_handler());
319 0 : if (base_type::m_ssNames != NULL)
320 0 : ssGrp.add(TokenizeString(base_type::m_ssNames));
321 : else
322 0 : ssGrp.add_all();
323 :
324 0 : double maxVal = 0.0;
325 :
326 : // loop subsets
327 0 : for (size_t i = 0; i < ssGrp.size(); ++i)
328 : {
329 : // get subset index
330 : const int si = ssGrp[i];
331 :
332 : // loop elements of subset and dim
333 0 : switch (ssGrp.dim(i))
334 : {
335 : case DIM_SUBSET_EMPTY_GRID: break;
336 0 : case 0: maxVal = std::max(maxVal, findFctMaxOnSubset<Vertex>(uFine, si)); break;
337 0 : case 1: maxVal = std::max(maxVal, findFctMaxOnSubset<Edge>(uFine, si)); break;
338 0 : case 2: maxVal = std::max(maxVal, findFctMaxOnSubset<Face>(uFine, si)); break;
339 0 : case 3: maxVal = std::max(maxVal, findFctMaxOnSubset<Volume>(uFine, si)); break;
340 0 : default: UG_THROW("SupErrorEvaluator::norm: Dimension " << ssGrp.dim(i) << " not supported.");
341 : }
342 : }
343 :
344 : #ifdef UG_PARALLEL
345 : // max over processes
346 : if (pcl::NumProcs() > 1)
347 : {
348 : pcl::ProcessCommunicator com;
349 : number local = maxVal;
350 : com.allreduce(&local, &maxVal, 1, PCL_DT_DOUBLE, PCL_RO_MAX);
351 : }
352 : #endif
353 :
354 : // return the result
355 0 : return maxVal;
356 : }
357 :
358 0 : double norm2(TGridFunction& uFine)
359 0 : { double norm1 = norm(uFine); return norm1*norm1;}
360 :
361 0 : double distance(TGridFunction& uFine, TGridFunction& uCoarse)
362 : {
363 0 : UG_COND_THROW(uFine.dof_distribution().get() != uCoarse.dof_distribution().get(),
364 : "Coarse and fine solutions do not have the same underlying dof distro.");
365 :
366 0 : SmartPtr<TGridFunction> uErr = uCoarse.clone();
367 0 : uErr->operator-=(uFine);
368 0 : return norm(*uErr);
369 : }
370 :
371 0 : double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
372 0 : { double dist = distance(uFine, uCoarse); return dist*dist;}
373 :
374 : protected:
375 : template <typename TBaseElem>
376 0 : number findFctMaxOnSubset(const TGridFunction& u, int si) const
377 : {
378 : ConstSmartPtr<DoFDistribution> dd = u.dof_distribution();
379 :
380 : size_t fct;
381 : try {fct = dd->fct_id_by_name(base_type::m_fctNames.c_str());}
382 0 : UG_CATCH_THROW("Function index could not be determined.\n"
383 : "Bear in mind that only one function can be evaluated in this error evaluator.");
384 :
385 0 : number maxVal = 0.0;
386 0 : typename DoFDistribution::traits<TBaseElem>::const_iterator it = dd->begin<TBaseElem>(si);
387 0 : typename DoFDistribution::traits<TBaseElem>::const_iterator itEnd = dd->end<TBaseElem>(si);
388 0 : for (; it != itEnd; ++it)
389 : {
390 : std::vector<DoFIndex> vInd;
391 :
392 : // we compare against all indices on the element
393 : // which means most indices will be compared against quite often
394 : // but as this is a sup norm, this is not a problem (only in terms of performance)
395 0 : size_t nInd = dd->dof_indices(*it, fct, vInd, false, false);
396 0 : for (size_t i = 0; i < nInd; ++i)
397 0 : maxVal = std::max(maxVal, fabs(DoFRef(u, vInd[i])));
398 : }
399 :
400 0 : return maxVal;
401 : }
402 : };
403 :
404 :
405 :
406 : /*!
407 : * For arbitrary $\rho$, this class defines the integrand $|\rho(u_1)- \rho(u_2)|^2$.
408 : * If the grid function $u_2$ is not provided, the class returns $|\rho(u_1)|^2$
409 : * */
410 : template <typename TDataIn, typename TGridFunction>
411 0 : class DeltaSquareIntegrand
412 : : public StdIntegrand<number, TGridFunction::dim, DeltaSquareIntegrand<TDataIn, TGridFunction> >
413 : {
414 : public:
415 : // world dimension of grid function
416 : static const int worldDim = TGridFunction::dim;
417 :
418 : // data type
419 : typedef TDataIn data_type;
420 :
421 : private:
422 :
423 : static number product(const number &x, const number &y)
424 0 : {return x*y;}
425 :
426 : static number product(const MathVector<worldDim> &x, const MathVector<worldDim> &y)
427 : {return VecDot(x,y);}
428 :
429 :
430 : // data to integrate
431 : SmartPtr<UserData<TDataIn, worldDim> > m_spData;
432 :
433 : // grid function
434 : const TGridFunction* m_pGridFct1;
435 : const TGridFunction* m_pGridFct2;
436 :
437 : // time
438 : number m_time;
439 :
440 : public:
441 : /// constructor
442 0 : DeltaSquareIntegrand(SmartPtr<UserData<TDataIn, worldDim> > spData,
443 : const TGridFunction* pGridFct1, const TGridFunction* pGridFct2,
444 : number time)
445 0 : : m_spData(spData), m_pGridFct1(pGridFct1), m_pGridFct2(pGridFct2), m_time(time)
446 : {
447 0 : m_spData->set_function_pattern(pGridFct1->function_pattern());
448 0 : };
449 :
450 : /// constructor
451 : DeltaSquareIntegrand(SmartPtr<UserData<TDataIn, worldDim> > spData, number time)
452 : : m_spData(spData), m_pGridFct1(NULL), m_pGridFct2(NULL), m_time(time)
453 : {
454 : if(m_spData->requires_grid_fct())
455 : UG_THROW("UserDataDeltaIntegrand: Missing GridFunction, but "
456 : " data requires grid function.")
457 : };
458 :
459 :
460 :
461 :
462 : template <int elemDim>
463 0 : void get_values(TDataIn vValue[],
464 : ConstSmartPtr<UserData<TDataIn, worldDim> > spData,
465 : const TGridFunction& gridFct,
466 : const MathVector<worldDim> vGlobIP[],
467 : GridObject* pElem,
468 : const MathVector<worldDim> vCornerCoords[],
469 : const MathVector<elemDim> vLocIP[],
470 : const MathMatrix<elemDim, worldDim> vJT[],
471 : const size_t numIP)
472 : {
473 : // collect local solution, if required
474 0 : if(spData->requires_grid_fct())
475 : {
476 : // create storage
477 : LocalIndices ind;
478 : LocalVector u;
479 :
480 : // get global indices
481 : gridFct.indices(pElem, ind);
482 :
483 : // adapt local algebra
484 0 : u.resize(ind);
485 :
486 : // read local values of u
487 0 : GetLocalVector(u, gridFct);
488 0 : std::cout << u << std::endl;
489 :
490 : // compute data
491 : try{
492 0 : (*spData)(vValue, vGlobIP, m_time, this->m_si, pElem,
493 : vCornerCoords, vLocIP, numIP, &u, &vJT[0]);
494 : }
495 0 : UG_CATCH_THROW("UserDataDeltaIntegrand: Cannot evaluate data.");
496 : }
497 : else
498 : {
499 : // compute data
500 : try{
501 0 : (*spData)(vValue, vGlobIP, m_time, this->m_si, numIP);
502 : }
503 0 : UG_CATCH_THROW("UserDataDeltaIntegrand: Cannot evaluate data.");
504 : }
505 :
506 0 : };
507 :
508 :
509 : /// \copydoc IIntegrand::values
510 : template <int elemDim>
511 0 : void evaluate(number vValue[],
512 : const MathVector<worldDim> vGlobIP[],
513 : GridObject* pElem,
514 : const MathVector<worldDim> vCornerCoords[],
515 : const MathVector<elemDim> vLocIP[],
516 : const MathMatrix<elemDim, worldDim> vJT[],
517 : const size_t numIP)
518 : {
519 0 : std::vector<TDataIn> v1(numIP);
520 :
521 0 : get_values<elemDim>(&v1[0], m_spData, *m_pGridFct1, vGlobIP, pElem, vCornerCoords, vLocIP, vJT, numIP);
522 : std::cout << "--- got v1!" << std::endl;
523 :
524 0 : if (m_pGridFct2 != NULL)
525 : {
526 0 : std::vector<TDataIn> v2(numIP);
527 : /* m_spGridFct->set(0.5);
528 : m_spGridFct2->set(0.5);*/
529 0 : get_values<elemDim>(&v2[0], m_spData, *m_pGridFct2, vGlobIP, pElem, vCornerCoords, vLocIP, vJT, numIP);
530 : std::cout << "--- got v2!" << std::endl;
531 :
532 0 : for (size_t ip=0; ip<numIP; ++ip)
533 : {
534 :
535 0 : std::cout << std::setprecision(12) << v1[ip] <<" "<< std::setprecision(12) << v2[ip] << std::endl;
536 0 : v1[ip] -= v2[ip];
537 0 : vValue[ip] = this->product(v1[ip], v1[ip]);
538 :
539 : }
540 0 : } else {
541 :
542 0 : for (size_t ip=0; ip<numIP; ++ip)
543 0 : { vValue[ip] = this->product(v1[ip], v1[ip]); }
544 :
545 : }
546 0 : };
547 : };
548 :
549 : //! Evaluate the difference for a (dependent) UserData type induced by different grid functions
550 : /*! UserData maybe of type TDataInput, i.e., number/vector/matrix/... */
551 : template <typename TGridFunction, typename TDataInput>
552 : class UserDataSpace : public IComponentSpace<TGridFunction>
553 : {
554 : public:
555 : typedef IComponentSpace<TGridFunction> base_type;
556 : typedef UserData<TDataInput, TGridFunction::dim> input_user_data_type;
557 :
558 0 : UserDataSpace(const char *fctNames) : base_type(fctNames) {};
559 0 : UserDataSpace(const char *fctNames, int order) : base_type(fctNames, order) {};
560 0 : ~UserDataSpace() {};
561 :
562 0 : void set_user_data(SmartPtr<input_user_data_type> spData)
563 0 : { m_userData = spData; }
564 :
565 :
566 : //! per default
567 : using IComponentSpace<TGridFunction>::norm;
568 : using IComponentSpace<TGridFunction>::distance;
569 :
570 0 : double norm2(TGridFunction& uFine)
571 : {
572 0 : DeltaSquareIntegrand<TDataInput, TGridFunction> spIntegrand(m_userData, &uFine, NULL, 0.0);
573 0 : return IntegrateSubsets(spIntegrand, uFine, base_type::m_ssNames, base_type::m_quadorder, "best");
574 : }
575 :
576 0 : double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
577 : {
578 0 : DeltaSquareIntegrand<TDataInput, TGridFunction> spIntegrand(m_userData, &uFine, &uCoarse, 0.0);
579 : std::cerr << "uFine="<<(void*) (&uFine) << ", uCoarse="<< (void*) (&uCoarse) << std::endl;
580 0 : return IntegrateSubsets(spIntegrand, uFine, base_type::m_ssNames, base_type::m_quadorder, "best");
581 : }
582 :
583 : protected:
584 : SmartPtr<input_user_data_type> m_userData;
585 : };
586 :
587 :
588 : /*
589 : /// Evaluate difference between two functions (w.r.t various norms)
590 : template <class TDomain, class TAlgebra>
591 : class GridFunctionEstimator :
592 : public ISubDiagErrorEst<typename TAlgebra::vector_type>
593 : {
594 : protected:
595 : typedef typename TAlgebra::vector_type TVector;
596 : typedef GridFunction<TDomain, TAlgebra> grid_function_type;
597 : typedef IErrorEvaluator<grid_function_type> evaluator_type;
598 :
599 : std::vector<SmartPtr<evaluator_type> > m_evaluators;
600 : number m_refNormValue;
601 :
602 : public:
603 : typedef ISubDiagErrorEst<TVector> base_type;
604 :
605 : // constructor
606 : ScaledGridFunctionEstimator() : m_refNormValue(0.0) {}
607 : ScaledGridFunctionEstimator(number ref) : m_refNormValue(ref) {}
608 :
609 : void add(SmartPtr<evaluator_type> eval)
610 : {
611 : m_evaluators.push_back(eval);
612 : }
613 :
614 : // apply w/ rel norm
615 : bool update(SmartPtr<TVector> vUpdate, number alpha, SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
616 : {
617 : // typedef ScaleAddLinker<number, TDomain::dim, number> linker_type;
618 : typedef GridFunctionNumberData<TGridFunction> TNumberData;
619 :
620 : // try upcast
621 : SmartPtr<grid_function_type> uFine = vFine.template cast_dynamic<grid_function_type>();
622 : SmartPtr<grid_function_type> uCoarse = vCoarse.template cast_dynamic<grid_function_type>();
623 : if (uFine.invalid() || uCoarse.invalid()) return false;
624 :
625 : // error estimate
626 : if (m_refNormValue<=0.0)
627 : {
628 : // relative error estimator
629 : //number unorm = L2Norm(uFine, m_fctNames.c_str(), m_quadorder);
630 : //number enorm = alpha*L2Error(uFine, m_fctNames.c_str(), uCoarse, m_fctNames.c_str() ,m_quadorder);
631 : number unorm = 0.0;
632 : number enorm = 0.0;
633 : double est = 0.0;
634 : for (typename std::vector<evaluator_type>::iterator it = m_evaluators.begin(); it!= m_evaluators.end(); ++it)
635 : {
636 : enorm = alpha * it->distance(uFine, uCoarse);
637 : unorm = std::max(it->norm(uFine), 1e-10);
638 : est += (enorm*enorm)/(unorm*unorm);
639 : std::cerr << "unorm=" << unorm << "enorm=" << enorm << "est="<<est << std::endl;
640 : }
641 :
642 : base_type::m_est = sqrt(est)/m_evaluators.size();
643 : std::cerr << "eps="<< base_type::m_est << std::endl;
644 :
645 : }
646 : else
647 : {
648 : // weighted error estimator
649 : number enorm = 0.0;
650 : for (typename std::vector<evaluator_type>::iterator it = m_evaluators.begin(); it!= m_evaluators.end(); ++it)
651 : {
652 : enorm += alpha * it->distance(uFine, uCoarse);
653 : }
654 : base_type::m_est = enorm/m_refNormValue;
655 :
656 : std::cerr << "unorm (FIXED)=" << m_refNormValue << "enorm=" << enorm << "eps="<< base_type::m_est << std::endl;
657 :
658 : }
659 :
660 : // update
661 : VecScaleAdd(*vUpdate, 1.0+alpha, *vFine, -alpha, *vCoarse);
662 : return true;
663 : }
664 :
665 : void set_reference_norm(number norm)
666 : {m_refNormValue = norm; }
667 :
668 :
669 : /// print config string
670 : std::string config_string() const
671 : {
672 : std::stringstream ss;
673 : ss << "GridFunctionEstimator:\n";
674 : for (typename std::vector<GridFunctionEvaluator>::const_iterator it = m_evaluators.begin(); it!= m_evaluators.end(); ++it)
675 : {
676 : ss << it->config_string();
677 : }
678 : return ss.str();
679 : }
680 : };
681 : */
682 :
683 : //! Evaluate using continuous norms
684 : /*! Can provide subspaces (incl. scaling).
685 : *
686 : * sum_I \sigma_I \frac{\|e\|_I}{\|u\|_I}
687 : * For ref value <= 0 TOL is scaled relatively w.r.t. function norm.
688 : */
689 : template <class TDomain, class TAlgebra>
690 : class GridFunctionEstimator :
691 : public ISubDiagErrorEst<typename TAlgebra::vector_type>
692 : {
693 : protected:
694 : typedef typename TAlgebra::vector_type TVector;
695 : public:
696 : typedef ISubDiagErrorEst<TVector> base_type;
697 : typedef GridFunction<TDomain, TAlgebra> grid_function_type;
698 : typedef IGridFunctionSpace<grid_function_type> subspace_type ;
699 : //typedef CompositeSpace<grid_function_type> composite_type;
700 :
701 : protected:
702 : typedef std::pair<SmartPtr<subspace_type>, number> weighted_obj_type;
703 :
704 : number m_refNormValue;
705 : std::vector<weighted_obj_type> m_spWeightedSubspaces;
706 :
707 :
708 : public:
709 0 : GridFunctionEstimator() : ISubDiagErrorEst<TVector>(), m_refNormValue(0.0)
710 : {}
711 :
712 :
713 0 : GridFunctionEstimator(double ref) : ISubDiagErrorEst<TVector>(), m_refNormValue(ref)
714 : {}
715 :
716 : /// add sub-space component
717 0 : void add(SmartPtr<subspace_type> spSubspace)
718 0 : { m_spWeightedSubspaces.push_back(std::make_pair(spSubspace, 1.0)); }
719 :
720 0 : void add(SmartPtr<subspace_type> spSubspace, number sigma)
721 0 : { m_spWeightedSubspaces.push_back(std::make_pair(spSubspace, sigma)); }
722 :
723 :
724 : /// apply w/ rel norm
725 0 : bool update(SmartPtr<TVector> vUpdate, number alpha, SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
726 : {
727 : // typedefs
728 : typedef ScaleAddLinker<number, TDomain::dim, number> linker_type;
729 : typedef GridFunction<TDomain, TAlgebra> grid_function_type;
730 :
731 : // upcast to GridFunction
732 0 : SmartPtr<grid_function_type> uFine = vFine.template cast_dynamic<grid_function_type>();
733 0 : SmartPtr<grid_function_type> uCoarse = vCoarse.template cast_dynamic<grid_function_type>();
734 0 : if (uFine.invalid() || uCoarse.invalid()) return false;
735 :
736 : // error estimate
737 0 : if (m_refNormValue<=0.0)
738 : {
739 : // relative error estimator
740 : number unorm2 = 0.0;
741 : number enorm2 = 0.0;
742 0 : for (typename std::vector<weighted_obj_type>::iterator it = m_spWeightedSubspaces.begin(); it!= m_spWeightedSubspaces.end(); ++it)
743 : {
744 0 : const double sigma = it->second;
745 0 : const double norm2 = it->first->norm2(*uFine);
746 0 : const double dist2 = alpha * alpha * (it->first->distance2(*uFine, *uCoarse));
747 0 : unorm2 += sigma * norm2;
748 0 : enorm2 += sigma * dist2;
749 :
750 : UG_LOG("unorm=" << norm2 << "\tenorm=" << dist2 << "\tsigma=" << sigma << std::endl);
751 : }
752 :
753 0 : base_type::m_est = sqrt(enorm2/unorm2);
754 0 : UG_LOG(">>> unorm2=" << unorm2 << "\tenorm2=" << enorm2 << "\teps="<< base_type::m_est << std::endl);
755 : }
756 : else
757 : {
758 : // weighted error estimator
759 : number enorm2 = 0.0;
760 0 : for (typename std::vector<weighted_obj_type>::iterator it = m_spWeightedSubspaces.begin(); it!= m_spWeightedSubspaces.end(); ++it)
761 : {
762 0 : const double sigma = it->second;
763 0 : const double dist2 = alpha * alpha * (it->first->distance2(*uFine, *uCoarse));
764 0 : enorm2 += sigma*dist2;
765 : }
766 :
767 0 : base_type::m_est = sqrt(enorm2)/m_refNormValue;
768 0 : UG_LOG("unorm (FIXED)=" << m_refNormValue << "enorm2=" << enorm2 << "eps="<< base_type::m_est << std::endl);
769 : }
770 :
771 : // update
772 0 : VecScaleAdd(*vUpdate, 1.0+alpha, *vFine, -alpha, *vCoarse);
773 0 : return true;
774 : }
775 :
776 0 : void set_reference_norm(number norm)
777 0 : {m_refNormValue = norm; }
778 :
779 :
780 : /// print config string
781 0 : std::string config_string() const
782 : {
783 0 : std::stringstream ss;
784 0 : ss << "GridFunctionEstimator:\n";
785 0 : for (typename std::vector<weighted_obj_type>::const_iterator it = m_spWeightedSubspaces.begin(); it!= m_spWeightedSubspaces.end(); ++it)
786 : {
787 0 : ss << it->second << "*" << it->first->config_string();
788 : }
789 0 : return ss.str();
790 0 : }
791 : };
792 :
793 : /// Evaluate difference between two functions (w.r.t various norms)
794 : template <class TDomain, class TAlgebra>
795 : class ScaledGridFunctionEstimator :
796 : public ISubDiagErrorEst<typename TAlgebra::vector_type>
797 : {
798 : protected:
799 : typedef typename TAlgebra::vector_type TVector;
800 :
801 : public:
802 : typedef ISubDiagErrorEst<TVector> base_type;
803 : typedef GridFunction<TDomain, TAlgebra> grid_function_type;
804 : typedef IComponentSpace<grid_function_type> subspace_type;
805 : typedef CompositeSpace<grid_function_type> composite_type;
806 :
807 : // constructor
808 0 : ScaledGridFunctionEstimator() : base_type() {}
809 :
810 : // add (single subspace)
811 0 : void add(SmartPtr<subspace_type> spSubspace)
812 0 : { m_spSubspaces.push_back(spSubspace); }
813 :
814 : // add subspaces (from container)
815 0 : void add(SmartPtr<composite_type> spCompositeSpace)
816 : {
817 : typedef typename composite_type::weighted_obj_type weighted_obj_type;
818 : const std::vector<weighted_obj_type> &spaces = spCompositeSpace->get_subspaces();
819 0 : for (typename std::vector<weighted_obj_type>::const_iterator it = spaces.begin(); it != spaces.end(); ++it)
820 : {
821 0 : m_spSubspaces.push_back(it->first);
822 : }
823 0 : }
824 :
825 :
826 : // apply w/ rel norm
827 0 : bool update(SmartPtr<TVector> vUpdate, number alpha, SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
828 : {
829 :
830 : // upcast
831 0 : SmartPtr<grid_function_type> uFine = vFine.template cast_dynamic<grid_function_type>();
832 0 : SmartPtr<grid_function_type> uCoarse = vCoarse.template cast_dynamic<grid_function_type>();
833 0 : if (uFine.invalid() || uCoarse.invalid()) return false;
834 :
835 : // error estimate
836 : number est = 0.0;
837 0 : for (typename std::vector<SmartPtr<subspace_type> >::iterator it = m_spSubspaces.begin();
838 0 : it!= m_spSubspaces.end(); ++it)
839 : {
840 : // use sub-diagonal error estimator (i.e. multiply with alpha)
841 0 : double enorm2 = (alpha*alpha) * (*it)->distance2(*uFine, *uCoarse);
842 0 : double unorm2 = std::max((*it)->norm2(*uFine), 1e-10*enorm2);
843 0 : est += (enorm2)/(unorm2);
844 : UG_LOGN("unorm2=" << unorm2 << "\tenorm2=" << enorm2 << "\tratio2="<< (enorm2)/(unorm2) << "est2=" << est);
845 : }
846 :
847 0 : base_type::m_est = sqrt(est)/m_spSubspaces.size();
848 0 : UG_LOGN("eps="<< base_type::m_est);
849 :
850 : // update
851 0 : VecScaleAdd(*vUpdate, 1.0+alpha, *vFine, -alpha, *vCoarse);
852 0 : return true;
853 : }
854 :
855 :
856 : /// print config string
857 0 : std::string config_string() const
858 : {
859 0 : std::stringstream ss;
860 0 : ss << "ScaledGridFunctionEstimator:\n";
861 0 : for (typename std::vector<SmartPtr<subspace_type> >::const_iterator it = m_spSubspaces.begin();
862 0 : it!= m_spSubspaces.end(); ++it)
863 : {
864 0 : ss << (*it)->config_string();
865 : }
866 0 : return ss.str();
867 0 : }
868 :
869 : protected:
870 :
871 : std::vector<SmartPtr<subspace_type> > m_spSubspaces;
872 :
873 : };
874 :
875 :
876 : //! Evaluate difference between two functions (w.r.t various norms)
877 : /*! \| u \|^2 = \sum |u|_i^2 */
878 : template <class TDomain, class TAlgebra>
879 : class CompositeGridFunctionEstimator :
880 : public ISubDiagErrorEst<typename TAlgebra::vector_type>
881 : {
882 : protected:
883 : typedef typename TAlgebra::vector_type TVector;
884 :
885 : public:
886 : typedef ISubDiagErrorEst<TVector> base_type;
887 : typedef GridFunction<TDomain, TAlgebra> grid_function_type;
888 : typedef IComponentSpace<grid_function_type> subspace_type;
889 : typedef CompositeSpace<grid_function_type> composite_type;
890 :
891 : // constructor
892 0 : CompositeGridFunctionEstimator() : base_type(), m_strictRelativeError(0) {}
893 :
894 0 : void use_strict_relative_norms(int b)//bool b)
895 0 : { m_strictRelativeError = b; }
896 :
897 0 : void use_strict_relative_norms(bool b){
898 0 : if(b){
899 0 : m_strictRelativeError=1;
900 : }else{
901 0 : m_strictRelativeError=0;
902 : }
903 0 : }
904 :
905 : // add (single subspace)
906 0 : void add(SmartPtr<subspace_type> spSubspace)
907 0 : { m_spSubspaces.push_back(spSubspace); }
908 :
909 : // add subspaces (from container)
910 0 : void add(SmartPtr<composite_type> spCompositeSpace)
911 : {
912 : typedef typename composite_type::weighted_obj_type weighted_obj_type;
913 : const std::vector<weighted_obj_type> &spaces = spCompositeSpace->get_subspaces();
914 0 : for (typename std::vector<weighted_obj_type>::const_iterator it = spaces.begin(); it != spaces.end(); ++it)
915 : {
916 0 : m_spSubspaces.push_back(it->first);
917 : }
918 0 : }
919 :
920 :
921 : // apply w/ rel norm
922 0 : bool update(SmartPtr<TVector> vUpdate, number alpha, SmartPtr<TVector> vFine, SmartPtr<TVector> vCoarse)
923 : {
924 :
925 : // upcast
926 0 : SmartPtr<grid_function_type> uFine = vFine.template cast_dynamic<grid_function_type>();
927 0 : SmartPtr<grid_function_type> uCoarse = vCoarse.template cast_dynamic<grid_function_type>();
928 0 : if (uFine.invalid() || uCoarse.invalid()) return false;
929 :
930 : size_t numFct=uFine->num_fct();
931 : std::vector<double> vcmp_e2, vcmp_u2;
932 0 : vcmp_e2.resize(numFct, 0);
933 0 : vcmp_u2.resize(numFct, 0);
934 :
935 : // error estimate
936 : double enorm2 = 0.0;
937 0 : double unorm2 = 0.0;
938 : const double SMALL = 1e-10;
939 :
940 0 : double cmp_rel=0.0;
941 0 : double cmp_rel4=0.0;
942 0 : double max_rel = 0.0;
943 0 : double max_rel_e4=0.0;
944 :
945 0 : for (typename std::vector<SmartPtr<subspace_type> >::iterator it = m_spSubspaces.begin();
946 0 : it!= m_spSubspaces.end(); ++it)
947 : {
948 : // use sub-diagonal error estimator (i.e. multiply with alpha)
949 : std::string cmp=(*it)->function_name();
950 : size_t cmp_id= uFine->fct_id_by_name(cmp.c_str());
951 0 : double cmp_e2 = (alpha*alpha) * (*it)->distance2(*uFine, *uCoarse);
952 0 : double cmp_u2 = (*it)->norm2(*uFine);
953 0 : vcmp_e2[cmp_id]+=cmp_e2;
954 0 : vcmp_u2[cmp_id]+=cmp_u2;
955 0 : cmp_rel=std::max(cmp_rel,(cmp_e2==0)?0:cmp_e2/(cmp_u2 + SMALL/ (1.0+cmp_e2+cmp_u2)));
956 0 : cmp_rel4=std::max(cmp_rel4,(cmp_e2==0)?0:cmp_e2/std::sqrt(cmp_u2*cmp_u2+SMALL*SMALL));
957 0 : UG_LOGN("cmp=" << cmp << "(id=" << cmp_id << "):\tui-2=" << cmp_u2 << "\tei-2=" << cmp_e2<<
958 : "\tratio2="<< (vcmp_e2[cmp_id])/(vcmp_u2[cmp_id]) <<
959 : "\tmax.rel(squared)="<< cmp_rel<<"\tmax.rel(ei4)="<<cmp_rel4
960 : );
961 : }
962 :
963 0 : for(size_t cmp_id=0; cmp_id<numFct; ++cmp_id){
964 0 : double cmp_e2=vcmp_e2[cmp_id];
965 0 : double cmp_u2=vcmp_u2[cmp_id];
966 : // |delta_i|/|u_i|
967 0 : max_rel = std::max(max_rel, (cmp_e2==0)?0:cmp_e2/
968 0 : (cmp_u2 + SMALL/ (1.0+cmp_e2+cmp_u2)));
969 : //std::max(cmp_u2, SMALL*cmp_e2));
970 0 : max_rel_e4=std::max(max_rel_e4, (cmp_e2==0)?0:cmp_e2/std::sqrt(cmp_u2*cmp_u2+SMALL*SMALL));
971 :
972 0 : enorm2 += cmp_e2;
973 0 : unorm2 += cmp_u2;
974 :
975 0 : UG_LOGN("summe by cmp_id=" << cmp_id << ":\tui-2=" << cmp_u2 << "\tei-2=" << cmp_e2<<
976 : "\tratio2="<< (enorm2)/(unorm2) <<
977 : "\tmax.rel(squared)="<< max_rel<<"\tmax.rel(ei4)="<<max_rel_e4);
978 : }
979 :
980 : //prevent division by zero
981 : /*base_type::m_est = (m_strictRelativeError) ? sqrt(max_rel) :
982 : sqrt(enorm2/std::max(unorm2, SMALL*enorm2));*/
983 :
984 0 : switch(m_strictRelativeError)
985 : {
986 0 : case 0:
987 0 : base_type::m_est=sqrt(enorm2/std::max(unorm2, SMALL*enorm2)); /* absolute: */
988 0 : break;
989 0 : case 1:
990 0 : base_type::m_est=sqrt(cmp_rel); /* relative max: max(e2/(u2+SMALL/(1+e2+u2))) for each cmp on each subset */
991 0 : break;
992 0 : case 2:
993 0 : base_type::m_est=sqrt(cmp_rel4); /* relative max: max(e2/sqrt(u2*u2+SMALL*SMALL)) for each cmp on each subset*/
994 0 : break;
995 0 : case 3:
996 0 : base_type::m_est=sqrt(max_rel); /* relative max based on cmp: max(sum_si e2/( sum_si u2+SMALL/(1+sum_si e2+ sum_si u2))) for each cmp on whole domain.*/
997 0 : break;
998 0 : case 4:
999 0 : base_type::m_est=sqrt(max_rel_e4); /* relative max based on cmp: max(sum_si e2/sqrt( sum_si u2* sum_si u2+SMALL*SMALL)) for each cmp on whole domain.*/
1000 0 : break;
1001 : }
1002 0 : UG_LOGN("eps="<< base_type::m_est<<" with strictRelativeError="<<m_strictRelativeError);
1003 :
1004 : // update
1005 0 : VecScaleAdd(*vUpdate, 1.0+alpha, *vFine, -alpha, *vCoarse);
1006 : return true;
1007 0 : }
1008 :
1009 :
1010 : /// print config string
1011 0 : std::string config_string() const
1012 : {
1013 0 : std::stringstream ss;
1014 0 : ss << "CompositeGridFunctionEstimator:\n";
1015 0 : for (typename std::vector<SmartPtr<subspace_type> >::const_iterator it = m_spSubspaces.begin();
1016 0 : it!= m_spSubspaces.end(); ++it)
1017 : {
1018 0 : ss << (*it)->config_string();
1019 : }
1020 0 : return ss.str();
1021 0 : }
1022 :
1023 : protected:
1024 :
1025 : std::vector<SmartPtr<subspace_type> > m_spSubspaces;
1026 : //bool m_strictRelativeError;
1027 : int m_strictRelativeError;
1028 : };
1029 :
1030 :
1031 :
1032 : template <typename TVector>
1033 : class AitkenNevilleTimex
1034 : {
1035 : public:
1036 : /// vector type of solutions
1037 : typedef TVector vector_type;
1038 :
1039 : public:
1040 :
1041 : /** Aitken Neville scheme with a given number */
1042 0 : AitkenNevilleTimex(std::vector<size_t> nsteps)
1043 0 : : m_stepsize(0.0),
1044 0 : m_subdiag (make_sp(new Norm2Estimator<TVector>())),
1045 0 : m_num_steps(nsteps),
1046 0 : m_solution(nsteps.size()),
1047 0 : m_subdiag_error_est(nsteps.size(), INFINITY)
1048 0 : {};
1049 :
1050 0 : AitkenNevilleTimex(std::vector<size_t> nsteps, SmartPtr<ISubDiagErrorEst<vector_type> > error)
1051 0 : : m_stepsize(0.0),
1052 : m_subdiag(error),
1053 0 : m_num_steps(nsteps),
1054 0 : m_solution(nsteps.size()),
1055 0 : m_subdiag_error_est(nsteps.size(), INFINITY)
1056 0 : {};
1057 :
1058 0 : virtual ~AitkenNevilleTimex() {}
1059 :
1060 0 : void set_global_stepsize(number H) {m_stepsize=H;}
1061 0 : number get_global_stepsize() {return m_stepsize;}
1062 :
1063 : //! set solution (for stage i)
1064 0 : void set_solution(SmartPtr<vector_type> soli, int i)
1065 0 : { m_solution[i] = soli; }
1066 :
1067 : //! get solution (on stage i)
1068 0 : SmartPtr<vector_type> get_solution(size_t i)
1069 0 : { return m_solution[i]; }
1070 :
1071 : //! set error estimator
1072 0 : void set_error_estimate(SmartPtr<ISubDiagErrorEst<vector_type> > subdiag)
1073 : {
1074 0 : m_subdiag = subdiag;
1075 0 : }
1076 :
1077 :
1078 0 : const std::vector<number>& get_error_estimates() const
1079 0 : {return m_subdiag_error_est; }
1080 :
1081 : //! error estimate on stage k
1082 0 : number get_error_estimate(int k) const
1083 0 : { return m_subdiag_error_est[k];}
1084 :
1085 :
1086 : //! best error estimate
1087 : /*number get_error_estimate()
1088 : {
1089 : std::vector<number>::iterator best = std::min_element(m_subdiag_error_est.begin(), m_subdiag_error_est.end());
1090 : return *best;
1091 : }
1092 :
1093 : //! return best index
1094 : int get_best_index() const
1095 : {
1096 : std::vector<number>::iterator best = std::min_element(m_subdiag_error_est.begin(), m_subdiag_error_est.end());
1097 : //std::cout << "min element at: " << std::distance(std::begin(m_subdiag_error_est), best);
1098 : return std::distance(m_subdiag_error_est.begin(), best);
1099 : }
1100 :
1101 : */
1102 :
1103 : /**
1104 : * Triangular Aitken Neville extrapolation:
1105 : *
1106 : * T_11
1107 : * T_21 T_22
1108 : *
1109 : * T_N1 T_NN
1110 : *
1111 : * for T_ik
1112 : * */
1113 0 : void apply(size_t nstages, bool with_error=true)
1114 : {
1115 : UG_ASSERT(nstages <= m_solution.size(),
1116 : "Dimensions do not match:" << nstages << ">" << m_solution.size());
1117 :
1118 0 : if (with_error)
1119 : {
1120 : // reset (for safety reasons...)
1121 0 : for (size_t k=1; k<m_solution.size(); ++k)
1122 0 : { m_subdiag_error_est[k] = INFINITY; }
1123 : }
1124 :
1125 : //m_subdiag_error_est[0] = ;
1126 : // process columns (left to right)
1127 0 : for (size_t k=1; k<nstages; ++k)
1128 : {
1129 :
1130 : // process rows (bottom-up -> recycling memory)
1131 0 : for (size_t i=nstages-1; i>=k; --i)
1132 : {
1133 : UG_ASSERT(m_solution[i].valid(), "Invalid SmarPtr!");
1134 : UG_ASSERT(m_solution[i-1].valid(), "Invalid SmarPtr!");
1135 :
1136 0 : SmartPtr<vector_type> solcoarse = m_solution[i-1];
1137 : SmartPtr<vector_type> solfine = m_solution[i];
1138 :
1139 : // (2^p -1)
1140 : // m_solution[i] += (1.0/scal)*(m_solution[i]- m_solution[i-1]);
1141 0 : const number scaling = ((1.0*m_num_steps[i])/(1.0*m_num_steps[i-k])-1.0);
1142 0 : UG_LOG("scaling="<<i << ","<< k <<
1143 : ": ns["<<i<<"]="<< m_num_steps[i] <<
1144 : "ns2["<<i-k<<"]=" << m_num_steps[i-k] <<"=" << scaling << std::endl);
1145 :
1146 0 : if (with_error && (i==k))
1147 : {
1148 : // compute subdiagonal error estimate
1149 0 : m_subdiag->update(solfine, (1.0/scaling), solfine, solcoarse);
1150 : number subdiag_error_est=m_subdiag->get_current_estimate();
1151 :
1152 : //m_subdiag_error_est[k] = sqrt(subdiag_error_est*scaling);
1153 0 : m_subdiag_error_est[k]=subdiag_error_est; //*scaling;
1154 :
1155 0 : UG_LOG(" ErrorEst["<< k<<"]=" << m_subdiag_error_est[k] << ";" << std::endl);
1156 : }
1157 : else
1158 : {
1159 : // standard case
1160 0 : VecScaleAdd(*solfine, (1.0+1.0/scaling), *solfine,
1161 0 : -(1.0/scaling), *solcoarse);
1162 :
1163 : }
1164 : } // rows i
1165 :
1166 : } // columns k
1167 :
1168 0 : }
1169 :
1170 : /// apply for all stages
1171 0 : void apply()
1172 : {
1173 : //UG_ASSERT(m_num_steps.size() == m_solution.size(), "Dimensions do not match");
1174 : const size_t nstages = m_num_steps.size();
1175 0 : apply(nstages);
1176 0 : }
1177 :
1178 :
1179 :
1180 : protected:
1181 0 : number substep(size_t i) {return m_stepsize/m_num_steps[i];}
1182 :
1183 : private:
1184 : /** time step */
1185 : number m_stepsize;
1186 : static const int m_order=1;
1187 :
1188 : /** evaluation on last stage */
1189 : SmartPtr<ISubDiagErrorEst<vector_type> > m_subdiag;
1190 :
1191 : /** n_i: number of intermediate steps (per stage)*/
1192 : std::vector<size_t> m_num_steps;
1193 :
1194 : /** vector of solutions (per stage)*/
1195 : std::vector<SmartPtr<vector_type> > m_solution;
1196 :
1197 : /** sub-diagonal error estimate (per stage)*/
1198 : std::vector<number> m_subdiag_error_est;
1199 :
1200 :
1201 : };
1202 :
1203 :
1204 :
1205 : /** Estimate*/
1206 : /*
1207 : class TimeStepEstimator{
1208 :
1209 : public:
1210 : TimeStepEstimator(double tol) :
1211 : m_rho(0.9), m_gamma(2.0), m_tol(tol)
1212 : {
1213 :
1214 : }
1215 :
1216 : bool check(double eps, double &factor)
1217 : {
1218 : if (eps>=tol) return false;
1219 :
1220 : double val = (m_rho*m_tol)/eps;
1221 : factor = pow(val, m_gamma);
1222 :
1223 : }
1224 :
1225 :
1226 : private:
1227 : double m_rho;
1228 : double m_gamma;
1229 : double m_tol;
1230 : };
1231 : */
1232 :
1233 : }
1234 : #endif
|
