Skip to content
Snippets Groups Projects
Commit 21419a18 authored by Dominic Kempf's avatar Dominic Kempf
Browse files

Code uptade to work with Release 2.6

parent 189c0f45
No related branches found
No related tags found
1 merge request!38Piotr's Tutorial 7 on hyperbolic problems
...@@ -24,9 +24,9 @@ void driver (const GV& gv, const FEMDG& femdg, NUMFLUX& numflux, Dune::Parameter ...@@ -24,9 +24,9 @@ void driver (const GV& gv, const FEMDG& femdg, NUMFLUX& numflux, Dune::Parameter
typedef Dune::PDELab::NoConstraints CON; typedef Dune::PDELab::NoConstraints CON;
//necessary to pass proper block size for a component //necessary to pass proper block size for a component
using VBE0 = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::none,Dune::QkStuff::QkSize<FEMDG::order(),dim>::value>; //using VBE0 = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::none,Dune::QkStuff::QkSize<FEMDG::order(),dim>::value>;
//2.6 //2.6
//using VBE0 = Dune::PDELab::ISTL::VectorBackend<>; using VBE0 = Dune::PDELab::ISTL::VectorBackend<>;
using VBE = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed>; using VBE = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed>;
using OrderingTag = Dune::PDELab::EntityBlockedOrderingTag; using OrderingTag = Dune::PDELab::EntityBlockedOrderingTag;
...@@ -96,7 +96,9 @@ void driver (const GV& gv, const FEMDG& femdg, NUMFLUX& numflux, Dune::Parameter ...@@ -96,7 +96,9 @@ void driver (const GV& gv, const FEMDG& femdg, NUMFLUX& numflux, Dune::Parameter
// prepare VTK writer and write first file // prepare VTK writer and write first file
int subsampling=ptree.get("output.subsampling",(int)0); int subsampling=ptree.get("output.subsampling",(int)0);
using VTKWRITER = Dune::SubsamplingVTKWriter<GV>; using VTKWRITER = Dune::SubsamplingVTKWriter<GV>;
VTKWRITER vtkwriter(gv,subsampling); //VTKWRITER vtkwriter(gv,subsampling);
//2.6 todo
VTKWRITER vtkwriter(gv,Dune::refinementLevels(subsampling));
std::string filename=ptree.get("output.filename","output"); std::string filename=ptree.get("output.filename","output");
......
...@@ -7,9 +7,7 @@ ...@@ -7,9 +7,7 @@
#include<dune/common/exceptions.hh> #include<dune/common/exceptions.hh>
#include<dune/common/fvector.hh> #include<dune/common/fvector.hh>
#include<dune/geometry/referenceelements.hh> #include<dune/geometry/referenceelements.hh>
#include<dune/pdelab/common/quadraturerules.hh> #include<dune/pdelab/common/quadraturerules.hh>
#include<dune/pdelab/common/referenceelements.hh>
#include<dune/pdelab/common/geometrywrapper.hh> #include<dune/pdelab/common/geometrywrapper.hh>
#include<dune/pdelab/common/function.hh> #include<dune/pdelab/common/function.hh>
#include<dune/pdelab/localoperator/pattern.hh> #include<dune/pdelab/localoperator/pattern.hh>
......
...@@ -166,120 +166,83 @@ public: ...@@ -166,120 +166,83 @@ public:
{ {
// check for discontinuity // check for discontinuity
if ( Dune::FloatCmp::eq( model().problem.c(inside,x_inside), model().problem.c(outside,x_outside) ) ) if ( Dune::FloatCmp::eq( model().problem.c(inside,x_inside), model().problem.c(outside,x_outside) ) )
{ {
fluxVectorSplitting_.numericalFlux(inside, x_inside, outside, x_outside, n_F, u_s, u_n, f); fluxVectorSplitting_.numericalFlux(inside, x_inside, outside, x_outside, n_F, u_s, u_n, f);
}
/*
// standard case
Dune::FieldMatrix<DF,m,m> D(0.0);
// fetch eigenvalues
model().diagonal(inside,x_inside,D);
Dune::FieldMatrix<DF,m,m> Dplus(0.0);
Dune::FieldMatrix<DF,m,m> Dminus(0.0);
for (size_t i =0 ; i<m;i++)
(D[i][i] > 0) ? Dplus[i][i] = D[i][i] : Dminus[i][i] = D[i][i];
// fetch eigenvectors
Dune::FieldMatrix<DF,m,m> R;
model().eigenvectors(inside,x_inside,n_F,R);
// compute B+ = RD+R^-1 and B- = RD-R^-1
Dune::FieldMatrix<DF,m,m> Bplus(R);
Dune::FieldMatrix<DF,m,m> Bminus(R);
//multiply by D+-
Bplus.rightmultiply(Dplus);
Bminus.rightmultiply(Dminus);
//multiply by R^-1
R.invert();
Bplus.rightmultiply(R);
Bminus.rightmultiply(R);
// Compute numerical flux at the integration point
f = 0.0;
// f = Bplus*u_s + Bminus*u_n
Bplus.umv(u_s,f);
Bminus.umv(u_n,f);
*/
}
else // discontinuous coefficient case else // discontinuous coefficient case
{ {
// fetch eigenvalues // fetch eigenvalues
Dune::FieldMatrix<DF,m,m> D_s(0.0), D_n(0.0); Dune::FieldMatrix<DF,m,m> D_s(0.0), D_n(0.0);
model().diagonal(inside,x_inside,D_s); model().diagonal(inside,x_inside,D_s);
model().diagonal(outside,x_outside,D_n); model().diagonal(outside,x_outside,D_n);
// split positive and negative eigenvalues // split positive and negative eigenvalues
Dune::FieldMatrix<DF,m,m> Dplus_s(0.0); Dune::FieldMatrix<DF,m,m> Dplus_s(0.0);
Dune::FieldMatrix<DF,m,m> Dminus_n(0.0); Dune::FieldMatrix<DF,m,m> Dminus_n(0.0);
for (size_t i=0 ; i<m; i++) for (size_t i=0 ; i<m; i++)
(D_s[i][i] > 0) ? Dplus_s[i][i] = D_s[i][i] : Dminus_n[i][i] = D_n[i][i]; (D_s[i][i] > 0) ? Dplus_s[i][i] = D_s[i][i] : Dminus_n[i][i] = D_n[i][i];
// fetch eigenvectors // fetch eigenvectors
Dune::FieldMatrix<DF,m,m> R_s, R_n; Dune::FieldMatrix<DF,m,m> R_s, R_n;
model().eigenvectors(inside,x_inside,n_F,R_s); model().eigenvectors(inside,x_inside,n_F,R_s);
model().eigenvectors(outside,x_outside,n_F,R_n); model().eigenvectors(outside,x_outside,n_F,R_n);
// compute B+ = RD+R^-1 and B- = RD-R^-1 // compute B+ = RD+R^-1 and B- = RD-R^-1
Dune::FieldMatrix<DF,m,m> Bplus_s(R_s); Dune::FieldMatrix<DF,m,m> Bplus_s(R_s);
Dune::FieldMatrix<DF,m,m> Bminus_n(R_n); Dune::FieldMatrix<DF,m,m> Bminus_n(R_n);
//multiply by D+- //multiply by D+-
Bplus_s.rightmultiply(Dplus_s); Bplus_s.rightmultiply(Dplus_s);
Bminus_n.rightmultiply(Dminus_n); Bminus_n.rightmultiply(Dminus_n);
//multiply by R^-1 //multiply by R^-1
Dune::FieldMatrix<DF,m,m> Rinv_s(R_s), Rinv_n(R_n); Dune::FieldMatrix<DF,m,m> Rinv_s(R_s), Rinv_n(R_n);
Rinv_s.invert(); Rinv_n.invert(); Rinv_s.invert(); Rinv_n.invert();
Bplus_s.rightmultiply(Rinv_s); Bplus_s.rightmultiply(Rinv_s);
Bminus_n.rightmultiply(Rinv_n); Bminus_n.rightmultiply(Rinv_n);
// compute rectangular S matrix // compute rectangular S matrix
Dune::FieldMatrix<DF,m,mstar> S(0.0); Dune::FieldMatrix<DF,m,mstar> S(0.0);
for (int j=0; j<mstar; j++)
if (D_s[j][j]>0.0)
{
for (int i=0; i<m; i++) S[i][j] = R_n[i][j]*D_n[j][j];
}
else
{
for (int i=0; i<m; i++) S[i][j] = -R_s[i][j]*D_s[j][j];
}
// compute square normal matrix
Dune::FieldMatrix<DF,mstar,mstar> SS(0.0);
for (int i=0; i<mstar; i++)
for (int j=0; j<mstar; j++) for (int j=0; j<mstar; j++)
if (D_s[j][j]>0.0) for (int k=0; k<m; k++)
{ SS[i][j] += S[k][i]*S[k][j];
for (int i=0; i<m; i++) S[i][j] = R_n[i][j]*D_n[j][j];
} // compute right hand side of interface problem
else Dune::FieldVector<RF,m> fd(0.0);
{ Bplus_s.umv(u_s,fd);
for (int i=0; i<m; i++) S[i][j] = -R_s[i][j]*D_s[j][j]; Bminus_n.usmv(-1.0,u_n,fd);
} Dune::FieldVector<RF,mstar> rhs(0.0);
for (int i=0; i<mstar; i++)
// compute square normal matrix for (int j=0; j<m; j++)
Dune::FieldMatrix<DF,mstar,mstar> SS(0.0); rhs[i] += S[j][i]*fd[j];
for (int i=0; i<mstar; i++)
for (int j=0; j<mstar; j++) // Solve interface system
for (int k=0; k<m; k++) Dune::FieldVector<RF,mstar> alpha(0.0);
SS[i][j] += S[k][i]*S[k][j]; SS.solve(alpha,rhs);
// compute right hand side of interface problem // extend alpha
Dune::FieldVector<RF,m> fd(0.0); Dune::FieldVector<RF,m> alpha_ext(0.0);
Bplus_s.umv(u_s,fd); for (int i=0; i<mstar; i++)
Bminus_n.usmv(-1.0,u_n,fd); if (D_s[i][i]<0.0) alpha_ext[i] = D_s[i][i]*alpha[i];
Dune::FieldVector<RF,mstar> rhs(0.0);
for (int i=0; i<mstar; i++) // compute flux
for (int j=0; j<m; j++) f = 0.0;
rhs[i] += S[j][i]*fd[j]; Bplus_s.umv(u_s,f);
R_s.umv(alpha_ext,f);
// Solve interface system }
Dune::FieldVector<RF,mstar> alpha(0.0);
SS.solve(alpha,rhs);
// extend alpha
Dune::FieldVector<RF,m> alpha_ext(0.0);
for (int i=0; i<mstar; i++)
if (D_s[i][i]<0.0) alpha_ext[i] = D_s[i][i]*alpha[i];
// compute flux
f = 0.0;
Bplus_s.umv(u_s,f);
R_s.umv(alpha_ext,f);
}
} }
const MODEL& model() const const MODEL& model() const
......
...@@ -33,10 +33,10 @@ public: ...@@ -33,10 +33,10 @@ public:
RF alpha_s(0.0); RF alpha_s(0.0);
RF alpha_n(0.0); RF alpha_n(0.0);
alpha_s = std::abs(u_s[1]/u_s[0]) + sqrt(g*u_s[0]); alpha_s = std::abs(u_s[1]/u_s[0]) + sqrt(g*u_s[0]);
alpha_n = std::abs(u_n[1]/u_n[0]) + sqrt(g*u_n[0]); alpha_n = std::abs(u_n[1]/u_n[0]) + sqrt(g*u_n[0]);
alpha = std::max(alpha_s, alpha_n); alpha = std::max(alpha_s, alpha_n);
} }
//Flux function //Flux function
...@@ -82,14 +82,15 @@ public: ...@@ -82,14 +82,15 @@ public:
RF alpha_n(0.0); RF alpha_n(0.0);
for(size_t k=0;k<dim;++k) for(size_t k=0;k<dim;++k)
{ {
alpha_s += n_F[k]*u_s[k+1]/u_s[0]; alpha_s += n_F[k]*u_s[k+1]/u_s[0];
alpha_n += -n_F[k]*u_n[k+1]/u_n[0]; alpha_n += -n_F[k]*u_n[k+1]/u_n[0];
} }
alpha_s = std::abs(alpha_s) + sqrt(g*u_s[0]);
alpha_n = std::abs(alpha_n) + sqrt(g*u_n[0]); alpha_s = std::abs(alpha_s) + sqrt(g*u_s[0]);
alpha_n = std::abs(alpha_n) + sqrt(g*u_n[0]);
alpha = std::max(alpha_s, alpha_n);
alpha = std::max(alpha_s, alpha_n);
} }
......
...@@ -24,6 +24,9 @@ ...@@ -24,6 +24,9 @@
#if HAVE_UG #if HAVE_UG
#include<dune/grid/uggrid.hh> #include<dune/grid/uggrid.hh>
#endif #endif
#include<dune/geometry/virtualrefinement.hh>
#include<dune/istl/bvector.hh> #include<dune/istl/bvector.hh>
#include<dune/istl/operators.hh> #include<dune/istl/operators.hh>
#include<dune/istl/solvers.hh> #include<dune/istl/solvers.hh>
...@@ -32,9 +35,8 @@ ...@@ -32,9 +35,8 @@
#include<dune/istl/superlu.hh> #include<dune/istl/superlu.hh>
#include<dune/pdelab/finiteelementmap/qkdg.hh> #include<dune/pdelab/finiteelementmap/qkdg.hh>
#include<dune/pdelab/gridfunctionspace/subspace.hh> #include<dune/pdelab/gridfunctionspace/subspace.hh>
#include <dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh> #include<dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh> #include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh> #include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
#include<dune/pdelab/gridfunctionspace/genericdatahandle.hh> #include<dune/pdelab/gridfunctionspace/genericdatahandle.hh>
...@@ -48,9 +50,7 @@ ...@@ -48,9 +50,7 @@
#include<dune/pdelab/gridoperator/onestep.hh> #include<dune/pdelab/gridoperator/onestep.hh>
#include<dune/pdelab/backend/istl.hh> #include<dune/pdelab/backend/istl.hh>
#include<dune/pdelab/instationary/onestep.hh> #include<dune/pdelab/instationary/onestep.hh>
#include<dune/pdelab/function/callableadapter.hh> #include<dune/pdelab/function/callableadapter.hh>
#include<dune/pdelab/gridoperator/onestep.hh> #include<dune/pdelab/gridoperator/onestep.hh>
#include<dune/pdelab/instationary/onestep.hh> #include<dune/pdelab/instationary/onestep.hh>
#include<dune/pdelab/gridoperator/gridoperator.hh> #include<dune/pdelab/gridoperator/gridoperator.hh>
...@@ -107,9 +107,9 @@ int main(int argc, char** argv) ...@@ -107,9 +107,9 @@ int main(int argc, char** argv)
std::bitset<dim> periodic(false); std::bitset<dim> periodic(false);
int overlap=1; int overlap=1;
Dune::array<double,dim> lower_left; for (int i=0; i<dim; i++) lower_left[i]=0.0; std::array<double,dim> lower_left; for (int i=0; i<dim; i++) lower_left[i]=0.0;
auto upper_right = ptree.get<Dune::FieldVector<double,dim> >("grid.L"); auto upper_right = ptree.get<Dune::FieldVector<double,dim> >("grid.L");
auto cells = ptree.get<Dune::array<int,dim> >("grid.N"); auto cells = ptree.get<std::array<int,dim> >("grid.N");
// make grid // make grid
using GM = Dune::YaspGrid<dim>; using GM = Dune::YaspGrid<dim>;
...@@ -130,7 +130,7 @@ int main(int argc, char** argv) ...@@ -130,7 +130,7 @@ int main(int argc, char** argv)
//create numerical flux //create numerical flux
using NUMFLUX = VariableFluxVectorSplitting<MODEL>; using NUMFLUX = VariableFluxVectorSplitting<MODEL>;
// using NUMFLUX = FluxVectorSplitting<MODEL>; //using NUMFLUX = FluxVectorSplitting<MODEL>;
//using NUMFLUX = LLFflux<MODEL>; //using NUMFLUX = LLFflux<MODEL>;
NUMFLUX numflux(model); NUMFLUX numflux(model);
......
...@@ -31,10 +31,11 @@ ...@@ -31,10 +31,11 @@
#include<dune/istl/io.hh> #include<dune/istl/io.hh>
#include<dune/istl/superlu.hh> #include<dune/istl/superlu.hh>
#include<dune/pdelab/finiteelementmap/qkdg.hh> #include<dune/geometry/virtualrefinement.hh>
#include<dune/pdelab/finiteelementmap/qkdg.hh>
#include<dune/pdelab/gridfunctionspace/subspace.hh> #include<dune/pdelab/gridfunctionspace/subspace.hh>
#include <dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh> #include<dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh> #include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh> #include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
#include<dune/pdelab/gridfunctionspace/genericdatahandle.hh> #include<dune/pdelab/gridfunctionspace/genericdatahandle.hh>
...@@ -48,9 +49,7 @@ ...@@ -48,9 +49,7 @@
#include<dune/pdelab/gridoperator/onestep.hh> #include<dune/pdelab/gridoperator/onestep.hh>
#include<dune/pdelab/backend/istl.hh> #include<dune/pdelab/backend/istl.hh>
#include<dune/pdelab/instationary/onestep.hh> #include<dune/pdelab/instationary/onestep.hh>
#include<dune/pdelab/function/callableadapter.hh> #include<dune/pdelab/function/callableadapter.hh>
#include<dune/pdelab/gridoperator/onestep.hh> #include<dune/pdelab/gridoperator/onestep.hh>
#include<dune/pdelab/instationary/onestep.hh> #include<dune/pdelab/instationary/onestep.hh>
#include<dune/pdelab/gridoperator/gridoperator.hh> #include<dune/pdelab/gridoperator/gridoperator.hh>
...@@ -107,9 +106,9 @@ int main(int argc, char** argv) ...@@ -107,9 +106,9 @@ int main(int argc, char** argv)
std::bitset<dim> periodic(false); std::bitset<dim> periodic(false);
int overlap=1; int overlap=1;
Dune::array<double,dim> lower_left; for (int i=0; i<dim; i++) lower_left[i]=0.0; std::array<double,dim> lower_left; for (int i=0; i<dim; i++) lower_left[i]=0.0;
auto upper_right = ptree.get<Dune::FieldVector<double,dim> >("grid.L"); auto upper_right = ptree.get<Dune::FieldVector<double,dim> >("grid.L");
auto cells = ptree.get<Dune::array<int,dim> >("grid.N"); auto cells = ptree.get<std::array<int,dim> >("grid.N");
// make grid // make grid
using GM = Dune::YaspGrid<dim>; using GM = Dune::YaspGrid<dim>;
......
...@@ -19,6 +19,7 @@ ...@@ -19,6 +19,7 @@
#include<dune/grid/io/file/vtk.hh> #include<dune/grid/io/file/vtk.hh>
#include<dune/grid/io/file/vtk/subsamplingvtkwriter.hh> #include<dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include<dune/grid/io/file/gmshreader.hh> #include<dune/grid/io/file/gmshreader.hh>
#include<dune/grid/yaspgrid.hh> #include<dune/grid/yaspgrid.hh>
#include<dune/grid/onedgrid.hh> #include<dune/grid/onedgrid.hh>
...@@ -35,7 +36,7 @@ ...@@ -35,7 +36,7 @@
#include<dune/pdelab/finiteelementmap/qkdg.hh> #include<dune/pdelab/finiteelementmap/qkdg.hh>
#include<dune/pdelab/gridfunctionspace/subspace.hh> #include<dune/pdelab/gridfunctionspace/subspace.hh>
#include <dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh> #include<dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh> #include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh> #include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
#include<dune/pdelab/gridfunctionspace/genericdatahandle.hh> #include<dune/pdelab/gridfunctionspace/genericdatahandle.hh>
...@@ -49,9 +50,7 @@ ...@@ -49,9 +50,7 @@
#include<dune/pdelab/gridoperator/onestep.hh> #include<dune/pdelab/gridoperator/onestep.hh>
#include<dune/pdelab/backend/istl.hh> #include<dune/pdelab/backend/istl.hh>
#include<dune/pdelab/instationary/onestep.hh> #include<dune/pdelab/instationary/onestep.hh>
#include<dune/pdelab/function/callableadapter.hh> #include<dune/pdelab/function/callableadapter.hh>
#include<dune/pdelab/gridoperator/onestep.hh> #include<dune/pdelab/gridoperator/onestep.hh>
#include<dune/pdelab/instationary/onestep.hh> #include<dune/pdelab/instationary/onestep.hh>
#include<dune/pdelab/gridoperator/gridoperator.hh> #include<dune/pdelab/gridoperator/gridoperator.hh>
...@@ -108,9 +107,9 @@ int main(int argc, char** argv) ...@@ -108,9 +107,9 @@ int main(int argc, char** argv)
// read grid parameters from input file // read grid parameters from input file
using DF = Dune::OneDGrid::ctype; using DF = Dune::OneDGrid::ctype;
auto a = 0; auto a = 0;
Dune::array<double,dim> lower_left; for (int i=0; i<dim; i++) lower_left[i]=0.0; std::array<double,dim> lower_left; for (int i=0; i<dim; i++) lower_left[i]=0.0;
auto upper_right = ptree.get<Dune::FieldVector<double,dim> >("grid.L"); auto upper_right = ptree.get<Dune::FieldVector<double,dim> >("grid.L");
auto cells = ptree.get<Dune::array<int,dim> >("grid.N"); auto cells = ptree.get<std::array<int,dim> >("grid.N");
auto b = upper_right[0]; auto b = upper_right[0];
auto N = cells[0]; auto N = cells[0];
...@@ -175,9 +174,9 @@ int main(int argc, char** argv) ...@@ -175,9 +174,9 @@ int main(int argc, char** argv)
std::bitset<dim> periodic(false); std::bitset<dim> periodic(false);
int overlap=1; int overlap=1;
Dune::array<double,dim> lower_left; for (int i=0; i<dim; i++) lower_left[i]=0.0; std::array<double,dim> lower_left; for (int i=0; i<dim; i++) lower_left[i]=0.0;
auto upper_right = ptree.get<Dune::FieldVector<double,dim> >("grid.L"); auto upper_right = ptree.get<Dune::FieldVector<double,dim> >("grid.L");
auto cells = ptree.get<Dune::array<int,dim> >("grid.N"); auto cells = ptree.get<std::array<int,dim> >("grid.N");
// make grid // make grid
using GM = Dune::YaspGrid<dim>; using GM = Dune::YaspGrid<dim>;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment