1 // include configurations options
2
3 #ifdef HAVE_CONFIG_H
4 #include <config.h>
5 #endif
6
7 // Includes from the IOStream Library
8 // ----------------------------------
9
10 #include <iostream>
11 #include <sstream>
12
13 //#define USE_OLD_RANNACHERTUREK_SPACE
14
15 // Includes from DUNE-FEM
16 // ----------------------
17
18 // include Lagrange discrete function space
19 #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
20 #include <dune/fem/space/lagrange.hh>
21 #include <dune/fem/space/common/adaptationmanager.hh>
22 #include <dune/fem/function/localfunction/bindable.hh>
23
24 #if HAVE_PETSC && defined USE_PETSCDISCRETEFUNCTION
25 #include <dune/fem/function/petscdiscretefunction.hh>
26 #include <dune/fem/operator/linear/petscoperator.hh>
27 #include <dune/fem/solver/petscinverseoperators.hh>
28 #elif HAVE_DUNE_ISTL && defined USE_BLOCKVECTORFUNCTION
29 #include <dune/fem/function/blockvectorfunction.hh>
30 #include <dune/fem/operator/linear/istloperator.hh>
31 #include <dune/fem/solver/istlinverseoperators.hh>
32 #else
33 #include <dune/fem/function/adaptivefunction.hh>
34 #include <dune/fem/operator/linear/spoperator.hh>
35 #if HAVE_SUITESPARSE_UMFPACK && defined USE_UMFPACK
36 #include <dune/fem/solver/umfpacksolver.hh>
37 #else
38 #include <dune/fem/solver/krylovinverseoperators.hh>
39 #endif
40 #endif
41
42 #include <dune/fem/misc/l2norm.hh>
43 #include <dune/fem/misc/h1norm.hh>
44
45 #include "massoperator.hh"
46 #include "testgrid.hh"
47
48 // Global Type Definitions
49 // -----------------------
50
51 typedef Dune::GridSelector :: GridType GridType;
52
53 typedef Dune::Fem::AdaptiveLeafGridPart< GridType, Dune::InteriorBorder_Partition > GridPartType;
54 typedef Dune::Fem::FunctionSpace< typename GridType::ctype, typename GridType::ctype, GridType::dimensionworld, 1 > SpaceType;
55 typedef Dune::Fem::DynamicLagrangeDiscreteFunctionSpace< SpaceType, GridPartType > DiscreteSpaceType;
56
57 #if HAVE_PETSC && defined USE_PETSCDISCRETEFUNCTION
58 typedef Dune::Fem::PetscDiscreteFunction< DiscreteSpaceType > DiscreteFunctionType;
59 typedef Dune::Fem::PetscLinearOperator< DiscreteFunctionType, DiscreteFunctionType > LinearOperatorType;
60 typedef Dune::Fem::PetscInverseOperator< DiscreteFunctionType, LinearOperatorType > InverseOperatorType;
61 #elif HAVE_DUNE_ISTL && defined USE_BLOCKVECTORFUNCTION
62 typedef Dune::Fem::ISTLBlockVectorDiscreteFunction< DiscreteSpaceType > DiscreteFunctionType;
63 typedef Dune::Fem::ISTLLinearOperator< DiscreteFunctionType, DiscreteFunctionType > LinearOperatorType;
64 typedef Dune::Fem::ISTLInverseOperator< DiscreteFunctionType, Dune::Fem::SolverParameter::cg > InverseOperatorType;
65 #else
66 typedef Dune::Fem::AdaptiveDiscreteFunction< DiscreteSpaceType > DiscreteFunctionType;
67 typedef Dune::Fem::SparseRowLinearOperator< DiscreteFunctionType, DiscreteFunctionType > LinearOperatorType;
68 #if HAVE_SUITESPARSE_UMFPACK && defined USE_UMFPACK
69 typedef Dune::Fem::UMFPACKInverseOperator< DiscreteFunctionType > InverseOperatorType;
70 #else
71 typedef Dune::Fem::KrylovInverseOperator< DiscreteFunctionType > InverseOperatorType;
72 #endif
73 #endif
74
75 typedef MassOperator< DiscreteFunctionType, LinearOperatorType > MassOperatorType;
76
77 // Function to Project
78 // -------------------
79 template <class GridPart, class RangeType>
80 struct Function : public Dune::Fem::BindableGridFunction< GridPart, RangeType >
81 {
82 typedef Dune::Fem::BindableGridFunction<GridPart, RangeType > Base;
83 using Base::Base;
84
85 template <class Point>
evaluateFunction86 void evaluate ( const Point &p, RangeType &value ) const
87 {
88 auto x = Base::global(p);
89 value = 1.;
90 for( int k = 0; k < x.dimension; ++k )
91 value *= sin( M_PI * x[k] );
92 }
93 template <class Point>
jacobianFunction94 void jacobian ( const Point &p, typename Base::JacobianRangeType &jacobian ) const
95 {
96 auto x = Base::global(p);
97 for( int j = 0; j < x.dimension; ++j )
98 {
99 // jacobian has only one row, calc j-th column
100 jacobian[0][j] = M_PI;
101 for( int k = 0; k < x.dimension; ++k )
102 jacobian[0][j] *= (j == k ? cos( M_PI*x[k] ) : sin( M_PI*x[k] ));
103 }
104 }
orderFunction105 unsigned int order() const { return 4; }
nameFunction106 std::string name() const { return "MyFunction"; } // needed for output
107 };
108
109 // Algorithm
110 // ---------
111
112 struct Algorithm
113 {
114 typedef Dune::FieldVector< double, 2 > ErrorType;
115 typedef Function< GridPartType, SpaceType::RangeType > FunctionType;
116
117 explicit Algorithm ( GridType &grid );
118 ErrorType operator() ( int step, int polOrder );
119 ErrorType finalize ( DiscreteFunctionType &u );
120 DiscreteSpaceType &space ();
121 void nextMesh ();
122 void resetMesh (const int steps);
123
124 private:
125 GridPartType gridPart_;
126 FunctionType function_;
127 };
128
Algorithm(GridType & grid)129 inline Algorithm::Algorithm ( GridType &grid )
130 : gridPart_( grid ), function_(gridPart_)
131 {
132 }
133
nextMesh()134 inline void Algorithm :: nextMesh ()
135 {
136 Dune::Fem::GlobalRefine::apply(gridPart_.grid(), Dune::DGFGridInfo< GridType >::refineStepsForHalf() );
137 }
138
resetMesh(const int steps)139 inline void Algorithm :: resetMesh (const int steps)
140 {
141 gridPart_.grid().globalRefine( -steps * Dune::DGFGridInfo< GridType >::refineStepsForHalf() );
142 }
143
operator ()(int step,int polOrder)144 inline Algorithm::ErrorType Algorithm::operator() ( int step, int polOrder )
145 {
146 DiscreteSpaceType space( gridPart_, polOrder );
147 DiscreteFunctionType solution( "solution", space );
148
149 // get operator
150 MassOperatorType massOperator( space );
151
152 // assemble RHS
153 DiscreteFunctionType rhs( "rhs", space );
154 massOperator.assembleRHS( function_, rhs );
155
156 unsigned long maxIter = space.size();
157 maxIter = space.gridPart().comm().sum( maxIter );
158
159 // clear result
160 solution.clear();
161
162 // apply solver
163 InverseOperatorType inverseOperator;
164 inverseOperator.bind( massOperator );
165
166 inverseOperator.parameter().setTolerance( 1e-10 );
167 inverseOperator.parameter().setMaxIterations( maxIter );
168
169 inverseOperator( rhs, solution );
170 return finalize( solution );
171 }
172
173
finalize(DiscreteFunctionType & solution)174 inline Algorithm::ErrorType Algorithm::finalize ( DiscreteFunctionType &solution )
175 {
176 const GridPartType &gridPart = solution.space().gridPart();
177 ErrorType error;
178
179 Dune::Fem::L2Norm< GridPartType > l2norm( gridPart );
180 Dune::Fem::H1Norm< GridPartType > h1norm( gridPart );
181
182 error[ 0 ] = l2norm.distance( function_, solution );
183 error[ 1 ] = h1norm.distance( function_, solution );
184 std::cout << error <<std::endl;
185 return error;
186 }
187
run(GridType & grid,const int polOrder)188 void run( GridType& grid, const int polOrder )
189 {
190 const int nrSteps = 4;
191
192 std::cout<< "testing with polorder "<< polOrder <<std::endl;
193 Algorithm algorithm( grid );
194 std::vector< typename Algorithm::ErrorType > error( nrSteps );
195 for( int step = 0; step<nrSteps; ++step )
196 {
197 error[ step ] = algorithm( step, polOrder );
198 algorithm.nextMesh();
199 }
200
201 for( int step = 1; step <nrSteps; ++step )
202 {
203 double l2eoc = log( error[ step ][ 0 ] / error[ step -1 ][ 0 ] ) / log( 0.5 );
204 double h1eoc = log( error[ step ][ 1 ] / error[ step -1 ][ 1 ] ) / log( 0.5 );
205
206 if( std::abs( l2eoc -1 - polOrder ) > 0.2 )
207 {
208 DUNE_THROW(Dune::InvalidStateException,"EOC check of solving mass matrix system failed L2eoc " << l2eoc << " " << polOrder);
209 }
210
211 if( std::abs( h1eoc - polOrder ) > 0.2 )
212 {
213 //note: This will fail with Yaspgrid, bug in Geometry JacobianInverse
214 DUNE_THROW(Dune::InvalidStateException,"EOC check of solving mass matrix system failed H1eoc " << h1eoc << " " << polOrder);
215 }
216 }
217
218 algorithm.resetMesh(nrSteps);
219 }
220
221
222 // Main Program
223 // ------------
224
main(int argc,char ** argv)225 int main ( int argc, char **argv )
226 try
227 {
228 // initialize MPI manager and PETSc
229 Dune::Fem::MPIManager::initialize( argc, argv );
230
231 // add command line parameters to global parameter table
232 Dune::Fem::Parameter::append( argc, argv );
233 // append parameters from the parameter file
234 Dune::Fem::Parameter::append( (argc < 2) ? "parameter" : argv[ 1 ] );
235
236 GridType &grid = Dune::Fem::TestGrid :: grid();
237 for( int p=1; p<2; ++p )
238 run( grid, p );
239
240 Dune::Fem::Parameter::write( "parameter.log" );
241
242 return 0;
243 }
244 catch( const Dune::Exception &exception )
245 {
246 // display the exception message on the console
247 std::cerr << exception << std::endl;
248 return 1;
249 }
250