KASKADE 7 development version
errorEstimator_efficient.hh
Go to the documentation of this file.
1/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2/* */
3/* This file is part of the library KASKADE 7 */
4/* https://www.zib.de/research/projects/kaskade7-finite-element-toolbox */
5/* */
6/* Copyright (C) 2002-2012 Zuse Institute Berlin */
7/* */
8/* KASKADE 7 is distributed under the terms of the ZIB Academic License. */
9/* see $KASKADE/academic.txt */
10/* */
11/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
12
13#ifndef ALGORITHM_3_ERROR_ESTIMATOR_EE
14#define ALGORITHM_3_ERROR_ESTIMATOR_EE
15
16#include <algorithm>
17#include <iostream>
18#include <utility>
19#include <vector>
20
21#include <boost/timer/timer.hpp>
22#include <boost/fusion/include/at_c.hpp>
23
24#include "algorithm/newton_bridge.hh"
25#include "algorithm/adaptationStrategy.hh"
26#include "algorithm/errorDistribution.hh"
27#include "linalg/jacobiPreconditioner.hh"
28#include "fem/forEach.hh"
29#include "fem/hierarchicErrorEstimator.hh"
30#include "fem/variables.hh"
31#include "utilities/enums.hh"
32#include "errorEstimationTraits.hh"
33#include "lagrangeLinearization.hh"
34
35// forward declarations
36namespace Dune
37{
38 struct InverseOperatorResult;
39 template <class,int> class FieldVector;
40}
41
42namespace Kaskade
43{
44 struct Dummy{};
45
46 template <template <class,class,class,bool> class Functional, class VariableSetDescription, class ExtensionVariableSetDescription, class ExtensionSpace,
47 class NormFunctional, template <class> class RefinementStrategy = Adaptivity::ErrorEquilibration, bool lump=false, int components=1, class ReferenceSolution = Dummy, class ReferenceOperator=Dummy>
48 class YetAnotherHBErrorEstimator : public ErrorEstimatorBase<VariableSetDescription,typename ErrorDistribution<NormFunctional,ExtensionVariableSetDescription>::AnsatzVars::VariableSet,ExtensionSpace,RefinementStrategy>
49 {
50 template <class AnsatzVars, class TestVars, class OriginVars> using EstimatorFunctional = Functional<AnsatzVars,TestVars,OriginVars,false>;
51 template <class AnsatzVars, class TestVars, class OriginVars> using LumpedFunctional = Functional<AnsatzVars,TestVars,OriginVars,lump>;
52 typedef ErrorEstimationTraits<EstimatorFunctional,VariableSetDescription,ExtensionVariableSetDescription> Traits;
53 typedef ErrorEstimationTraits<LumpedFunctional,VariableSetDescription,ExtensionVariableSetDescription> LumpedTraits;
54 typedef ErrorEstimatorBase<VariableSetDescription,typename ErrorDistribution<NormFunctional,ExtensionVariableSetDescription>::AnsatzVars::VariableSet,ExtensionSpace,RefinementStrategy> Base;
55 using Base::squaredError;
56 using Base::mde;
57 using Base::errorDistribution;
58 using Base::extensionSpace;
59 public:
60 // assemblers
61
62 typedef typename Traits::Scalar Scalar;
63 typedef typename VariableSetDescription::Grid Grid;
64 static constexpr int dim = VariableSetDescription::Grid::dimension;
65
66 typedef typename ExtensionVariableSetDescription::GridView::template Codim<0>::Iterator CellIterator;
67 typedef typename ExtensionVariableSetDescription::GridView::template Codim<0>::Entity Cell;
68 typedef ErrorDistribution<NormFunctional,ExtensionVariableSetDescription> EnergyError;
69 typedef typename EnergyError::AnsatzVars::VariableSet ErrorRepresentation;
70
71
72 YetAnotherHBErrorEstimator(NormFunctional& normFunctional_, VariableSetDescription& variableSetDescription_, ExtensionVariableSetDescription& extensionVariableSetDescription_,
73 ExtensionSpace& extensionSpace_, Scalar fraction=0.7, bool verbose_=false)
74 : Base(extensionSpace_,fraction,verbose), normFunctional(normFunctional_), variableSetDescription(variableSetDescription_), extensionVariableSetDescription(extensionVariableSetDescription_),
75 squaredFraction(fraction*fraction), verbose(verbose_)
76 {
77 ass_EE.reset(new typename LumpedTraits::Assembler_EE(extensionVariableSetDescription.spaces));
78 ass_HE.reset(new typename Traits::Assembler_HE(variableSetDescription.spaces));
79 ass_EH.reset(new typename Traits::Assembler_EH(extensionVariableSetDescription.spaces));
80 ass_HH.reset(new typename Traits::Assembler_HH(variableSetDescription.spaces));
81 }
82
83 virtual ~YetAnotherHBErrorEstimator(){}
84
85 void operator()(AbstractLinearization const& lin, AbstractFunctionSpaceElement const& x_, AbstractFunctionSpaceElement const& dx_, int step, AbstractFunctionSpaceElement const&)
86 {
87 boost::timer::cpu_timer overallTimer;
88 using boost::fusion::at_c;
89 if(verbose) std::cout << "ERROR ESTIMATOR: Start." << std::endl;
90 Bridge::Vector<typename Traits::VarSet> const& xl = dynamic_cast<const Bridge::Vector<typename Traits::VarSet>&>(x_);
91 computeError(xl.get());
92 Bridge::Vector<typename Traits::ExtensionVarSet> xe(extensionVariableSetDescription);
93 boost::timer::cpu_timer atimer;
94 assembleAt(xl,xe);
95 if(verbose) std::cout << "ERROR ESTIMATOR: assembly time: " << boost::timer::format(atimer.elapsed()) << std::endl;
96 Bridge::Vector<typename Traits::VarSet> const& dxl = dynamic_cast<const Bridge::Vector<typename Traits::VarSet>&>(dx_);
97
98 /***************************************************************************************/
99 // state equation -> DLY
100 typename Traits::A_HH a_HH(*ass_HH,false);
101 typename Traits::A_HE a_HE(*ass_HE,false);
102 typename LumpedTraits::A_EE a_EE(*ass_EE,false);
103 typename Traits::B_HH b_HH(*ass_HH,false);
104 typename Traits::B_HE b_HE(*ass_HE,false);
105 typename Traits::B_EH b_EH(*ass_EH,false);
106 typename LumpedTraits::B_EE b_EE(*ass_EE,false);
107
108 typename Traits::ExtensionVectorP rhsPE(ass_EE->template rhs<Traits::adjointId,Traits::adjointId+1>());
109 rhsPE *= -1.0;
110 typename Traits::VectorY dy(at_c<Traits::stateId>(dxl.get().data).coefficients()), solYH = dy; solYH *= 0;
111 typename Traits::VectorU du(at_c<Traits::controlId>(dxl.get().data).coefficients());
112 typename Traits::ExtensionVectorY solYE(Traits::ExtensionVectorY_Initializer::init(extensionVariableSetDescription));
113
114 b_HE.applyscaleadd(-1.0,du,rhsPE);
115 a_HE.applyscaleadd(-1.0,dy,rhsPE);
116
117 JacobiPreconditionerForTriplets<Scalar,typename LumpedTraits::ExtensionVectorY,typename LumpedTraits::ExtensionVectorP> jacobi(a_EE);
118 jacobi.apply(solYE,rhsPE);
119
120
121 // error in adjoint equation
122 typename Traits::Lyy_HH lyy_HH(*ass_HH,false);
123 typename Traits::Lyy_HE lyy_HE(*ass_HE,false);
124 typename Traits::Lyy_EH lyy_EH(*ass_EH,false);
125 typename LumpedTraits::Lyy_EE lyy_EE(*ass_EE,false);
126 TransposedOperator<typename Traits::A_HE> at_EH(a_HE);//(*ass_EH,false);
127
128 typename Traits::VectorY rhsYH(ass_HH->template rhs<Traits::stateId,Traits::stateId+1>());
129 rhsYH *= -1.0;
130 typename Traits::ExtensionVectorY rhsYE(ass_EE->template rhs<Traits::stateId, Traits::stateId+1>());
131 rhsYE *= -1.0;
132
133 lyy_HH.applyscaleadd(-1.0,dy,rhsYH);
134 lyy_HE.applyscaleadd(-1.0,dy,rhsYE);
135 lyy_EH.applyscaleadd(-1.0,solYE,rhsYH);
136 lyy_EE.applyscaleadd(-1.0,solYE,rhsYE);
137
138 typename Traits::ExtensionVectorP solPE(Traits::ExtensionVectorP_Initializer::init(extensionVariableSetDescription));
139 typename Traits::VectorP solPH(Traits::VectorP_Initializer::init(variableSetDescription));
140
141 jacobi.apply(solPE,rhsYE);
142 at_EH.applyscaleadd(-1.0,solPE,rhsYH);
143
144 boost::timer::cpu_timer mgTimer;
145 MultiGridSolver<Grid,components>( a_HH, boost::fusion::at_c<0>(ass_HH->spaces())->gridManager().grid(),
146 typename MultiGridSolver<Grid,components>::Parameter(mgSteps,mgSmoothingSteps,relativeAccuracy), true ).apply(at_c<0>(solPH.data),
147 at_c<0>(rhsYH.data));
148 std::cout << "mg: " << boost::timer::format(mgTimer.elapsed()) << std::endl;
149
150 // variational equality
151 TransposedOperator<typename Traits::B_HH> bt_HH(b_HH);
152 TransposedOperator<typename Traits::B_EH> bt_HE(b_EH);
153 TransposedOperator<typename Traits::B_HE> bt_EH(b_HE);
154 TransposedOperator<typename LumpedTraits::B_EE> bt_EE(b_EE);
155 typename Traits::Luu_HH luu_HH(*ass_HH,false);
156 typename Traits::Luu_EH luu_EH(*ass_EH,false);
157 typename Traits::Luu_HE luu_HE(*ass_HE,false);
158 typename LumpedTraits::Luu_EE luu_EE(*ass_EE,false);
159
160 typename Traits::VectorU rhsUH(ass_HH->template rhs<Traits::controlId,Traits::controlId+1>());
161 rhsUH *= -1;
162 typename Traits::ExtensionVectorU rhsUE(ass_EE->template rhs<Traits::controlId,Traits::controlId+1>());
163 rhsUE *= -1;
164
165 luu_HH.applyscaleadd(-1.0,du,rhsUH);
166 luu_HE.applyscaleadd(-1.0,du,rhsUE);
167
168 bt_HH.applyscaleadd(-1.0,solPH,rhsUH);
169 bt_HE.applyscaleadd(-1.0,solPH,rhsUE);
170 bt_EH.applyscaleadd(-1.0,solPE,rhsUH);
171 bt_EE.applyscaleadd(-1.0,solPE,rhsUE);
172
173 typename Traits::VectorU solUH(Traits::VectorU_Initializer::init(variableSetDescription));
174 typename Traits::ExtensionVectorU solUE(Traits::ExtensionVectorU_Initializer::init(extensionVariableSetDescription));
175
176 JacobiPreconditionerForTriplets<Scalar,typename LumpedTraits::ExtensionVectorU,typename LumpedTraits::ExtensionVectorU>(luu_EE).apply(solUE,rhsUE);
177 luu_EH.applyscaleadd(-1.0,solUE,rhsUH);
178 boost::timer::cpu_timer chebTimer;
179 ChebyshevPreconditioner<typename Traits::Luu_HH> cheb(luu_HH,chebySteps);
180 cheb.initForMassMatrix_TetrahedralQ1Elements();
181 cheb.apply(solUH,rhsUH);
182 std::cout << "cheb: " << boost::timer::format(chebTimer.elapsed()) << std::endl;
183
184 typename ExtensionVariableSetDescription::VariableSet errorEstimate_E(extensionVariableSetDescription);
185 typename VariableSetDescription::VariableSet errorEstimate_H(variableSetDescription);
186 solYH *= 0.0; // DLY in state equation -> no error transport into space of linear finite elements
187 at_c<Traits::stateId>(errorEstimate_H.data) = at_c<0>(solYH.data);
188 at_c<Traits::controlId>(errorEstimate_H.data) = at_c<0>(solUH.data);
189 at_c<Traits::stateId>(errorEstimate_E.data) = at_c<0>(solYE.data);
190 at_c<Traits::controlId>(errorEstimate_E.data) = at_c<0>(solUE.data);
191// std::string savefilename = createFileName("estimate_E",".vtu",false);
192// std::cout << "Error estimator: writing estimates" << std::endl;
193// writeVTKFile(extensionVariableSetDescription.gridView,errorEstimate_E,savefilename,IoOptions(),2);
194// std::cout << "Error estimator: writing lower order estimates" << std::endl;
195// writeVTKFile(variableSetDescription.gridView,errorEstimate_H,savefilename,IoOptions(),1);
196 /***************************************************************************************/
197 // Transfer error indicators to cells.
198 auto const& is = extensionSpace.gridManager().grid().leafIndexSet();
199 errorDistribution.clear();
200 errorDistribution.resize(is.size(0),std::make_pair(0.0,0));
201
202 EnergyError energyError(normFunctional,xl.get(),errorEstimate_H,errorEstimate_E);
203 energyError.considerStateVariable(true);
204 energyError.considerControlVariable(true);
205 energyError.considerAdjointVariable(false);
206 typedef VariationalFunctionalAssembler<LinearizationAt<EnergyError> > EnergyErrorAssembler;
207 EnergyErrorAssembler eeAssembler(extensionSpace.gridManager(), energyError.getSpaces());
208
209 eeAssembler.assemble(linearization(energyError,xl.get()), EnergyErrorAssembler::RHS);
210 typename EnergyError::ErrorVector distError( eeAssembler.rhs() );
211 mde.reset(new ErrorRepresentation(energyError.getVariableSetDescription()) );
212 *mde = distError;
213
214// if(verbose)
215// {
216// std::string name = "errorDistribution_";
217// name += std::to_string(step);
218// writeVTKFile(mde->descriptions.gridView,*mde,name);
219// }
220
221 CellIterator cend = extensionVariableSetDescription.gridView.template end<0>();
222 for (CellIterator ci=extensionVariableSetDescription.gridView.template begin<0>(); ci!=cend; ++ci)
223 errorDistribution[is.index(*ci)] = std::make_pair( fabs(boost::fusion::at_c<0>(mde->data).value(*ci,Dune::FieldVector<Scalar,dim>(0.3))) , is.index(*ci));
224
225 squaredError = std::accumulate(errorDistribution.begin(), errorDistribution.end(), 0.0, ErrorEstimator_Detail::add);
226
227 if(verbose) std::cout << "overall error estimation time: " << boost::timer::format(overallTimer.elapsed()) << std::endl;
228 }
229
230 template <typename... Args>
231 void initFunctionals(const Args&... args)
232 {
233 F_HH.reset(new typename Traits::Functional_HH(args...));
234 F_HE.reset(new typename Traits::Functional_HE(args...));
235 }
236
237 template <typename... Args>
238 void initExtensionFunctionals(const Args&... args)
239 {
240 F_EH.reset(new typename Traits::Functional_EH(args...));
241 F_EE.reset(new typename LumpedTraits::Functional_EE(args...));
242 }
243
244 void setReference(ReferenceOperator const& Aref, ReferenceSolution const& ref, typename Traits::Vector& refcv)
245 {
246 referenceSolution = &ref;
247 referenceOperator = &Aref;
248 referenceCoeffVec = &refcv;
249 }
250
251 void computeError(typename Traits::VarSet const& x)
252 {
253 if(referenceSolution == nullptr) return;
254 ReferenceSolution tmp(*referenceSolution);
255 interpolateGloballyFromFunctor<PlainAverage>(boost::fusion::at_c<0>(tmp.data),[&x](Cell const& cell, Dune::FieldVector<Scalar,dim> const& xLocal)
256 {
257 return boost::fusion::at_c<0>(x.data).value(cell.geometry().global(xLocal));
258 });
259
260 interpolateGloballyFromFunctor<PlainAverage>(boost::fusion::at_c<1>(tmp.data),[&x](Cell const& cell, Dune::FieldVector<Scalar,dim> const& xLocal)
261 {
262 return boost::fusion::at_c<1>(x.data).value(cell.geometry().global(xLocal));
263 });
264// interpolateGloballyWeak<PlainAverage>(boost::fusion::at_c<2>(tmp.data),boost::fusion::at_c<2>(x.data));
265 tmp -= *referenceSolution;
266
267 boost::fusion::at_c<0>(referenceCoeffVec->data) = boost::fusion::at_c<0>(tmp.data).coefficients();
268 boost::fusion::at_c<1>(referenceCoeffVec->data) = boost::fusion::at_c<1>(tmp.data).coefficients();
269 // boost::fusion::at_c<2>(referenceCoeffVec->data) = boost::fusion::at_c<2>(tmp.data).coefficients();
270
271 auto tmp2 = *referenceCoeffVec;
272 referenceOperator->apply(*referenceCoeffVec,tmp2);
273 std::cout << "REAL ERROR: " << (*referenceCoeffVec)*tmp2 << std::endl;
274 }
275
276 private:
277 void assembleAt(const Bridge::Vector<typename Traits::VarSet>& xl, const Bridge::Vector<typename Traits::ExtensionVarSet>& xe)
278 {
279 ass_HH->assemble(linearization(*F_HH,xl.get()));
280 ass_HE->assemble(linearization(*F_HE,xl.get()));
281 ass_EH->assemble(linearization(*F_EH,xe.get()));
282 ass_EE->assemble(linearization(*F_EE,xe.get()));
283 }
284
285 std::unique_ptr<typename Traits::Functional_HH> F_HH;
286 std::unique_ptr<typename Traits::Functional_HE> F_HE;
287 std::unique_ptr<typename Traits::Functional_EH> F_EH;
288 std::unique_ptr<typename LumpedTraits::Functional_EE> F_EE;
289 std::unique_ptr<typename Traits::Assembler_HH> ass_HH;
290 //typename Traits::Assembler_HH* ass_HH;
291 std::unique_ptr<typename Traits::Assembler_HE> ass_HE;
292 std::unique_ptr<typename Traits::Assembler_EH> ass_EH;
293 std::unique_ptr<typename LumpedTraits::Assembler_EE> ass_EE;
294 NormFunctional& normFunctional;
295 VariableSetDescription& variableSetDescription;
296 ExtensionVariableSetDescription& extensionVariableSetDescription;
297 Scalar squaredFraction;
298 bool verbose;
299 size_t chebySteps = 30, mgSteps = 25, mgSmoothingSteps = 20;
300 Scalar relativeAccuracy = 1e-3;
301 ReferenceSolution const* referenceSolution = nullptr;
302 ReferenceOperator const* referenceOperator = nullptr;
303 typename Traits::Vector* referenceCoeffVec = nullptr;
304 //std::function<Scalar(typename Traits::VarSet)> compareError;
305 };
306
307 template <template <class,class,class,bool> class Functional, class VariableSetDescription, class ExtensionVariableSetDescription, class ExtensionSpace,
308 class NormFunctional, template <class> class RefinementStrategy = Adaptivity::ErrorEquilibration, bool lump=false, int components=1>
309 class YetAnotherHBErrorEstimator_Elasticity : public ErrorEstimatorBase<VariableSetDescription,typename ErrorDistribution<NormFunctional,ExtensionVariableSetDescription>::AnsatzVars::VariableSet,ExtensionSpace,RefinementStrategy>
310 {
311 template <class AnsatzVars, class TestVars, class OriginVars> using EstimatorFunctional = Functional<AnsatzVars,TestVars,OriginVars,false>;
312 template <class AnsatzVars, class TestVars, class OriginVars> using LumpedFunctional = Functional<AnsatzVars,TestVars,OriginVars,lump>;
313 typedef ErrorEstimationTraits<EstimatorFunctional,VariableSetDescription,ExtensionVariableSetDescription> Traits;
314 typedef ErrorEstimationTraits<LumpedFunctional,VariableSetDescription,ExtensionVariableSetDescription> LumpedTraits;
315 typedef ErrorEstimatorBase<VariableSetDescription,typename ErrorDistribution<NormFunctional,ExtensionVariableSetDescription>::AnsatzVars::VariableSet,ExtensionSpace,RefinementStrategy> Base;
316 using Base::squaredError;
317 using Base::mde;
318 using Base::errorDistribution;
319 using Base::extensionSpace;
320 public:
321 // assemblers
322
323 typedef typename Traits::Scalar Scalar;
324 typedef typename VariableSetDescription::Grid Grid;
325 static constexpr int dim = VariableSetDescription::Grid::dimension;
326
327 typedef typename ExtensionVariableSetDescription::GridView::template Codim<0>::Iterator CellIterator;
328 typedef typename ExtensionVariableSetDescription::GridView::template Codim<0>::Entity Cell;
329 typedef ErrorDistribution<NormFunctional,ExtensionVariableSetDescription> EnergyError;
330 typedef typename EnergyError::AnsatzVars::VariableSet ErrorRepresentation;
331
332
333 YetAnotherHBErrorEstimator_Elasticity(NormFunctional& normFunctional_, VariableSetDescription& variableSetDescription_, ExtensionVariableSetDescription& extensionVariableSetDescription_,
334 ExtensionSpace& extensionSpace_, Scalar fraction=0.7, bool verbose_=false)
335 : Base(extensionSpace_,fraction,verbose), normFunctional(normFunctional_), variableSetDescription(variableSetDescription_), extensionVariableSetDescription(extensionVariableSetDescription_),
336 squaredFraction(fraction*fraction), verbose(verbose_)
337 {
338 ass_EE.reset(new typename LumpedTraits::Assembler_EE(extensionVariableSetDescription.spaces));
339 ass_HE.reset(new typename Traits::Assembler_HE(variableSetDescription.spaces));
340 ass_EH.reset(new typename Traits::Assembler_EH(extensionVariableSetDescription.spaces));
341 ass_HH.reset(new typename Traits::Assembler_HH(variableSetDescription.spaces));
342 }
343
344 virtual ~YetAnotherHBErrorEstimator_Elasticity(){}
345
346 void operator()(AbstractLinearization const& lin, AbstractFunctionSpaceElement const& x_, AbstractFunctionSpaceElement const& dx_, int step, AbstractFunctionSpaceElement const&)
347 {
348 boost::timer::cpu_timer overallTimer;
349 using boost::fusion::at_c;
350 if(verbose) std::cout << "ERROR ESTIMATOR: Start." << std::endl;
351 Bridge::Vector<typename Traits::VarSet> const& xl = dynamic_cast<const Bridge::Vector<typename Traits::VarSet>&>(x_);
352 Bridge::Vector<typename Traits::ExtensionVarSet> xe(extensionVariableSetDescription);
353 boost::timer::cpu_timer atimer;
354 assembleAt(xl,xe);
355 if(verbose) std::cout << "ERROR ESTIMATOR: assembly time: " << boost::timer::format(atimer.elapsed()) << std::endl;
356 Bridge::Vector<typename Traits::VarSet> const& dxl = dynamic_cast<const Bridge::Vector<typename Traits::VarSet>&>(dx_);
357
358 /***************************************************************************************/
359 // state equation -> DLY
360 typename Traits::A_HH a_HH(*ass_HH,false);
361 typename Traits::A_HE a_HE(*ass_HE,false);
362 typename LumpedTraits::A_EE a_EE(*ass_EE,false);
363 typename Traits::B_HH b_HH(*ass_HH,false);
364 typename Traits::B_HE b_HE(*ass_HE,false);
365 typename Traits::B_EH b_EH(*ass_EH,false);
366 typename LumpedTraits::B_EE b_EE(*ass_EE,false);
367
368 typename Traits::ExtensionVectorP rhsPE(ass_EE->template rhs<Traits::adjointId,Traits::adjointId+1>());
369 rhsPE *= -1.0;
370 typename Traits::VectorY dy(at_c<Traits::stateId>(dxl.get().data).coefficients()), solYH = dy; solYH *= 0;
371 typename Traits::VectorU du(at_c<Traits::controlId>(dxl.get().data).coefficients());
372 typename Traits::ExtensionVectorY solYE(Traits::ExtensionVectorY_Initializer::init(extensionVariableSetDescription));
373
374 b_HE.applyscaleadd(-1.0,du,rhsPE);
375 a_HE.applyscaleadd(-1.0,dy,rhsPE);
376
377 JacobiPreconditionerForTriplets<Scalar,typename LumpedTraits::ExtensionVectorY,typename LumpedTraits::ExtensionVectorP> jacobi(a_EE);
378 jacobi.apply(solYE,rhsPE);
379
380
381 // error in adjoint equation
382 typename Traits::Lyy_HH lyy_HH(*ass_HH,false);
383 typename Traits::Lyy_HE lyy_HE(*ass_HE,false);
384 typename Traits::Lyy_EH lyy_EH(*ass_EH,false);
385 typename LumpedTraits::Lyy_EE lyy_EE(*ass_EE,false);
386 typename Traits::Luy_HH luy_HH(*ass_HH,false);
387 typename Traits::Luy_HE luy_HE(*ass_HE,false);
388 typename Traits::Luy_EH luy_EH(*ass_EH,false);
389 typename LumpedTraits::Luy_EE luy_EE(*ass_EE,false);
390 TransposedOperator<typename Traits::Luy_HH> lyu_HH(luy_HH);
391 TransposedOperator<typename Traits::Luy_EH> lyu_HE(luy_EH);
392// typename Traits::Lyu_HH lyu_HH(*ass_HH,false);
393// typename Traits::Lyu_HE lyu_HE(*ass_HE,false);
394 TransposedOperator<typename Traits::A_HE> at_EH(a_HE);//(*ass_EH,false);
395
396 typename Traits::VectorY rhsYH(ass_HH->template rhs<Traits::stateId,Traits::stateId+1>());
397 rhsYH *= -1.0;
398 typename Traits::ExtensionVectorY rhsYE(ass_EE->template rhs<Traits::stateId, Traits::stateId+1>());
399 rhsYE *= -1.0;
400
401 lyy_HH.applyscaleadd(-1.0,dy,rhsYH);
402 lyy_HE.applyscaleadd(-1.0,dy,rhsYE);
403 lyy_EH.applyscaleadd(-1.0,solYE,rhsYH);
404 lyy_EE.applyscaleadd(-1.0,solYE,rhsYE);
405 lyu_HH.applyscaleadd(-1.0,du,rhsYH);
406 lyu_HE.applyscaleadd(-1.0,du,rhsYE);
407
408 typename Traits::ExtensionVectorP solPE(Traits::ExtensionVectorP_Initializer::init(extensionVariableSetDescription));
409 typename Traits::VectorP solPH(Traits::VectorP_Initializer::init(variableSetDescription));
410
411 jacobi.apply(solPE,rhsYE);
412 at_EH.applyscaleadd(-1.0,solPE,rhsYH);
413
414 boost::timer::cpu_timer mgTimer;
415 MultiGridSolver<Grid,components>( a_HH, boost::fusion::at_c<0>(ass_HH->spaces())->gridManager().grid(),
416 typename MultiGridSolver<Grid,components>::Parameter(mgSteps,mgSmoothingSteps,relativeAccuracy), true ).apply(at_c<0>(solPH.data),
417 at_c<0>(rhsYH.data));
418 std::cout << "mg: " << boost::timer::format(mgTimer.elapsed()) << std::endl;
419
420
421 // variational equality
422 TransposedOperator<typename Traits::B_HH> bt_HH(b_HH);
423 TransposedOperator<typename Traits::B_EH> bt_HE(b_EH);
424 TransposedOperator<typename Traits::B_HE> bt_EH(b_HE);
425 TransposedOperator<typename LumpedTraits::B_EE> bt_EE(b_EE);
426 typename Traits::Luu_HH luu_HH(*ass_HH,false);
427 typename Traits::Luu_EH luu_EH(*ass_EH,false);
428 typename Traits::Luu_HE luu_HE(*ass_HE,false);
429 typename LumpedTraits::Luu_EE luu_EE(*ass_EE,false);
430
431 typename Traits::VectorU rhsUH(ass_HH->template rhs<Traits::controlId,Traits::controlId+1>());
432 rhsUH *= -1;
433 typename Traits::ExtensionVectorU rhsUE(ass_EE->template rhs<Traits::controlId,Traits::controlId+1>());
434 rhsUE *= -1;
435
436 luu_HH.applyscaleadd(-1.0,du,rhsUH);
437 luu_HE.applyscaleadd(-1.0,du,rhsUE);
438
439 bt_HH.applyscaleadd(-1.0,solPH,rhsUH);
440 bt_HE.applyscaleadd(-1.0,solPH,rhsUE);
441 bt_EH.applyscaleadd(-1.0,solPE,rhsUH);
442 bt_EE.applyscaleadd(-1.0,solPE,rhsUE);
443
444 luy_HH.applyscaleadd(-1.0,dy,rhsUH);
445 luy_HE.applyscaleadd(-1.0,dy,rhsUE);
446 luy_EH.applyscaleadd(-1.0,solYE,rhsUH);
447 luy_EE.applyscaleadd(-1.0,solYE,rhsUE);
448
449 typename Traits::VectorU solUH(Traits::VectorU_Initializer::init(variableSetDescription));
450 typename Traits::ExtensionVectorU solUE(Traits::ExtensionVectorU_Initializer::init(extensionVariableSetDescription));
451
452 JacobiPreconditionerForTriplets<Scalar,typename LumpedTraits::ExtensionVectorU,typename LumpedTraits::ExtensionVectorU>(luu_EE).apply(solUE,rhsUE);
453 luu_EH.applyscaleadd(-1.0,solUE,rhsUH);
454 boost::timer::cpu_timer chebTimer;
455 ChebyshevPreconditioner<typename Traits::Luu_HH> cheb(luu_HH,chebySteps);
456 cheb.initForMassMatrix_TetrahedralQ1Elements();
457 cheb.apply(solUH,rhsUH);
458 std::cout << "cheb: " << boost::timer::format(chebTimer.elapsed()) << std::endl;
459
460 typename ExtensionVariableSetDescription::VariableSet errorEstimate_E(extensionVariableSetDescription);
461 typename VariableSetDescription::VariableSet errorEstimate_H(variableSetDescription);
462 solYH *= 0.0; // DLY in state equation -> no error transport into space of linear finite elements
463 at_c<Traits::stateId>(errorEstimate_H.data) = at_c<0>(solYH.data);
464 at_c<Traits::controlId>(errorEstimate_H.data) = at_c<0>(solUH.data);
465 at_c<Traits::stateId>(errorEstimate_E.data) = at_c<0>(solYE.data);
466 at_c<Traits::controlId>(errorEstimate_E.data) = at_c<0>(solUE.data);
467// std::string savefilename = createFileName("estimate_E",".vtu",false);
468// std::cout << "Error estimator: writing estimates" << std::endl;
469// writeVTKFile(extensionVariableSetDescription.gridView,errorEstimate_E,savefilename,IoOptions(),2);
470// std::cout << "Error estimator: writing lower order estimates" << std::endl;
471// writeVTKFile(variableSetDescription.gridView,errorEstimate_H,savefilename,IoOptions(),1);
472 /***************************************************************************************/
473 // Transfer error indicators to cells.
474 auto const& is = extensionSpace.gridManager().grid().leafIndexSet();
475 errorDistribution.clear();
476 errorDistribution.resize(is.size(0),std::make_pair(0.0,0));
477
478 EnergyError energyError(normFunctional,xl.get(),errorEstimate_H,errorEstimate_E);
479 energyError.considerStateVariable(true);
480 energyError.considerControlVariable(true);
481 energyError.considerAdjointVariable(false);
482 typedef VariationalFunctionalAssembler<LinearizationAt<EnergyError> > EnergyErrorAssembler;
483 EnergyErrorAssembler eeAssembler(extensionSpace.gridManager(), energyError.getSpaces());
484
485 eeAssembler.assemble(linearization(energyError,xl.get()), EnergyErrorAssembler::RHS);
486 typename EnergyError::ErrorVector distError( eeAssembler.rhs() );
487 mde.reset(new ErrorRepresentation(energyError.getVariableSetDescription()) );
488 *mde = distError;
489
490// if(verbose)
491// {
492// std::string name = "errorDistribution_";
493// name += std::to_string(step);
494// writeVTKFile(mde->descriptions.gridView,*mde,name);
495// }
496
497 CellIterator cend = extensionVariableSetDescription.gridView.template end<0>();
498 for (CellIterator ci=extensionVariableSetDescription.gridView.template begin<0>(); ci!=cend; ++ci)
499 errorDistribution[is.index(*ci)] = std::make_pair( fabs(boost::fusion::at_c<0>(mde->data).value(*ci,Dune::FieldVector<Scalar,dim>(0.3))) , is.index(*ci));
500
501 squaredError = std::accumulate(errorDistribution.begin(), errorDistribution.end(), 0.0, ErrorEstimator_Detail::add);
502
503 if(verbose) std::cout << "overall error estimation time: " << boost::timer::format(overallTimer.elapsed()) << std::endl;
504 }
505
506 template <typename... Args>
507 void initFunctionals(const Args&... args)
508 {
509 F_HH.reset(new typename Traits::Functional_HH(args...));
510 F_HE.reset(new typename Traits::Functional_HE(args...));
511 }
512
513 template <typename... Args>
514 void initExtensionFunctionals(const Args&... args)
515 {
516 F_EH.reset(new typename Traits::Functional_EH(args...));
517 F_EE.reset(new typename LumpedTraits::Functional_EE(args...));
518 }
519
520 private:
521 void assembleAt(const Bridge::Vector<typename Traits::VarSet>& xl, const Bridge::Vector<typename Traits::ExtensionVarSet>& xe)
522 {
523 ass_HH->assemble(linearization(*F_HH,xl.get()));
524 ass_HE->assemble(linearization(*F_HE,xl.get()));
525 ass_EH->assemble(linearization(*F_EH,xe.get()));
526 ass_EE->assemble(linearization(*F_EE,xe.get()));
527 }
528
529 std::unique_ptr<typename Traits::Functional_HH> F_HH;
530 std::unique_ptr<typename Traits::Functional_HE> F_HE;
531 std::unique_ptr<typename Traits::Functional_EH> F_EH;
532 std::unique_ptr<typename LumpedTraits::Functional_EE> F_EE;
533 std::unique_ptr<typename Traits::Assembler_HH> ass_HH;
534 //typename Traits::Assembler_HH* ass_HH;
535 std::unique_ptr<typename Traits::Assembler_HE> ass_HE;
536 std::unique_ptr<typename Traits::Assembler_EH> ass_EH;
537 std::unique_ptr<typename LumpedTraits::Assembler_EE> ass_EE;
538 NormFunctional& normFunctional;
539 VariableSetDescription& variableSetDescription;
540 ExtensionVariableSetDescription& extensionVariableSetDescription;
541 Scalar squaredFraction;
542 bool verbose;
543 size_t chebySteps = 30, mgSteps = 25, mgSmoothingSteps = 20;
544 Scalar relativeAccuracy = 1e-3;
545 };
546}
547
548#endif
Kaskade::AbstractFunctionSpaceElement
Abstract Vector for function space algorithms.
Definition: abstract_interface.hh:26
Kaskade::AssembledGalerkinOperator
An adaptor presenting a Dune::LinearOperator <domain_type,range_type> interface to a contiguous sub-b...
Definition: istlinterface.hh:758
Kaskade::AssembledGalerkinOperator::applyscaleadd
virtual void applyscaleadd(Scalar alpha, Domain const &x, Range &b) const
Compute .
Definition: istlinterface.hh:938
Kaskade::Bridge::Vector
Mathematical Vector that supports copy-on-write, implements AbstractFunctionSpaceElement.
Definition: newton_bridge.hh:196
Kaskade::Bridge::Vector::get
Implementation const & get() const
Access to the data.
Definition: newton_bridge.hh:216
Kaskade::ChebyshevPreconditioner
Preconditioner based on the Chebyshev semi-iteration. In order to provide a linear preconditioner the...
Definition: chebyshev.hh:160
Kaskade::ChebyshevSemiIteration::initForMassMatrix_TetrahedralQ1Elements
void initForMassMatrix_TetrahedralQ1Elements()
Init spectral bounds for the mass matrix arising from tetrahedral discretization of the domain and li...
Definition: chebyshev.hh:133
Kaskade::ChebyshevSemiIteration::apply
virtual void apply(Domain &x, Range const &y)
Definition: chebyshev.hh:68
Kaskade::ErrorDistribution::ErrorVector
AnsatzVars::template CoefficientVectorRepresentation< 0, 1 >::type ErrorVector
Definition: errorDistribution.hh:115
Kaskade::ErrorDistribution::getSpaces
AnsatzSpaces const & getSpaces() const
Definition: errorDistribution.hh:363
Kaskade::ErrorDistribution::getVariableSetDescription
AnsatzVars const & getVariableSetDescription() const
Definition: errorDistribution.hh:365
Kaskade::ErrorDistribution::considerStateVariable
void considerStateVariable(bool consider)
Definition: errorDistribution.hh:368
Kaskade::ErrorDistribution::considerAdjointVariable
void considerAdjointVariable(bool consider)
Definition: errorDistribution.hh:370
Kaskade::ErrorDistribution::considerControlVariable
void considerControlVariable(bool consider)
Definition: errorDistribution.hh:369
Kaskade::ErrorEstimatorBase::mde
std::unique_ptr< ErrorRepresentation > mde
Definition: errorEstimationTraits.hh:151
Kaskade::ErrorEstimatorBase::errorDistribution
std::vector< std::pair< double, int > > errorDistribution
Definition: errorEstimationTraits.hh:150
Kaskade::JacobiPreconditionerForTriplets::apply
void apply(Domain &solution, Range const &rhs)
Definition: jacobiSolver.hh:218
Kaskade::MultigridSolver::apply
void apply(CoeffVector &solution, CoeffVector &rightHand)
Definition: multiGridSolver.hh:406
Kaskade::TransposedOperator::applyscaleadd
virtual void applyscaleadd(Scalar alpha, Domain const &x, Range &b) const
Definition: istlinterface.hh:1447
Kaskade::VariableSetDescription
A class that stores information about a set of variables.
Definition: variables.hh:537
Kaskade::VariableSetDescription< AnsatzSpaces, VariableDescriptions >::VariableSet
Kaskade::VariableSet< VariableSetDescription< Spaces, VariableDescriptions > > VariableSet
Type that contains a set of variable values with some functionality, such as simple vector arithmetic...
Definition: variables.hh:566
Kaskade::VariableSetDescription::spaces
Spaces spaces
A heterogeneous sequence of pointers to (const) spaces involved.
Definition: variables.hh:807
Kaskade::VariableSetDescription::Grid
SpaceType< Spaces, 0 >::type::Grid Grid
Grid type.
Definition: variables.hh:571
Kaskade::VariableSet
A class for storing a heterogeneous collection of FunctionSpaceElement s.
Definition: variables.hh:341
Kaskade::VariationalFunctionalAssembler
A class for assembling Galerkin representation matrices and right hand sides for variational function...
Definition: assemblerCore.hh:1432
Kaskade::YetAnotherHBErrorEstimator_Elasticity
Definition: errorEstimator_efficient.hh:310
Kaskade::YetAnotherHBErrorEstimator_Elasticity::ErrorRepresentation
EnergyError::AnsatzVars::VariableSet ErrorRepresentation
Definition: errorEstimator_efficient.hh:330
Kaskade::YetAnotherHBErrorEstimator_Elasticity::Scalar
Traits::Scalar Scalar
Definition: errorEstimator_efficient.hh:323
Kaskade::YetAnotherHBErrorEstimator_Elasticity::Grid
VariableSetDescription::Grid Grid
Definition: errorEstimator_efficient.hh:324
Kaskade::YetAnotherHBErrorEstimator_Elasticity::initFunctionals
void initFunctionals(const Args &... args)
Definition: errorEstimator_efficient.hh:507
Kaskade::YetAnotherHBErrorEstimator_Elasticity::CellIterator
ExtensionVariableSetDescription::GridView::template Codim< 0 >::Iterator CellIterator
Definition: errorEstimator_efficient.hh:327
Kaskade::YetAnotherHBErrorEstimator_Elasticity::dim
static constexpr int dim
Definition: errorEstimator_efficient.hh:325
Kaskade::YetAnotherHBErrorEstimator_Elasticity::Cell
ExtensionVariableSetDescription::GridView::template Codim< 0 >::Entity Cell
Definition: errorEstimator_efficient.hh:328
Kaskade::YetAnotherHBErrorEstimator_Elasticity::~YetAnotherHBErrorEstimator_Elasticity
virtual ~YetAnotherHBErrorEstimator_Elasticity()
Definition: errorEstimator_efficient.hh:344
Kaskade::YetAnotherHBErrorEstimator_Elasticity::initExtensionFunctionals
void initExtensionFunctionals(const Args &... args)
Definition: errorEstimator_efficient.hh:514
Kaskade::YetAnotherHBErrorEstimator_Elasticity::YetAnotherHBErrorEstimator_Elasticity
YetAnotherHBErrorEstimator_Elasticity(NormFunctional &normFunctional_, VariableSetDescription &variableSetDescription_, ExtensionVariableSetDescription &extensionVariableSetDescription_, ExtensionSpace &extensionSpace_, Scalar fraction=0.7, bool verbose_=false)
Definition: errorEstimator_efficient.hh:333
Kaskade::YetAnotherHBErrorEstimator_Elasticity::operator()
void operator()(AbstractLinearization const &lin, AbstractFunctionSpaceElement const &x_, AbstractFunctionSpaceElement const &dx_, int step, AbstractFunctionSpaceElement const &)
Definition: errorEstimator_efficient.hh:346
Kaskade::YetAnotherHBErrorEstimator_Elasticity::EnergyError
ErrorDistribution< NormFunctional, ExtensionVariableSetDescription > EnergyError
Definition: errorEstimator_efficient.hh:329
Kaskade::YetAnotherHBErrorEstimator::ErrorRepresentation
EnergyError::AnsatzVars::VariableSet ErrorRepresentation
Definition: errorEstimator_efficient.hh:69
Kaskade::YetAnotherHBErrorEstimator::setReference
void setReference(ReferenceOperator const &Aref, ReferenceSolution const &ref, typename Traits::Vector &refcv)
Definition: errorEstimator_efficient.hh:244
Kaskade::YetAnotherHBErrorEstimator::YetAnotherHBErrorEstimator
YetAnotherHBErrorEstimator(NormFunctional &normFunctional_, VariableSetDescription &variableSetDescription_, ExtensionVariableSetDescription &extensionVariableSetDescription_, ExtensionSpace &extensionSpace_, Scalar fraction=0.7, bool verbose_=false)
Definition: errorEstimator_efficient.hh:72
Kaskade::YetAnotherHBErrorEstimator::EnergyError
ErrorDistribution< NormFunctional, ExtensionVariableSetDescription > EnergyError
Definition: errorEstimator_efficient.hh:68
Kaskade::YetAnotherHBErrorEstimator::initExtensionFunctionals
void initExtensionFunctionals(const Args &... args)
Definition: errorEstimator_efficient.hh:238
Kaskade::YetAnotherHBErrorEstimator::computeError
void computeError(typename Traits::VarSet const &x)
Definition: errorEstimator_efficient.hh:251
Kaskade::YetAnotherHBErrorEstimator::Cell
ExtensionVariableSetDescription::GridView::template Codim< 0 >::Entity Cell
Definition: errorEstimator_efficient.hh:67
Kaskade::YetAnotherHBErrorEstimator::CellIterator
ExtensionVariableSetDescription::GridView::template Codim< 0 >::Iterator CellIterator
Definition: errorEstimator_efficient.hh:66
Kaskade::YetAnotherHBErrorEstimator::operator()
void operator()(AbstractLinearization const &lin, AbstractFunctionSpaceElement const &x_, AbstractFunctionSpaceElement const &dx_, int step, AbstractFunctionSpaceElement const &)
Definition: errorEstimator_efficient.hh:85
hierarchicErrorEstimator.hh
Error estimation via hierachic FEM.
Kaskade::Assembler::RHS
@ RHS
Definition: assemblerCore.hh:1402
Kaskade::ErrorEstimator_Detail::add
double add(const double &x1, const std::pair< double, int > &x2)
Definition: errorEstimationTraits.hh:17
newton_bridge.hh
Bridge classes that connect low level FE algorithms to higher level algorithms.
Kaskade::ErrorEstimationTraits::Functional_EH
Functional< ExtensionVariableSetDescription, VariableSetDescription, ExtensionVariableSetDescription > Functional_EH
Definition: errorEstimationTraits.hh:31
Kaskade::ErrorEstimationTraits::ExtensionVectorY
ExtensionVariableSetDescription::template CoefficientVectorRepresentation< stateId, stateId+1 >::type ExtensionVectorY
Definition: errorEstimationTraits.hh:116
Kaskade::ErrorEstimationTraits::Functional_EE
Functional< ExtensionVariableSetDescription, ExtensionVariableSetDescription, ExtensionVariableSetDescription > Functional_EE
Definition: errorEstimationTraits.hh:32
Kaskade::ErrorEstimationTraits::ExtensionVectorP
ExtensionVariableSetDescription::template CoefficientVectorRepresentation< adjointId, adjointId+1 >::type ExtensionVectorP
Definition: errorEstimationTraits.hh:118
Kaskade::ErrorEstimationTraits::Functional_HH
Functional< VariableSetDescription, VariableSetDescription, VariableSetDescription > Functional_HH
Definition: errorEstimationTraits.hh:29
Kaskade::ErrorEstimationTraits::VectorU
VariableSetDescription::template CoefficientVectorRepresentation< controlId, controlId+1 >::type VectorU
Definition: errorEstimationTraits.hh:113
Kaskade::ErrorEstimationTraits::Functional_HE
Functional< VariableSetDescription, ExtensionVariableSetDescription, VariableSetDescription > Functional_HE
Definition: errorEstimationTraits.hh:30
Kaskade::ErrorEstimationTraits::ExtensionVectorU
ExtensionVariableSetDescription::template CoefficientVectorRepresentation< controlId, controlId+1 >::type ExtensionVectorU
Definition: errorEstimationTraits.hh:117
Kaskade::ErrorEstimationTraits::VectorY
VariableSetDescription::template CoefficientVectorRepresentation< stateId, stateId+1 >::type VectorY
Definition: errorEstimationTraits.hh:112
Kaskade::ErrorEstimationTraits::VectorP
VariableSetDescription::template CoefficientVectorRepresentation< adjointId, adjointId+1 >::type VectorP
Definition: errorEstimationTraits.hh:114
Kaskade::ErrorEstimationTraits::Vector
VariableSetDescription::template CoefficientVectorRepresentation ::type Vector
Definition: errorEstimationTraits.hh:111
variables.hh
Variables and their descriptions.
Generated on Fri Sep 6 2024 13:12:58 for KASKADE 7 development version by doxygen 1.9.4