Skip to content
Snippets Groups Projects
Commit 5942092e authored by Peter Bastian's avatar Peter Bastian
Browse files

added support and example with discontinuous coefficients

parent ce947833
No related branches found
No related tags found
1 merge request!38Piotr's Tutorial 7 on hyperbolic problems
Hide whitespace changes
Inline Side-by-side
tutorial07/src/hyperbolicdg.hh
+ 3
− 12
View file @ 5942092e
......@@ -90,7 +90,6 @@ namespace Dune {
// evaluate speed of sound (assumed constant per element)
auto ref_el = referenceElement(geo);
auto localcenter = ref_el.position(0,0);
auto c = numflux.model.problem.c(cell,localcenter);
// Transformation
typename EG::Geometry::JacobianInverseTransposed jac;
......@@ -124,7 +123,7 @@ namespace Dune {
Dune::FieldMatrix<RF,m,dim> F;
numflux.model.flux(u,F);
numflux.model.flux(eg,ip.position(),u,F);
// integrate
auto factor = ip.weight() * geo.integrationElement(ip.position());
......@@ -191,12 +190,9 @@ namespace Dune {
auto ref_el_outside = referenceElement(geo_outside);
auto local_inside = ref_el_inside.position(0,0);
auto local_outside = ref_el_outside.position(0,0);
//auto c_s = numflux.model.problem.c(cell_inside,local_inside);
//auto c_n = numflux.model.problem.c(cell_outside,local_outside);
// for now assume that c is constant
// the case that non-homogenious coefficient we leave for the future
//auto c = c_s;
// Initialize vectors outside for loop
Dune::FieldVector<RF,m> u_s(0.0);
......@@ -229,7 +225,7 @@ namespace Dune {
u_n[i] += x_n(lfsv_n.child(i),k)*phi_n[k];
// Compute numerical flux at the integration point
numflux.numericalFlux(ig,x_s,u_s,u_n,f);
numflux.numericalFlux(cell_inside,iplocal_s,cell_outside,iplocal_n,ig.centerUnitOuterNormal(),u_s,u_n,f);
// Integrate
auto factor = ip.weight() * geo.integrationElement(ip.position());
......@@ -282,13 +278,8 @@ namespace Dune {
// Evaluate speed of sound (assumed constant per element)
auto ref_el_inside = referenceElement(geo_inside);
auto local_inside = ref_el_inside.position(0,0);
auto c_s = numflux.model.problem.c(cell_inside,local_inside);
// for now assume that c is constant
// the case that non-homogenious coefficient we leave for the future
auto c = c_s;
// Initialize vectors outside for loop
Dune::FieldVector<RF,m> u_s(0.0);
Dune::FieldVector<RF,m> f(0.0);
......@@ -314,7 +305,7 @@ namespace Dune {
Dune::FieldVector<RF,m> u_n(numflux.model.problem.g(ig.intersection(),ip.position(),u_s));
// Compute numerical flux at integration point
numflux.numericalFlux(ig,x_s,u_s,u_n,f);
numflux.numericalFlux(cell_inside,iplocal_s,cell_inside,iplocal_s,ig.centerUnitOuterNormal(),u_s,u_n,f);
// Integrate
auto factor = ip.weight() * geo.integrationElement(ip.position());
......
tutorial07/src/linearacoustics.hh
+ 27
− 57
View file @ 5942092e
......@@ -16,10 +16,12 @@ class Model ;
template<typename PROBLEM>
class Model<2,PROBLEM>
{
public:
static constexpr int dim = 2;
static constexpr int m = 3;
static constexpr int dim = 2; // space dimension
static constexpr int m = dim+1; // system size
static constexpr int mplus = 1; // number of positive eigenvectors
static constexpr int mminus = 1; // number of negative eigenvectors
static constexpr int mstar = mplus+mminus; // number of nonzero eigenvalues
using RangeField = typename PROBLEM::RangeField;
......@@ -28,74 +30,42 @@ public:
{
}
//TODO discontinuos coefficients
/*
template<typename E, typename X, typename T2, typename T3>
static void eigenvectors (
const E& e, const X& x,
const Dune::FieldVector<T2,dim>& n,
Dune::FieldMatrix<T3,m,m>& RT
)
void eigenvectors (const E& e, const X& x, const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,m,m>& RT) const
{
auto c = problem.c(e,x);
RT[0][0] = 0; RT[1][0] = -n[1]; RT[2][0] = n[0];
RT[0][1] = 1; RT[1][1] = c*n[0]; RT[2][1] = c*n[1];
RT[0][2] = -1; RT[1][2] = c*n[0]; RT[2][2] = c*n[1];
RT[0][0] = 1.0/c; RT[0][1] = -1.0/c; RT[0][2] = 0.0;
RT[1][0] = n[0]; RT[1][1] = n[0]; RT[1][2] = -n[1];
RT[2][0] = n[1]; RT[2][1] = n[1]; RT[2][2] = n[0];
}
*/
template<typename T2, typename T3>
static void eigenvectors ( const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,m,m>& RT)
{
//auto c = problem.c(e,x);
int c = 1;
// where do we need this matrix for?
// template<typename RF>
// static void coefficients (Dune::FieldMatrix<RF,m,m>& A)
// {
// int c2 = 1;
RT[0][0] = 0; RT[1][0] = -n[1]; RT[2][0] = n[0];
RT[0][1] = 1; RT[1][1] = c*n[0]; RT[2][1] = c*n[1];
RT[0][2] = -1; RT[1][2] = c*n[0]; RT[2][2] = c*n[1];
}
//one can also provide eigenvectors inverse
// A[0][0] = 0.0; A[0][1] = 1.0; A[0][2] = 1.0;
// A[1][0] = c2; A[1][1] = 0.0; A[1][2] = 0.0;
// A[2][0] = c2; A[2][1] = 0.0; A[2][2] = 0.0;
// }
template<typename RF>
static void coefficients (Dune::FieldMatrix<RF,m,m>& A)
template<typename E, typename X, typename RF>
void diagonal (const E& e, const X& x, Dune::FieldMatrix<RF,m,m>& D) const
{
int c2 = 1;
A[0][0] = 0.0; A[0][1] = 1.0; A[0][2] = 1.0;
A[1][0] = c2; A[1][1] = 0.0; A[1][2] = 0.0;
A[2][0] = c2; A[2][1] = 0.0; A[2][2] = 0.0;
}
template<typename RF>
static void diagonal (Dune::FieldMatrix<RF,m,m>& D)
{
int c = 1;
auto c = problem.c(e,x);
D[0][0] = 0.0; D[0][1] = 0.0; D[0][2] = 0.0;
D[1][0] = 0.0; D[1][1] = c ; D[1][2] = 0.0;
D[2][0] = 0.0; D[2][1] = 0.0; D[2][2] = -c;
D[0][0] = c; D[0][1] = 0.0; D[0][2] = 0.0;
D[1][0] = 0.0; D[1][1] = -c ; D[1][2] = 0.0;
D[2][0] = 0.0; D[2][1] = 0.0; D[2][2] = 0.0;
}
//Flux function
//template<typename RF,typename EG>
template<typename RF>
static void flux (const Dune::FieldVector<RF,m>& u, Dune::FieldMatrix<RF,m,dim>& F)
template<typename E, typename X, typename RF>
void flux (const E& e, const X& x, const Dune::FieldVector<RF,m>& u, Dune::FieldMatrix<RF,m,dim>& F) const
{
//fetch parameters
/*
// Reference to cell
const auto& cell = eg.entity();
// Get geometry
auto geo = eg.geometry();
// evaluate speed of sound (assumed constant per element)
auto ref_el = referenceElement(geo);
auto localcenter = ref_el.position(0,0);
auto c = param.c(cell,localcenter);
*/
int c = 1;
auto c = problem.c(e,x);
F[0][0] = u[1] ; F[0][1] = u[2];
F[1][0] = c*c*u[0]; F[1][1] = 0.0;
......
tutorial07/src/linearacousticsproblem.hh
+ 7
− 5
View file @ 5942092e
......@@ -9,7 +9,7 @@ public:
//problem specification depends on dimension
static constexpr int dim = 2;
static constexpr int m = 3;
static constexpr int m = dim+1;
using Range = Dune::FieldVector<NUMBER,m>;
......@@ -22,10 +22,12 @@ public:
template<typename E, typename X>
NUMBER c (const E& e, const X& x) const
{
X xglobal = e.geometry().global(x);
if ( xglobal[1] < 1-(0.6/0.9)*(xglobal[0]-0.1) ) return 1.0;
return 0.5;
auto xglobal = e.geometry().center();
// if (xglobal[0]>1.0 && xglobal[0]<2.0 && xglobal[1]>0.375 && xglobal[1]<0.625)
if (xglobal[1]>0.625)
return 0.33333;
else
return 1.0;
}
//! Neumann boundary condition
......
tutorial07/src/maxwell.hh
+ 6
− 6
View file @ 5942092e
......@@ -29,8 +29,8 @@ public:
{
}
template<typename T1, typename T2>
static void eigenvectors (const Dune::FieldVector<T1,dim>& n, Dune::FieldMatrix<T2,m,m>& R)
template<typename E, typename X, typename T1, typename T2>
void eigenvectors (const E& e, const X& x, const Dune::FieldVector<T1,dim>& n, Dune::FieldMatrix<T2,m,m>& R) const
{
int eps = 1.0;
int mu = 1.0;
......@@ -102,8 +102,8 @@ public:
A[5][0] =-1./mu; A[5][1] = 1./mu; A[5][2] = 0.0; A[5][3] = 0.0; A[5][4] = 1.0; A[5][5] = 1.0;
}
template<typename RF>
static void diagonal (Dune::FieldMatrix<RF,m,m>& D)
template<typename E, typename X, typename RF>
void diagonal (const E& e, const X& x, Dune::FieldMatrix<RF,m,m>& D) const
{
RF c(1.0);
......@@ -116,8 +116,8 @@ public:
}
//Flux function
template<typename RF>
static void flux (Dune::FieldVector<RF,m>& u, Dune::FieldMatrix<RF,m,dim>& F)
template<typename E, typename X, typename RF>
void flux (const E& e, const X& x, const Dune::FieldVector<RF,m>& u, Dune::FieldMatrix<RF,m,dim>& F) const
{
RF mu(1.0);
RF ep(1.0);
......
tutorial07/src/numericalflux.hh
+ 155
− 17
View file @ 5942092e
......@@ -36,21 +36,18 @@ public:
{
}
template<typename IG, typename X>
static void numericalFlux(const IG& ig, const X& x,
template<typename E, typename X>
void numericalFlux(const E& inside, const X& x_inside, const E& outside, const X& x_outside, const Dune::FieldVector<DF,dim> n_F,
const Dune::FieldVector<RF,m>& u_s,
const Dune::FieldVector<RF,m>& u_n,Dune::FieldVector<RF,m>& f)
const Dune::FieldVector<RF,m>& u_n,Dune::FieldVector<RF,m>& f) const
{
// Normal: assume faces are planar
const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
// fetch flux
Dune::FieldMatrix<RF,m,dim> Fs;
Dune::FieldMatrix<RF,m,dim> Fn;
//evaluate flux
MODEL::flux(u_s,Fs);
MODEL::flux(u_n,Fn);
model.flux(inside,x_inside,u_s,Fs);
model.flux(outside,x_outside,u_n,Fn);
//Fs*n_F + Fn*n_F
Fs.umv(n_F,f);
......@@ -59,7 +56,7 @@ public:
Dune::FieldMatrix<DF,m,m> D(0.0);
// fetch eigenvalues
MODEL::diagonal(D);
model.diagonal(inside,x_inside,D);
//max eigenvalue
RF alpha(0.0);
......@@ -94,19 +91,16 @@ public:
{
}
template<typename IG, typename X>
static void numericalFlux(const IG& ig, const X& x,
template<typename E, typename X>
void numericalFlux(const E& inside, const X& x_inside, const E& outside, const X& x_outside, const Dune::FieldVector<DF,dim> n_F,
const Dune::FieldVector<RF,m>& u_s,
const Dune::FieldVector<RF,m>& u_n,Dune::FieldVector<RF,m>& f)
const Dune::FieldVector<RF,m>& u_n,Dune::FieldVector<RF,m>& f) const
//const std::vector<Dune::FieldVector<RF,dim> >& gradu_s,
//const std::vector<Dune::FieldVector<RF,dim> >& gradu_n)
{
// Normal: assume faces are planar
const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
Dune::FieldMatrix<DF,m,m> D(0.0);
// fetch eigenvalues
MODEL::diagonal(D);
model.diagonal(inside,x_inside,D);
Dune::FieldMatrix<DF,m,m> Dplus(0.0);
Dune::FieldMatrix<DF,m,m> Dminus(0.0);
......@@ -116,7 +110,7 @@ public:
// fetch eigenvectors
Dune::FieldMatrix<DF,m,m> Rot;
MODEL::eigenvectors(n_F,Rot);
model.eigenvectors(inside,x_inside,n_F,Rot);
// compute B+ = RD+R^-1 and B- = RD-R^-1
Dune::FieldMatrix<DF,m,m> Bplus(Rot);
......@@ -136,10 +130,154 @@ public:
// f = Bplus*u_s + Bminus*u_n
Bplus.umv(u_s,f);
Bminus.umv(u_n,f);
}
const MODEL& model;
};// FVS
//Flux Vector splitting for discontinuous coefficients
template<typename MODEL>
class VariableFluxVectorSplitting
{
public:
static constexpr int dim = MODEL::dim;
static constexpr int m = MODEL::Model::m;
static constexpr int mstar = MODEL::mstar;
using RF = typename MODEL::RangeField; // type for computations
using DF = RF;
VariableFluxVectorSplitting (const MODEL& model_)
: model(model_)
{
}
template<typename E, typename X>
void numericalFlux(const E& inside, const X& x_inside,
const E& outside, const X& x_outside,
const Dune::FieldVector<DF,dim> n_F,
const Dune::FieldVector<RF,m>& u_s,
const Dune::FieldVector<RF,m>& u_n,Dune::FieldVector<RF,m>& f) const
{
// check for discontinuity
if (model.problem.c(inside,x_inside)==model.problem.c(outside,x_outside))
{
// 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
{
// fetch eigenvalues
Dune::FieldMatrix<DF,m,m> D_s(0.0), D_n(0.0);
model.diagonal(inside,x_inside,D_s);
model.diagonal(outside,x_outside,D_n);
// split positive and negative eigenvalues
Dune::FieldMatrix<DF,m,m> Dplus_s(0.0);
Dune::FieldMatrix<DF,m,m> Dminus_n(0.0);
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];
// fetch eigenvectors
Dune::FieldMatrix<DF,m,m> R_s, R_n;
model.eigenvectors(inside,x_inside,n_F,R_s);
model.eigenvectors(outside,x_outside,n_F,R_n);
// compute B+ = RD+R^-1 and B- = RD-R^-1
Dune::FieldMatrix<DF,m,m> Bplus_s(R_s);
Dune::FieldMatrix<DF,m,m> Bminus_n(R_n);
//multiply by D+-
Bplus_s.rightmultiply(Dplus_s);
Bminus_n.rightmultiply(Dminus_n);
//multiply by R^-1
Dune::FieldMatrix<DF,m,m> Rinv_s(R_s), Rinv_n(R_n);
Rinv_s.invert(); Rinv_n.invert();
Bplus_s.rightmultiply(Rinv_s);
Bminus_n.rightmultiply(Rinv_n);
// compute rectangular S matrix
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 k=0; k<m; k++)
SS[i][j] += S[k][i]*S[k][j];
// compute right hand side of interface problem
Dune::FieldVector<RF,m> fd(0.0);
Bplus_s.umv(u_s,fd);
Bminus_n.usmv(-1.0,u_n,fd);
Dune::FieldVector<RF,mstar> rhs(0.0);
for (int i=0; i<mstar; i++)
for (int j=0; j<m; j++)
rhs[i] += S[j][i]*fd[j];
// 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;
};// FVS
#endif //NUMERICALFLUX
tutorial07/src/tutorial07-acoustics.cc
+ 3
− 2
View file @ 5942092e
......@@ -129,8 +129,9 @@ int main(int argc, char** argv)
MODEL model(problem);
//create numerical flux
//using NUMFLUX = FluxVectorSplitting<MODEL>;
using NUMFLUX = LLFflux<MODEL>;
using NUMFLUX = VariableFluxVectorSplitting<MODEL>;
// using NUMFLUX = FluxVectorSplitting<MODEL>;
//using NUMFLUX = LLFflux<MODEL>;
NUMFLUX numflux(model);
......
tutorial07/src/tutorial07-acoustics.ini
+ 5
− 5
View file @ 5942092e
......@@ -7,14 +7,14 @@ N = 5 2
refinement = 4
[fem]
degree=1
torder=1
dt=0.001
degree=2
torder=3
dt=0.0024
[problem]
T=2.0
[output]
filename = acoustics_dg
subsampling = 0
every = 1
subsampling = 2
every = 5
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