Commit 3dd68990 authored by Steffen Müthing's avatar Steffen Müthing
Browse files

[stationary-poisson] Update to work with and without instrumentation

parent ae56d958
......@@ -98,6 +98,55 @@ stationary_poisson_blocked_old_LDFLAGS = \
$(DUNE_LDFLAGS)
noinst_PROGRAMS += stationary-poisson-blocked-counted
stationary_poisson_blocked_counted_SOURCES = stationary-poisson.cc
stationary_poisson_blocked_counted_CPPFLAGS = \
$(AM_CPPFLAGS) \
$(TBB_CPPFLAGS) \
$(DUNEMPICPPFLAGS) \
-DSIMD_BLOCK_SIZE=16 \
-DBLOCK_SOLVER=1 \
-DLOP_TYPE=sumfactorized \
-DENABLE_HP_TIMERS \
-DENABLE_COUNTER
stationary_poisson_blocked_counted_LDADD = \
$(DUNE_LDFLAGS) $(DUNE_LIBS) \
$(DUNEMPILIBS) \
$(TBB_LIBS) $(TBB_LDFLAGS) \
$(LDADD)
stationary_poisson_blocked_counted_LDFLAGS = \
$(AM_LDFLAGS) \
$(DUNEMPILDFLAGS) \
$(TBB_LDFLAGS) \
$(DUNE_LDFLAGS)
noinst_PROGRAMS += stationary-poisson-blocked-timed
stationary_poisson_blocked_timed_SOURCES = stationary-poisson.cc
stationary_poisson_blocked_timed_CPPFLAGS = \
$(AM_CPPFLAGS) \
$(TBB_CPPFLAGS) \
$(DUNEMPICPPFLAGS) \
-DSIMD_BLOCK_SIZE=16 \
-DBLOCK_SOLVER=1 \
-DLOP_TYPE=sumfactorized \
-DENABLE_HP_TIMERS
stationary_poisson_blocked_timed_LDADD = \
$(DUNE_LDFLAGS) $(DUNE_LIBS) \
$(DUNEMPILIBS) \
$(TBB_LIBS) $(TBB_LDFLAGS) \
$(LDADD)
stationary_poisson_blocked_timed_LDFLAGS = \
$(AM_LDFLAGS) \
$(DUNEMPILDFLAGS) \
$(TBB_LDFLAGS) \
$(DUNE_LDFLAGS)
noinst_PROGRAMS += stationary-poisson-flat
stationary_poisson_flat_SOURCES = stationary-poisson.cc
......
......@@ -9,6 +9,10 @@
// enable basis function caching for faster assembly
#define USECACHE 1
#ifdef ENABLE_COUNTER
#include <dune/pdelab/common/opcounter.hh>
#endif
#include<iostream>
#include<vector>
#include<map>
......@@ -20,6 +24,9 @@
#include<dune/common/fvector.hh>
#include<dune/common/static_assert.hh>
#include<dune/common/timer.hh>
#include <dune/pdelab/common/timer.hh>
#include<dune/grid/yaspgrid.hh>
#include<dune/istl/bvector.hh>
#include<dune/istl/operators.hh>
......@@ -269,6 +276,85 @@ private:
};
inline void reset_timers()
{
HP_TIMER_RESET(alpha_volume);
HP_TIMER_RESET(alpha_skeleton);
HP_TIMER_RESET(alpha_boundary);
HP_TIMER_RESET(lambda_volume);
HP_TIMER_RESET(lambda_skeleton);
HP_TIMER_RESET(lambda_boundary);
HP_TIMER_RESET(jacobian_volume);
HP_TIMER_RESET(jacobian_skeleton);
HP_TIMER_RESET(jacobian_boundary);
}
#ifdef ENABLE_COUNTER
#define DUMP_TIMER(name,os,reset) \
os << "===== " #name " =====" << std::endl; \
os << "elapsed: " << HP_TIMER_ELAPSED(name) << std::endl; \
HP_TIMER_OPCOUNTERS(name).reportOperations(os,reset);
#define DUMP_AND_ACCUMULATE_TIMER(name,os,reset,time,ops) \
os << "===== " #name " =====" << std::endl; \
os << "elapsed: " << HP_TIMER_ELAPSED(name) << std::endl; \
time += HP_TIMER_ELAPSED(name); \
ops += HP_TIMER_OPCOUNTERS(name); \
HP_TIMER_OPCOUNTERS(name).reportOperations(os,reset);
#elif defined ENABLE_HP_TIMERS
#define DUMP_TIMER(name,os,reset) \
os << "===== " #name " =====" << std::endl; \
os << "elapsed: " << HP_TIMER_ELAPSED(name) << std::endl;
#define DUMP_AND_ACCUMULATE_TIMER(name,os,reset,time,ops) \
os << "===== " #name " =====" << std::endl; \
os << "elapsed: " << HP_TIMER_ELAPSED(name) << std::endl; \
time += HP_TIMER_ELAPSED(name);
#else
#define DUMP_TIMER(name,os,reset)
#define DUMP_AND_ACCUMULATE_TIMER(name,os,reset,time,ops)
#endif
template<typename Stream>
inline void dump_timers(Stream& os, bool reset)
{
double t = 0.0;
#ifdef ENABLE_COUNTER
auto counter = HP_TIMER_OPCOUNTERS(alpha_volume);
counter.reset();
#endif
DUMP_AND_ACCUMULATE_TIMER(alpha_volume,os,reset,t,counter);
DUMP_AND_ACCUMULATE_TIMER(alpha_skeleton,os,reset,t,counter);
DUMP_AND_ACCUMULATE_TIMER(alpha_boundary,os,reset,t,counter);
DUMP_AND_ACCUMULATE_TIMER(lambda_volume,os,reset,t,counter);
DUMP_AND_ACCUMULATE_TIMER(lambda_skeleton,os,reset,t,counter);
DUMP_AND_ACCUMULATE_TIMER(lambda_boundary,os,reset,t,counter);
DUMP_AND_ACCUMULATE_TIMER(jacobian_volume,os,reset,t,counter);
DUMP_AND_ACCUMULATE_TIMER(jacobian_skeleton,os,reset,t,counter);
DUMP_AND_ACCUMULATE_TIMER(jacobian_boundary,os,reset,t,counter);
os << std::endl
<< "=========== local operator ===========" << std::endl
<< "elapsed: " << t << std::endl;
#ifdef ENABLE_COUNTER
counter.reportOperations(os);
#endif
}
HP_DECLARE_TIMER(total);
//! solve problem with DG method
template<class GV, class PROBLEM, int degree, int blocksize, typename COUT, typename CERR>
void runDG ( const GV& gv,
......@@ -304,7 +390,12 @@ void runDG ( const GV& gv,
Dune::Timer timer;
// coordinate and result type
#ifdef ENABLE_COUNTER
using Real = typename GV::ctype;
#else
typedef double Real;
#endif
const int dim = GV::Grid::dimension;
std::stringstream fullname;
......@@ -465,10 +556,15 @@ void runDG ( const GV& gv,
timer.reset();
cout << "Evaluating residual... " << std::flush;
HP_TIMER_START(total);
go.residual(u,*r);
HP_TIMER_STOP(total);
cout << timer.elapsed() << std::endl;
dump_timers(cout,true);
DUMP_TIMER(total,cout,true);
timer.reset();
cout << "Creating matrix... " << std::flush;
......@@ -494,13 +590,22 @@ void runDG ( const GV& gv,
#endif
cout << nonzeros << " effective non-zero entries" << std::endl;
timer.reset();
cout << "Evaluating jacobian... " << std::flush;
#if not LEGACY_ISTL
raw(*mat).setChunkSize(chunk_size);
#endif
cout << "Evaluating jacobian... " << std::flush;
timer.reset();
HP_TIMER_START(total);
go.jacobian(u,*mat);
HP_TIMER_STOP(total);
cout << timer.elapsed() << std::endl;
dump_timers(cout,true);
DUMP_TIMER(total,cout,true);
cout << timer.elapsed() << std::endl;
......@@ -669,6 +774,7 @@ void runDG ( const GV& gv,
cout << "wallclock time / iteration: " << (solve_time.count() / stat.iterations) << std::endl;
}
#ifndef ENABLE_COUNTER
if(params.get<bool>("io.vtk",false)){
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> UDGF;
UDGF udgf(gfs,u);
......@@ -676,6 +782,7 @@ void runDG ( const GV& gv,
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<UDGF>(udgf,"u_h"));
vtkwriter.write(fullname.str(),Dune::VTK::appendedraw);
}
#endif // ENABLE_COUNTER
}
......@@ -791,24 +898,22 @@ int main(int argc, char** argv)
typedef Grid::LeafGridView GV;
const GV& gv=grid.leafGridView();
typedef Parameter<GV,double> PROBLEM;
typedef Parameter<GV,typename GV::ctype> PROBLEM;
PROBLEM problem(L);
/*
if (degree_dyn==1) {
const int degree=1;
//typedef Dune::PDELab::QkDGLocalFiniteElementMap<Grid::ctype,double,degree,dim> FEMDG;
//FEMDG femdg;
const int blocksize = Dune::QkStuff::QkSize<degree,dim>::value;
runDG<GV,PROBLEM,degree,blocksize>(gv,problem,"CUBE",0,"SIPG","ON",params);
}*/
runDG<GV,PROBLEM,degree,blocksize>(gv,problem,"CUBE",0,"SIPG","ON",params,cout,cerr);
}
if (degree_dyn==2) {
const int degree=2;
//typedef Dune::PDELab::QkDGLocalFiniteElementMap<Grid::ctype,double,degree,dim> FEMDG;
//FEMDG femdg;
const int blocksize = Dune::QkStuff::QkSize<degree,dim>::value;
runDG<GV,PROBLEM,degree,blocksize>(gv,problem,"CUBE",0,"SIPG","ON",params,cout,cerr);
}
}
if (degree_dyn==3) {
const int degree=3;
//typedef Dune::PDELab::QkDGLocalFiniteElementMap<Grid::ctype,double,degree,dim> FEMDG;
......@@ -816,6 +921,27 @@ int main(int argc, char** argv)
const int blocksize = Dune::QkStuff::QkSize<degree,dim>::value;
runDG<GV,PROBLEM,degree,blocksize>(gv,problem,"CUBE",0,"SIPG","ON",params,cout,cerr);
}
if (degree_dyn==4) {
const int degree=4;
//typedef Dune::PDELab::QkDGLocalFiniteElementMap<Grid::ctype,double,degree,dim> FEMDG;
//FEMDG femdg;
const int blocksize = Dune::QkStuff::QkSize<degree,dim>::value;
runDG<GV,PROBLEM,degree,blocksize>(gv,problem,"CUBE",0,"SIPG","ON",params,cout,cerr);
}
if (degree_dyn==5) {
const int degree=5;
//typedef Dune::PDELab::QkDGLocalFiniteElementMap<Grid::ctype,double,degree,dim> FEMDG;
//FEMDG femdg;
const int blocksize = Dune::QkStuff::QkSize<degree,dim>::value;
runDG<GV,PROBLEM,degree,blocksize>(gv,problem,"CUBE",0,"SIPG","ON",params,cout,cerr);
}
if (degree_dyn==6) {
const int degree=6;
//typedef Dune::PDELab::QkDGLocalFiniteElementMap<Grid::ctype,double,degree,dim> FEMDG;
//FEMDG femdg;
const int blocksize = Dune::QkStuff::QkSize<degree,dim>::value;
runDG<GV,PROBLEM,degree,blocksize>(gv,problem,"CUBE",0,"SIPG","ON",params,cout,cerr);
}
if (degree_dyn==7) {
const int degree=7;
//typedef Dune::PDELab::QkDGLocalFiniteElementMap<Grid::ctype,double,degree,dim> FEMDG;
......@@ -823,6 +949,13 @@ int main(int argc, char** argv)
const int blocksize = Dune::QkStuff::QkSize<degree,dim>::value;
runDG<GV,PROBLEM,degree,blocksize>(gv,problem,"CUBE",0,"SIPG","ON",params,cout,cerr);
}
if (degree_dyn==8) {
const int degree=8;
//typedef Dune::PDELab::QkDGLocalFiniteElementMap<Grid::ctype,double,degree,dim> FEMDG;
//FEMDG femdg;
const int blocksize = Dune::QkStuff::QkSize<degree,dim>::value;
runDG<GV,PROBLEM,degree,blocksize>(gv,problem,"CUBE",0,"SIPG","ON",params,cout,cerr);
}
}
catch (Dune::Exception &e)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment