KASKADE 7 development version
check_derivative.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-2011 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#ifndef CHECK_DERIVATIVE_HH
13#define CHECK_DERIVATIVE_HH
14
15#include <algorithm>
16#include <cassert>
17#include <cmath>
18#include <fstream>
19#include <string>
20#include <utility>
21#include <vector>
22
23#include <boost/mpl/size.hpp>
24#include <boost/lexical_cast.hpp>
25
26#include "algorithm/dune_bridge.hh"
27#include "utilities/linalg/scalarproducts.hh"
28
29namespace Kaskade{
31namespace XMLHelper{
32
33 static int const xml_increment=2;
34 static int xml_level;
35
36 std::string getXMLIndent()
37 {
38 std::string result;
39 for(int i=0; i<xml_level; ++i) result.append(" ");
40 return result;
41 }
42
43 struct XMLBlock{
44 XMLBlock(std::ofstream &f, const char* s, const char* d="") : file(f), str(s), data(d), offset("")
45 {
46 xml_level += xml_increment;
47 getOffset();
48 writeBegin();
49 }
50 XMLBlock(std::ofstream &f, std::string &s, std::string d="") : file(f), str(s), data(d), offset("")
51 {
52 xml_level += xml_increment;
53 getOffset();
54 writeBegin();
55 }
56 ~XMLBlock()
57 {
58 file << offset.c_str() << "</" << str.c_str() << ">" << std::endl;
59 xml_level -= xml_increment;
60 }
61
62 private:
63 void getOffset(){
64 for(int i=0; i<xml_level; ++i) offset.append(" ");
65 }
66
67 void writeBegin(){
68 file << offset.c_str() << "<" << str.c_str() << data.c_str() << ">" << std::endl;
69 }
70
71 std::ofstream& file;
72 std::string const str, data;
73 std::string offset;
74 };
75
76}
77
78namespace DerivativeCheck
79{
80 namespace Private
81 {
82 template <class SparseInt>
83 void writeConnectivity(std::ofstream &file, SparseInt row, SparseInt col, SparseInt numRows)
84 {
85 int const upperLeftVertexId = row+col*numRows;
86 file << upperLeftVertexId << " " << upperLeftVertexId+1 << " " << upperLeftVertexId+numRows << " " << upperLeftVertexId+numRows+1 << " ";
87 }
88 }
89
90 template <class Assembler, class Functional,int firstBlock=0, int lastBlock = Functional::AnsatzVars::noOfVariables>
91 typename Functional::Scalar d1(Assembler& assembler, Functional const& f, typename Functional::AnsatzVars::VariableSet& x, typename Functional::Scalar tolerance = 1e-6, typename Functional::Scalar increment = 1e-12, bool toFile=false, std::string const& filename = std::string("d1error"))
92 {
93 typedef typename Functional::Scalar Scalar;
94 typedef typename Functional::AnsatzVars::template CoefficientVectorRepresentation<firstBlock,lastBlock>::type CoefficientVector;
95
96 std::cout << "checking first derivative" << std::endl;
97 std::cout << "assembling...";
98
99 // storage for the evaluation of the derivative at x and the at x+increment*e_i
100 Scalar distorted;
101
102 // get reference solution
103 assembler.assemble(linearization(f,x), Assembler::RHS | Assembler::VALUE );
104
105 CoefficientVector analytic(assembler.template rhs<firstBlock,lastBlock>());
106 std::vector<Scalar> analyticD1(assembler.template nrows<firstBlock,lastBlock>());
107 analytic.write(analyticD1.begin());
108
109 // get undistorted rhs
110 Scalar const reference_value = assembler.functional();
111
112 std::vector<Scalar> xAsVector(assembler.nrows());
113 x.write(xAsVector.begin());
114
115 // reserve storage for finite difference solution
116 std::vector<Scalar> firstDerivativeError(analyticD1.size(),0.);
117
118 size_t offset = (0!=firstBlock) ? assembler.template nrows<0,firstBlock>() : 0;
119 size_t end = (lastBlock!=Functional::AnsatzVars::noOfVariables) ? (offset + assembler.template nrows<firstBlock,lastBlock>()) : xAsVector.size();
120
121 // get finite difference gradient
122 for(size_t i=offset; i<end; ++i)
123 {
124
125 xAsVector[i] += increment;
126 x.read(xAsVector.begin());
127
128 assembler.assemble(linearization(f,x),Assembler::VALUE);
129 distorted = assembler.functional();
130 xAsVector[i] -= increment;
131
132 // store finite difference
133 firstDerivativeError[i-offset] = (distorted - reference_value)/increment;
134 if(i%100==0) std::cout << ".";
135 }
136 // set x back to initial state
137 x.read(xAsVector.begin());
138 std::cout << "done." << std::endl;
139 for(size_t i=0; i<firstDerivativeError.size(); ++i)
140 firstDerivativeError[i] -= analyticD1[i];
141
142 Scalar infinityError = LinAlg::InfinityNorm()(firstDerivativeError);
143 bool errorBelowTolerance = infinityError < tolerance;
144 if(errorBelowTolerance)
145 std::cout << "\n first derivative seems to be trustworthy!\nInfinityError: " << infinityError << std::endl << std::endl;
146 else
147 std::cout << "\n first derivative does not seem to be trustworthy!\nInfinity error: " << infinityError << std::endl << std::endl;
148
149 if(toFile) vectorToVTK(firstDerivativeError);
150
151 return infinityError;
152 }
153
154 template <class Assembler, class Functional, int firstRow=0, int lastRow=Functional::TestVars::noOfVariables, int firstCol=0, int lastCol=Functional::AnsatzVars::noOfVariables>
155 typename Functional::Scalar d2(Assembler& assembler, Functional const& f, typename Functional::AnsatzVars::VariableSet const& x, typename Functional::Scalar increment = 1e-12, typename Functional::Scalar tolerance = 1e-6, bool toFile = false, std::string const& savefilename = std::string("d2error"))
156 {
157 typedef typename Functional::Scalar Scalar;
158 typedef typename Functional::TestVars::template CoefficientVectorRepresentation<firstRow,lastRow>::type RowVector;
159 typedef typename Functional::AnsatzVars::template CoefficientVectorRepresentation<firstCol,lastCol>::type ColumnVector;
160
161 std::cout << "checking second derivative:" << std::endl;
162 std::cout << "assembling...";
163
164 // sparse storage for the analytic solution
165 typedef MatrixAsTriplet<Scalar> SparseMatrix;
166 SparseMatrix secondDerivativeError;
167
168 // get reference solution
169 assembler.assemble(linearization(f,x), Assembler::RHS, Assembler::MATRIX);
170 SparseMatrix analytic_solution = assembler.template get<SparseMatrix,firstRow,lastRow,firstCol,lastCol>(false);
171 RowVector reference_value(assembler. template rhs<firstRow,lastRow>());
172 RowVector copyOfRefVal(reference_value);
173
174 // get evaluation point
175 std::vector<Scalar> xAsVector(assembler.nrows());
176 x.write(xAsVector.begin());
177
178 size_t offset = (firstRow != 0) ? assembler.template nrows<0,firstRow>() : 0;
179 size_t end = (lastRow!=Functional::TestVars::noOfVariables) ? (offset + assembler.template nrows<firstRow,lastRow>()) : xAsVector.size();
180
181 // get finite difference gradient matrix
182 for(size_t i=offset; i<end; ++i)
183 {
184 copyOfRefVal = reference_value;
185 xAsVector[i] += increment;
186 x.read(xAsVector.begin());
187
188 assembler.assemble(linearization(f,x), Assembler::RHS);
189 RowVector distorted_value(assembler.template rhs<firstRow,lastRow>());
190
191 xAsVector[i] -= increment;
192
193 std::vector<Scalar> column(end-offset);
194 copyOfRefVal -= distorted_value;
195 copyOfRefVal.write(column.begin());
196
197 secondDerivativeError.addColumn(column,i);
198 // store finite difference
199 // for(size_t k=0; k<numRows; ++k)
200 // {
201 // secondDerivativeError.addEntry(k, i, (distorted_value[k] - reference_value[k])/increment);
202 // }
203 if(i%100==0) std::cout << ".";
204 }
205 // set x back to initial state
206 x.read(xAsVector.begin());
207 std::cout << "done." << std::endl;
208 analytic_solution *= -1;
209 secondDerivativeError += analytic_solution;
210
211
212 Scalar infinityError = LinAlg::InfinityNorm()(secondDerivativeError);
213 if( (tolerance > infinityError) )
214 std::cout << "\nSecond derivative seems to be trustworthy!\nInfinity error: " << infinityError << std::endl << std::endl;
215 else
216 std::cout << "\nSecond derivative does not seem to be trustworthy!\nInfinity error: " << infinityError << std::endl << std::endl;
217
218 return infinityError;
219 }
220
221
222 template <class Vector>
223 void vectorToVTK(Vector const& vec, std::string savefilename = std::string("d1error"))
224 {
225 using std::endl;
226 using std::string;
227
228 size_t const numberOfRows = vec.size(),
229 numberOfPoints = (1+numberOfRows)*2;
230 size_t const cell_type = 8,
231 offset = 4;
232
233 XMLHelper::xml_level = 0;
234 typedef typename XMLHelper::XMLBlock XMLBlock;
235
236 savefilename.append(".vtu");
237 std::ofstream file(savefilename.c_str());
238 file << "<?xml version=\"1.0\"?>" << endl;
239 { XMLBlock xvtk(file, "VTKFile", " type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\"");
240 { XMLBlock xgrid(file, "UnstructuredGrid");
241 { string str = string(" NumberOfPoints=\"") + boost::lexical_cast<string>(numberOfPoints) + " \" numberOfRows=\"" +
242 boost::lexical_cast<string>(numberOfRows) + "\"";
243 XMLBlock xpiece(file, "Piece", str.c_str());
244 { XMLBlock xcd(file, "CellData");
245 { XMLBlock xda(file, "DataArray", " type=\"Float32\" Name=\"d2error\" NumberOfComponents=\"1\" format=\"ascii\"");
246 for(size_t cell=0; cell<numberOfRows; ++cell){
247 if(cell%12==0) file << " ";
248 // invert row order for visualization tools
249 file << fabs(vec[numberOfRows-1-cell]) << " ";
250 if(cell%12==11 || cell==(numberOfRows-1) ) file << endl;
251 }
252 }
253 }
254 { XMLBlock xpoints(file, "Points");
255 { XMLBlock xda(file, "DataArray", " type=\"Float32\" Name=\"Coordinates\" NumberOfComponents=\"3\" format=\"ascii\"");
256 for(size_t cell=0; cell<numberOfRows+1; ++cell){
257 if(cell%2==0) file << " ";
258 file << cell << " 0 0 " << cell << " 1 0 ";
259 if(cell%2==1 || cell==numberOfRows) file << endl;
260 }
261 }
262 }
263 { XMLBlock xcells(file, "Cells");
264 { XMLBlock xda(file, "DataArray", " type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\"");
265 for(size_t cell=0; cell<numberOfRows; ++cell){
266 if(cell%3==0) file << " ";
267 file << 2*cell << " " << 2*cell+1 << " " << 2*cell+2 << " " << 2*cell+3 << " " << endl;
268 if(cell%3==11 || cell==(numberOfRows-1)) file << endl;
269 }
270 }
271 { XMLBlock xda(file, "DataArray", " type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\"");
272 for(size_t cell=0; cell<numberOfRows; ++cell){
273 if(cell%12==0) file << " ";
274 file << cell*offset << " ";
275 if(cell%12==1 || cell==(numberOfRows-1)) file << endl;
276 }
277 }
278 { XMLBlock xda(file, "DataArray", " type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\"");
279 for(size_t cell=0; cell<numberOfRows; ++cell){
280 if(cell%12==0) file << " ";
281 file << cell_type << " ";
282 if(cell%12==11 || cell==(numberOfRows-1) ) file << endl;
283 }
284 }
285 }
286 }
287 }
288 }
289 }
290
291 template <class Scalar, class SparseInt>
292 void matrixToVTK(MatrixAsTriplet<Scalar,SparseInt> const& matrix, std::string savefilename = std::string("d2error"))
293 {
294 using std::cout;
295 using std::endl;
296 using std::string;
297
298 cout << "writing file";
299
301 std::vector<std::pair<SparseInt,SparseInt> > indices(matrix.data.size());
302 for(size_t i=0; i<indices.size(); ++i)
303 indices[i] = std::make_pair(matrix.ridx[i], matrix.cidx[i]);
304
305 SparseInt const numberOfRows = matrix.nrows()+1,
306 numberOfColumns = matrix.ncols()+1,
307 numberOfPoints = numberOfRows*numberOfColumns,
308 numberOfCells = (numberOfRows-1)*(numberOfColumns-1);
309 size_t const cell_type = 8, // = vtk_pixel
310 offset = 4;
311
312 XMLHelper::xml_level = 0;
313 string indent;
314 typedef typename XMLHelper::XMLBlock XMLBlock;
315
316 savefilename.append(".vtu");
317 std::ofstream file(savefilename.c_str());
318 file << "<?xml version=\"1.0\"?>" << endl;
319 { XMLBlock xvtk(file, "VTKFile", " type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\"");
320 { XMLBlock xgrid(file, "UnstructuredGrid");
321 { string str = string(" NumberOfPoints=\"") + boost::lexical_cast<string>(numberOfPoints) + " \" NumberOfCells=\"" +
322 boost::lexical_cast<string>(numberOfCells) + "\"";
323 XMLBlock xpiece(file, "Piece", str.c_str());
324 { XMLBlock xcd(file, "CellData");
325 { XMLBlock xda(file, "DataArray", " type=\"Float32\" Name=\"d2error\" NumberOfComponents=\"1\" format=\"ascii\"");
326 int cell = 0;
327 indent = XMLHelper::getXMLIndent();
328 for(SparseInt col=0; col<numberOfColumns-1; ++col){
329 for(SparseInt row=0; row<numberOfRows-1; ++row){
330 if(cell%12==0) file << indent.c_str();
331 typename std::vector<std::pair<SparseInt,SparseInt> >::const_iterator iter = std::find(indices.begin(), indices.end(), std::make_pair(numberOfRows-2-row,col)),
332 iterStart = indices.begin();
333 if(iter!=indices.end())
334 file << fabs(matrix.data[std::distance(iterStart, iter)]) << " ";
335 else
336 file << 0.0 << " ";
337 if(cell%12==11 || cell==numberOfCells-1) file << endl;
338 ++cell;
339 }
340 }
341 }
342 }
343 cout << ".";
344 { XMLBlock xp(file, "Points");
345 { XMLBlock xda(file, "DataArray", " type=\"Float32\" Name=\"Coordinates\" NumberOfComponents=\"3\" format=\"ascii\"");
346 int cell = 0;
347 for(SparseInt col=0; col<numberOfColumns; ++col){
348 for(SparseInt row=0; row<numberOfRows; ++row){
349 if(cell%4==0) file << indent.c_str();
350 file << row << " " << col << " 0 ";
351 if(cell%4==3 || cell==(numberOfCells-1)) file << endl;
352 ++cell;
353 }
354 }
355 }
356 }
357 cout << ".";
358 { XMLBlock xcell(file, "Cells");
359 { XMLBlock xda(file, "DataArray", " type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\"");
360 int cell = 0;
361 for(SparseInt col=0; col<numberOfColumns-1; ++col){
362 for(SparseInt row=0; row<numberOfRows-1; ++row){
363 if(cell%3==0) file << indent.c_str();
364 Private::writeConnectivity(file, row, col, numberOfRows);
365 if(cell%3==2 || cell==(numberOfCells-1)) file << endl;
366 ++cell;
367 }
368 }
369 }
370 cout << ".";
371 { XMLBlock xda(file, "DataArray", " type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\"");
372 for(SparseInt cell=0; cell<numberOfCells; ++cell){
373 if(cell%12==0) file << indent.c_str();
374 file << (cell+1)*offset << " ";
375 if(cell%12==11 || cell==(numberOfCells-1)) file << endl;
376 }
377 }
378 cout << ".";
379 { XMLBlock xda(file, "DataArray", " type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\"");
380 for(SparseInt cell=0; cell<numberOfCells; ++cell){
381 if(cell%12==0) file << indent.c_str();
382 file << cell_type << " ";
383 if(cell%12==11 || cell==(numberOfCells-1) ) file << endl;
384 }
385 }
386 }
387 }
388 }
389 }
390 cout << "done" << endl;
391 }
392
393} // end of namespace DerivativeCheck
394
396
406template <class Functional, bool checkD1=true, class SparseInt=int>
407class DerivativeChecker
408{
409 typedef VariationalFunctionalAssembler<LinearizationAt<Functional> > Assembler;
410 typedef typename Functional::AnsatzVars VariableSet;
411 typedef typename VariableSet::VariableSet EvaluationPoint;
412 typedef typename VariableSet::template CoefficientVectorRepresentation<> CoefficientVectorRepresentation;
413 typedef typename CoefficientVectorRepresentation::type CoefficientVector;
414
415public:
416 typedef typename Functional::Scalar Scalar;
417 enum{ idOfFirstBlock = 0,
418 idOfLastBlock = boost::mpl::size<typename Functional::AnsatzVars::Variables>::type::value };
419
421
427 DerivativeChecker(Assembler& assembler_, Functional const& f_, EvaluationPoint const& x_, VariableSet const& variableSet_, Scalar const tol = 1.0e-6, Scalar const incr = 1.0e-9) :
428 assembler(assembler_), f(f_), x(x_), variableSet(variableSet_), tolerance(tol), increment(incr)
429 {
430 if(checkD1) checkFirstDerivative();
431 checkSecondDerivative();
432 }
433
434
436
440 void checkFirstDerivative()
441 {
442 std::cout << "checking first derivative" << std::endl;
443 std::cout << "assembling...";
444
445 // storage for the evaluation of the derivative at x and the at x+increment*e_i
446 Scalar distorted;
447
448 // get reference solution
449 assembler.assemble(linearization(f,x), Assembler::RHS | Assembler::VALUE );
450
451 CoefficientVector analytic(assembler.rhs());
452 std::vector<Scalar> analyticD1(assembler.nrows());
453 analytic.write(analyticD1.begin());
454
455 // get undistorted rhs
456 Scalar const reference_value = assembler.functional();
457
458 std::vector<Scalar> xAsVector(analyticD1.size());
459 x.write(xAsVector.begin());
460
461 // reserve storage for finite difference solution
462 firstDerivativeError.resize(xAsVector.size(),0);
463
464 EvaluationPoint copyOfX(x);
465 // get finite difference gradient
466 for(size_t i=0; i<xAsVector.size(); ++i)
467 {
468 xAsVector[i] += increment;
469 copyOfX.read(xAsVector.begin());
470
471 assembler.assemble(linearization(f,copyOfX),Assembler::VALUE);
472 Scalar distorted_value = assembler.functional();
473 xAsVector[i] -= increment;
474
475 // store finite difference
476 firstDerivativeError[i] = (distorted_value - reference_value)/increment;
477 if(i%100==0) std::cout << ".";
478 }
479 std::cout << "done." << std::endl;
480 for(size_t i=0; i<firstDerivativeError.size(); ++i)
481 firstDerivativeError[i] -= analyticD1[i];
482
483 Scalar infinityError = LinAlg::InfinityNorm()(firstDerivativeError);
484 firstDerivativeOk_ = infinityError < tolerance;
485 if(firstDerivativeOk_)
486 std::cout << "\n first derivative seems to be trustworthy!\nInfinityError: " << infinityError << std::endl << std::endl;
487 else
488 std::cout << "\n first derivative does not seem to be trustworthy!\nInfinity error: " << infinityError << std::endl << std::endl;
489 }
490
491
493
497 void checkSecondDerivative()
498 {
499 std::cout << "checking second derivative:" << std::endl;
500 // Define linearization at x
501 std::cout << "assembling...";
502
503 // sparse storage for the analytic solution
504 typedef MatrixAsTriplet<Scalar,SparseInt> SparseMatrix;
505
506 // get reference solution
507 assembler.assemble(linearization(f,x), Assembler::RHS, Assembler::MATRIX);
508 SparseMatrix analytic_solution = assembler.template get<SparseMatrix>(false);
509 CoefficientVector reference_value(assembler.rhs());
510 CoefficientVector copyOfRefVal(reference_value);
511
512 // get evaluation point
513 std::vector<Scalar> xAsVector(assembler.nrows());
514 x.write(xAsVector.begin());
515
516 EvaluationPoint copyOfX(x);
517 // get finite difference gradient matrix
518 for(size_t i=0; i<xAsVector.size(); ++i)
519 {
520 copyOfRefVal = reference_value;
521 xAsVector[i] += increment;
522 copyOfX.read(xAsVector.begin());
523
524 assembler.assemble(linearization(f,copyOfX), Assembler::RHS);
525 CoefficientVector distorted_value(assembler.rhs());
526
527 xAsVector[i] -= increment;
528
529 std::vector<Scalar> column(xAsVector.size());
530 copyOfRefVal -= distorted_value;
531 copyOfRefVal.write(column.begin());
532
533 secondDerivativeError.addColumn(column,i);
534 // store finite difference
535// for(size_t k=0; k<numRows; ++k)
536// {
537// secondDerivativeError.addEntry(k, i, (distorted_value[k] - reference_value[k])/increment);
538// }
539 if(i%100==0) std::cout << ".";
540 }
541 std::cout << "done." << std::endl;
542 analytic_solution *= -1;
543 secondDerivativeError += analytic_solution;
544
545
546 Scalar infinityError = LinAlg::InfinityNorm()(secondDerivativeError);
547 if( (secondDerivativeOk_ = tolerance > infinityError) )
548 std::cout << "\nSecond derivative seems to be trustworthy!\nInfinity error: " << infinityError << std::endl << std::endl;
549 else
550 std::cout << "\nSecond derivative does not seem to be trustworthy!\nInfinity error: " << infinityError << std::endl << std::endl;
551 }
552
553 bool firstDerivativeOk() const { return firstDerivativeOk_; }
554
555 bool secondDerivativeOk() const { return secondDerivativeOk_; }
556
557// /// get block of d1Error
558// std::vector<Scalar> getBlocksD1Error(int const firstBlockId, int const lastBlockId) const
559// {
560// assert(firstBlockId < lastBlockId);
561// std::pair<int,int> ids = getRowIds(firstBlockId, lastBlockId);
562// std::vector<Scalar> result;
563// for(int row = ids.first; row<ids.second; ++row)
564// result.push_back(firstDerivativeError[row]);
565// return result;
566// }
567//
568// /// get block of d2Error
569// MatrixAsTriplet<Scalar> getBlocksD2Error(int const firstBlockId, int const lastBlockId) const
570// {
571// assert(firstBlockId < lastBlockId);
572// std::pair<int,int> ids = getRowIds(firstBlockId, lastBlockId);
573//
574// MatrixAsTriplet<Scalar> blocks;
575// for(int i=0; i<secondDerivativeError.size(); ++i)
576// {
577// if(secondDerivativeError.ridx[i] >= ids.first && secondDerivativeError.ridx[i] < ids.second &&
578// secondDerivativeError.cidx[i] >= ids.first && secondDerivativeError.cidx[i] < ids.second)
579// blocks.addEntry(secondDerivativeError.ridx[i], secondDerivativeError.cidx[i], secondDerivativeError.data[i]);
580// }
581// return blocks;
582// }
583
585 void setTolerance(Scalar const tol){ tolerance = tol; }
586
588 void setIncrement(Scalar const incr){ increment = incr; }
589
591 void printD1Error() const
592 {
593 std::cout << "components in D1 over tolerance:\n";
594 for(size_t i=0; i<firstDerivativeError.size(); ++i)
595 if(fabs(firstDerivativeError[i]) > tolerance)
596 std::cout << i << ": " << firstDerivativeError[i] << std::endl;
597 }
598
600 void printFullD1Error() const
601 {
602 std::cout << "components in D1:\n" << firstDerivativeError << std::endl;
603 }
604
606 void printD2Error() const
607 {
608 std::cout << "components in D2 over tolerance:\n";
609 secondDerivativeError.print(std::cout,tolerance);
610 }
611
613 void printFullD2Error() const
614 {
615 std::cout << "components in D2:\n";
616 secondDerivativeError.print(std::cout);
617 }
618
620 void d2ToVTK(){
621 assembler.assemble(linearization(f,x),Assembler::MATRIX);
622 // sparse storage for the analytic solution
623 MatrixAsTriplet<Scalar,SparseInt> analytic_solution = assembler.template get<MatrixAsTriplet<Scalar,SparseInt>>(false);
624 // get reference solution
625 matrixToVTK(analytic_solution);
626 }
627
629 void d1ErrorToVTK(int const steps) const { vectorToVTK(firstDerivativeError, steps); }
630
631// /// write block of error in first derivative to .vtu file
632// void d1BlockErrorToVTK(int const firstBlockId, int const lastBlockId) const { vectorToVTK(getBlocksD1Error(firstBlockId,lastBlockId)); }
633
635
638 template <class Vector>
639 void vectorToVTK(Vector &vector, int const steps=0) const {
640
641 using std::endl;
642 using std::string;
643
644 size_t const numberOfRows = vector.size(),
645 numberOfPoints = (1+numberOfRows)*2;
646 size_t const cell_type = 8,
647 offset = 4;
648
649 XMLHelper::xml_level = 0;
650 typedef typename XMLHelper::XMLBlock XMLBlock;
651
652 string filename("d1error");
653 filename.append(boost::lexical_cast<string>(steps)).append(".vtu");
654
655 std::ofstream file(filename.c_str());
656 file << "<?xml version=\"1.0\"?>" << endl;
657 { XMLBlock xvtk(file, "VTKFile", " type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\"");
658 { XMLBlock xgrid(file, "UnstructuredGrid");
659 { string str = string(" NumberOfPoints=\"") + boost::lexical_cast<string>(numberOfPoints) + " \" numberOfRows=\"" +
660 boost::lexical_cast<string>(numberOfRows) + "\"";
661 XMLBlock xpiece(file, "Piece", str.c_str());
662 { XMLBlock xcd(file, "CellData");
663 { XMLBlock xda(file, "DataArray", " type=\"Float32\" Name=\"d2error\" NumberOfComponents=\"1\" format=\"ascii\"");
664 for(size_t cell=0; cell<numberOfRows; ++cell){
665 if(cell%12==0) file << " ";
666 // invert row order for visualization tools
667 file << fabs(vector[numberOfRows-1-cell]) << " ";
668 if(cell%12==11 || cell==(numberOfRows-1) ) file << endl;
669 }
670 }
671 }
672 { XMLBlock xpoints(file, "Points");
673 { XMLBlock xda(file, "DataArray", " type=\"Float32\" Name=\"Coordinates\" NumberOfComponents=\"3\" format=\"ascii\"");
674 for(size_t cell=0; cell<numberOfRows+1; ++cell){
675 if(cell%2==0) file << " ";
676 file << cell << " 0 0 " << cell << " 1 0 ";
677 if(cell%2==1 || cell==numberOfRows) file << endl;
678 }
679 }
680 }
681 { XMLBlock xcells(file, "Cells");
682 { XMLBlock xda(file, "DataArray", " type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\"");
683 for(size_t cell=0; cell<numberOfRows; ++cell){
684 if(cell%3==0) file << " ";
685 file << 2*cell << " " << 2*cell+1 << " " << 2*cell+2 << " " << 2*cell+3 << " " << endl;
686 if(cell%3==11 || cell==(numberOfRows-1)) file << endl;
687 }
688 }
689 { XMLBlock xda(file, "DataArray", " type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\"");
690 for(size_t cell=0; cell<numberOfRows; ++cell){
691 if(cell%12==0) file << " ";
692 file << cell*offset << " ";
693 if(cell%12==1 || cell==(numberOfRows-1)) file << endl;
694 }
695 }
696 { XMLBlock xda(file, "DataArray", " type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\"");
697 for(size_t cell=0; cell<numberOfRows; ++cell){
698 if(cell%12==0) file << " ";
699 file << cell_type << " ";
700 if(cell%12==11 || cell==(numberOfRows-1) ) file << endl;
701 }
702 }
703 }
704 }
705 }
706 }
707 }
708
710 void d2ErrorToVTK(int const steps) const { matrixToVTK(secondDerivativeError,steps); }
711
712// /// write block of second derivative to .vtu-file
713// void d2BlockErrorToVTK(int firstBlockId, int lastBlockId) const { matrixToVTK(getBlocksD2Error(firstBlockId,lastBlockId)); }
714
716 void matrixToVTK(MatrixAsTriplet<Scalar,SparseInt> const& matrix, int const steps=0) const
717 {
718 using std::cout;
719 using std::endl;
720 using std::string;
721
722 cout << "writing file";
723
725 std::vector<std::pair<SparseInt,SparseInt> > indices(matrix.data.size());
726 for(size_t i=0; i<indices.size(); ++i)
727 indices[i] = std::make_pair(matrix.ridx[i], matrix.cidx[i]);
728
729 SparseInt const numberOfRows = matrix.nrows()+1,
730 numberOfColumns = matrix.ncols()+1,
731 numberOfPoints = numberOfRows*numberOfColumns,
732 numberOfCells = (numberOfRows-1)*(numberOfColumns-1);
733 size_t const cell_type = 8, // = vtk_pixel
734 offset = 4;
735
736 XMLHelper::xml_level = 0;
737 string indent;
738 typedef typename XMLHelper::XMLBlock XMLBlock;
739
740 string filename("d2error");
741 filename.append(boost::lexical_cast<string>(steps)).append(".vtu");
742
743 std::ofstream file(filename.c_str());
744 file << "<?xml version=\"1.0\"?>" << endl;
745 { XMLBlock xvtk(file, "VTKFile", " type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\"");
746 { XMLBlock xgrid(file, "UnstructuredGrid");
747 { string str = string(" NumberOfPoints=\"") + boost::lexical_cast<string>(numberOfPoints) + " \" NumberOfCells=\"" +
748 boost::lexical_cast<string>(numberOfCells) + "\"";
749 XMLBlock xpiece(file, "Piece", str.c_str());
750 { XMLBlock xcd(file, "CellData");
751 { XMLBlock xda(file, "DataArray", " type=\"Float32\" Name=\"d2error\" NumberOfComponents=\"1\" format=\"ascii\"");
752 int cell = 0;
753 indent = XMLHelper::getXMLIndent();
754 for(SparseInt col=0; col<numberOfColumns-1; ++col){
755 for(SparseInt row=0; row<numberOfRows-1; ++row){
756 if(cell%12==0) file << indent.c_str();
757 typename std::vector<std::pair<SparseInt,SparseInt> >::const_iterator iter = std::find(indices.begin(), indices.end(), std::make_pair(numberOfRows-2-row,col)),
758 iterStart = indices.begin();
759 if(iter!=indices.end())
760 file << fabs(matrix.data[std::distance(iterStart, iter)]) << " ";
761 if(cell%12==11 || cell==numberOfCells-1) file << endl;
762 ++cell;
763 }
764 }
765 }
766 }
767 cout << ".";
768 { XMLBlock xp(file, "Points");
769 { XMLBlock xda(file, "DataArray", " type=\"Float32\" Name=\"Coordinates\" NumberOfComponents=\"3\" format=\"ascii\"");
770 int cell = 0;
771 for(size_t col=0; col<numberOfColumns; ++col){
772 for(size_t row=0; row<numberOfRows; ++row){
773 if(cell%4==0) file << indent.c_str();
774 file << row << " " << col << " 0 ";
775 if(cell%4==3 || cell==(numberOfCells-1)) file << endl;
776 ++cell;
777 }
778 }
779 }
780 }
781 cout << ".";
782 { XMLBlock xcell(file, "Cells");
783 { XMLBlock xda(file, "DataArray", " type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\"");
784 int cell = 0;
785 for(size_t col=0; col<numberOfColumns-1; ++col){
786 for(size_t row=0; row<numberOfRows-1; ++row){
787 if(cell%3==0) file << indent.c_str();
788 writeConnectivity(file, row, col, numberOfRows);
789 if(cell%3==2 || cell==(numberOfCells-1)) file << endl;
790 ++cell;
791 }
792 }
793 }
794 cout << ".";
795 { XMLBlock xda(file, "DataArray", " type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\"");
796 for(size_t cell=0; cell<numberOfCells; ++cell){
797 if(cell%12==0) file << indent.c_str();
798 file << (cell+1)*offset << " ";
799 if(cell%12==11 || cell==(numberOfCells-1)) file << endl;
800 }
801 }
802 cout << ".";
803 { XMLBlock xda(file, "DataArray", " type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\"");
804 for(size_t cell=0; cell<numberOfCells; ++cell){
805 if(cell%12==0) file << indent.c_str();
806 file << cell_type << " ";
807 if(cell%12==11 || cell==(numberOfCells-1) ) file << endl;
808 }
809 }
810 }
811 }
812 }
813 }
814 cout << "done" << endl;
815 }
816
817private:
818// /// get ids of first and last row in [firstBlockId,lastBlockId]
819// std::pair<int,int> getRowIds(int const firstBlockId, int const lastBlockId) const
820// {
821// int const firstIdx = (firstBlockId > 0) ? linearization.rows(0,firstBlockId) : 0;
822// int const lastIdx = firstIdx + linearization.rows(firstBlockId, lastBlockId);
823// return std::make_pair(firstIdx, lastIdx);
824// }
825
826 void writeConnectivity(std::ofstream &file, SparseInt row, SparseInt col, SparseInt numRows) const
827 {
828 int const upperLeftVertexId = row+col*numRows;
829 file << upperLeftVertexId << " " << upperLeftVertexId+1 << " " << upperLeftVertexId+numRows << " " << upperLeftVertexId+numRows+1 << " ";
830 }
831
832 Assembler& assembler;
833 Functional const& f;
834 EvaluationPoint x;
835 VariableSet const& variableSet;
836 Scalar tolerance, increment;
837 std::vector<Scalar> firstDerivativeError;
838 MatrixAsTriplet<Scalar,SparseInt> secondDerivativeError;
839 bool firstDerivativeOk_, secondDerivativeOk_;
840};
841} // end of namespace Kaskade
842#endif /* CHECK_DERIVATIVE_HH */
Kaskade::DerivativeChecker
Class that checks the derivatives of a functional at a linearization point.
Definition: check_derivative.hh:408
Kaskade::DerivativeChecker::DerivativeChecker
DerivativeChecker(Assembler &assembler_, Functional const &f_, EvaluationPoint const &x_, VariableSet const &variableSet_, Scalar const tol=1.0e-6, Scalar const incr=1.0e-9)
Constructor.
Definition: check_derivative.hh:427
Kaskade::DerivativeChecker::checkFirstDerivative
void checkFirstDerivative()
Check first derivative.
Definition: check_derivative.hh:440
Kaskade::DerivativeChecker::printD2Error
void printD2Error() const
print error in second derivative (only values below tolerance)
Definition: check_derivative.hh:606
Kaskade::DerivativeChecker::printD1Error
void printD1Error() const
print error in first derivative (only values below tolerance)
Definition: check_derivative.hh:591
Kaskade::DerivativeChecker::setTolerance
void setTolerance(Scalar const tol)
Only errors below tolerance will be displayed.
Definition: check_derivative.hh:585
Kaskade::DerivativeChecker::setIncrement
void setIncrement(Scalar const incr)
Set increment for finite differences.
Definition: check_derivative.hh:588
Kaskade::DerivativeChecker::printFullD2Error
void printFullD2Error() const
print error in second derivative
Definition: check_derivative.hh:613
Kaskade::DerivativeChecker::matrixToVTK
void matrixToVTK(MatrixAsTriplet< Scalar, SparseInt > const &matrix, int const steps=0) const
write, possibly sparse, matrix in triplet format to .vtu-file
Definition: check_derivative.hh:716
Kaskade::DerivativeChecker::d2ToVTK
void d2ToVTK()
write d2 to .vtu file
Definition: check_derivative.hh:620
Kaskade::DerivativeChecker::checkSecondDerivative
void checkSecondDerivative()
Check second derivative.
Definition: check_derivative.hh:497
Kaskade::DerivativeChecker::d2ErrorToVTK
void d2ErrorToVTK(int const steps) const
write error in second derivative to .vtu-file
Definition: check_derivative.hh:710
Kaskade::DerivativeChecker::printFullD1Error
void printFullD1Error() const
print error in first derivative
Definition: check_derivative.hh:600
Kaskade::DerivativeChecker::d1ErrorToVTK
void d1ErrorToVTK(int const steps) const
write error in first derivative to .vtu file
Definition: check_derivative.hh:629
Kaskade::DerivativeChecker::vectorToVTK
void vectorToVTK(Vector &vector, int const steps=0) const
write Vector vector to .vtu file
Definition: check_derivative.hh:639
Kaskade::MatrixAsTriplet::nrows
SparseIndexInt nrows() const
Returns number of rows (computes them, if not known)
Definition: triplet.hh:202
Kaskade::MatrixAsTriplet::ncols
SparseIndexInt ncols() const
Returns number of cols (computes them, if not known)
Definition: triplet.hh:215
Kaskade::MatrixAsTriplet::ridx
std::vector< SparseIndexInt > ridx
row indices
Definition: triplet.hh:669
Kaskade::MatrixAsTriplet::cidx
std::vector< SparseIndexInt > cidx
column indices
Definition: triplet.hh:671
Kaskade::MatrixAsTriplet::data
std::vector< Scalar > data
data
Definition: triplet.hh:673
Kaskade::VariableSet::VariableSet
VariableSet(Self const &vs)
Copy constructor.
Definition: variables.hh:362
Kaskade::VariationalFunctionalAssembler
A class for assembling Galerkin representation matrices and right hand sides for variational function...
Definition: assemblerCore.hh:1432
Kaskade::VariationalFunctionalAssembler::VALUE
static unsigned int const VALUE
DEPRECATED, enums in the Assembler namespace.
Definition: assemblerCore.hh:1528
dune_bridge.hh
Kaskade::Assembler::VALUE
@ VALUE
Definition: assemblerCore.hh:1402
Kaskade::Assembler::RHS
@ RHS
Definition: assemblerCore.hh:1402
Kaskade::Assembler::MATRIX
@ MATRIX
Definition: assemblerCore.hh:1402
GeometricObject::distance
Scalar distance(Point< Scalar, dim > const &first, Point< Scalar, dim > const &second)
Kaskade::DerivativeCheck::Private::writeConnectivity
void writeConnectivity(std::ofstream &file, SparseInt row, SparseInt col, SparseInt numRows)
Definition: check_derivative.hh:83
Kaskade::DerivativeCheck::d1
Functional::Scalar d1(Assembler &assembler, Functional const &f, typename Functional::AnsatzVars::VariableSet &x, typename Functional::Scalar tolerance=1e-6, typename Functional::Scalar increment=1e-12, bool toFile=false, std::string const &filename=std::string("d1error"))
Definition: check_derivative.hh:91
Kaskade::DerivativeCheck::matrixToVTK
void matrixToVTK(MatrixAsTriplet< Scalar, SparseInt > const &matrix, std::string savefilename=std::string("d2error"))
Definition: check_derivative.hh:292
Kaskade::DerivativeCheck::vectorToVTK
void vectorToVTK(Vector const &vec, std::string savefilename=std::string("d1error"))
Definition: check_derivative.hh:223
Kaskade::DerivativeCheck::d2
Functional::Scalar d2(Assembler &assembler, Functional const &f, typename Functional::AnsatzVars::VariableSet const &x, typename Functional::Scalar increment=1e-12, typename Functional::Scalar tolerance=1e-6, bool toFile=false, std::string const &savefilename=std::string("d2error"))
Definition: check_derivative.hh:155
Kaskade::Simplex::Vector
Dune::FieldVector< Scalar, dim > Vector
Definition: boundaryInterpolation.hh:176
Kaskade::XMLHelper::getXMLIndent
std::string getXMLIndent()
Definition: check_derivative.hh:36
Kaskade::column
Dune::FieldVector< Scalar, n > column(Dune::FieldMatrix< Scalar, n, m > const &A, int j)
Extracts the j-th column of the given matrix A.
Definition: fixdune.hh:445
Kaskade::LinAlg::InfinityNorm
Infinity norm for stl-container, Dune::FieldVector and Dune::FieldMatrix.
Definition: scalarproducts.hh:113
Kaskade::XMLHelper::XMLBlock::XMLBlock
XMLBlock(std::ofstream &f, std::string &s, std::string d="")
Definition: check_derivative.hh:50
Kaskade::XMLHelper::XMLBlock::XMLBlock
XMLBlock(std::ofstream &f, const char *s, const char *d="")
Definition: check_derivative.hh:44
Generated on Fri Sep 6 2024 13:13:00 for KASKADE 7 development version by doxygen 1.9.4