Hi,
Very beginner question. The attached file has a code for 1D advection
equation using runge kutta DG with upwinding flux.
Equation: u_t + u_x = 0;
Actual solution: u = sin(x-t)
Initial contition u(x,0) = sin(x)
Boundary condition: u(0,t) = -sin(t)
I have two issues:
1. The error increases after order 5 or more(clearly visible through plots)
2. The error increases when I take more than "order +1" gauss lobatto points
(even for less than 5th order methods)
Looks like I am missing a very trivial thing. Can someone please help me
out?
Thanks a lot
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
/* to the file deal.II/doc/license.html for the text and */
/* further information on this license. */
#include <base/quadrature_lib.h>
#include <base/logstream.h>
#include <base/function.h>
#include <base/utilities.h>
#include <lac/vector.h>
#include <lac/full_matrix.h>
#include <lac/sparse_matrix.h>
#include <lac/solver_cg.h>
#include <lac/precondition.h>
#include <lac/constraint_matrix.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/grid_tools.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_renumbering.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <fe/fe_dgq.h>
#include <fe/fe_system.h>
#include <fe/fe_values.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
#include <fstream>
#include <sstream>
#include <iostream>
#include <grid/grid_tools.h>
#include <base/tensor_function.h>
#include <lac/solver_richardson.h>
#include <lac/precondition_block.h>
using namespace dealii;
const long double pi = 3.141592653589793238462643;
double a[5] = {0,-567301805773.0/1357537059087.0,
-2404267990393.0/2016746695238.0,
-3550918686646.0/2091501179385.0,
-1275806237668.0/842570457699.0};
double b[5] = {1432997174477.0/9575080441755.0,
5161836677717.0/13612068292357.0,
1720146321549.0/2090206949498.0,
3134564353537.0/4481467310338.0,
2277821191437.0/14882151754819.0};
// Defining Advection Class
template <int dim>
class Advection
{
public:
Advection(const unsigned int degree);
void run();
private:
void make_grid_and_dofs();
void assemble_system_f();
void solve();
void output_results();
const unsigned int degree;
Triangulation<dim> triangulation;
FE_DGQ<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution_f;
Vector<double> mid_solution; //Required for RK4 time marching
Vector<double> old_solution_f;
Vector<double> system_rhs;
FullMatrix<double> hyperbolic_matrix;
Vector<double> hyperbolic_rhs;
Vector<double> hyperbolic_solution;
double time_step;
unsigned int timestep_number;
double time;
};
template <int dim>
Advection<dim>::Advection(const unsigned int degree):
degree(degree),
fe(degree),
dof_handler(triangulation)
{}
class InitialValues:public Function<1>
{
public:
InitialValues():Function<1>() {}
virtual double value (const Point<1> &p,
const unsigned int component = 0) const;
};
double InitialValues::value(const Point <1> &p,
const unsigned int) const
{
return sin(p(0));
}
template <int dim>
void Advection<dim>::make_grid_and_dofs()
{
GridGenerator::hyper_cube(triangulation,0,2);
triangulation.refine_global(4);
dof_handler.distribute_dofs(fe);
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Total number of cells: "
<< triangulation.n_cells()
<< std::endl;
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
sparsity_pattern.reinit(dof_handler.n_dofs(),
dof_handler.n_dofs(),
(degree+1)*dof_handler.max_couplings_between_dofs());
DoFTools::make_flux_sparsity_pattern(dof_handler,sparsity_pattern);
sparsity_pattern.compress();
solution_f.reinit(dof_handler.n_dofs());
mid_solution.reinit(dof_handler.n_dofs());
old_solution_f.reinit(dof_handler.n_dofs());
system_matrix.reinit(sparsity_pattern);
system_rhs.reinit(dof_handler.n_dofs());
const unsigned int dofs_per_cell=fe.dofs_per_cell;
hyperbolic_matrix.reinit(dofs_per_cell,dofs_per_cell);
hyperbolic_rhs.reinit(dofs_per_cell);
hyperbolic_solution.reinit(dofs_per_cell);
}
template <int dim>
void Advection<dim>::assemble_system_f()
{
hyperbolic_matrix = 0;
hyperbolic_rhs =0;
QGaussLobatto<1> quadrature_formula(degree+1);
FEValues<dim> fe_values(fe,quadrature_formula,
update_values | update_gradients |
update_quadrature_points| update_JxW_values);
FEValues<dim> neighbor_fe_values(fe,quadrature_formula,
update_values |update_gradients|
update_quadrature_points
| update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix(dofs_per_cell,dofs_per_cell);
FullMatrix<double> stiff_matrix(dofs_per_cell,dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<unsigned int> local_dof_indices(dofs_per_cell);
std::vector<double> old_solution_values(n_q_points);
std::vector<double> old_solution_values_neighbor(n_q_points);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
const Point<dim> normal = Point<dim>(1.0);
for(;cell!=endc;++cell)
{
cell_matrix=0.;
cell_rhs=0.;
stiff_matrix =0.;
fe_values.reinit(cell);
fe_values.get_function_values(old_solution_f,old_solution_values);
for(unsigned int i=0;i<dofs_per_cell;++i)
for(unsigned int j=0;j<dofs_per_cell;++j)
{
for(unsigned int q_point=0;q_point<n_q_points;++q_point)
cell_matrix(i,j) += (fe_values.shape_value(i,q_point)*
fe_values.shape_value(j,q_point)*
fe_values.JxW(q_point));
for(unsigned int q_point=0;q_point<n_q_points;++q_point)
stiff_matrix(i,j) += (normal*fe_values.shape_grad(i,q_point)*
fe_values.shape_value(j,q_point)*
fe_values.JxW(q_point));
}
for(unsigned int i=0;i<dofs_per_cell;++i)
for(unsigned int j=0;j<dofs_per_cell;++j)
{
double u = old_solution_values[j];
cell_rhs(i) += stiff_matrix(i,j)*u;
}
// Left Neighbor boundary point required for Flux calculation
if(cell->at_boundary(0) == false)
{
typename DoFHandler<dim>::cell_iterator
left_neighbor= cell->neighbor(0);
while (left_neighbor->has_children())
left_neighbor = left_neighbor->child(1);
neighbor_fe_values.reinit (left_neighbor);
neighbor_fe_values.get_function_values(old_solution_f,
old_solution_values_neighbor);
}
else
{
old_solution_values_neighbor[n_q_points-1] = -sin(time);
}
for(unsigned int i=0;i<dofs_per_cell;++i)
{
double u_n = old_solution_values_neighbor[n_q_points-1];
cell_rhs(i) += (u_n)*fe_values.shape_value(i,0);
}
for(unsigned int i=0;i<dofs_per_cell;++i)
{
double u = old_solution_values[n_q_points-1];
cell_rhs(i) -= (u)*fe_values.shape_value(i,n_q_points-1);
}
for(unsigned int i=0;i<dofs_per_cell;++i)
for(unsigned int j=0;j<dofs_per_cell;++j)
hyperbolic_matrix(i,j) = cell_matrix(i,j);
for(unsigned int i=0;i<dofs_per_cell;++i)
hyperbolic_rhs(i) = cell_rhs(i);
solve();
cell->get_dof_indices(local_dof_indices);
for(unsigned int i=0;i<dofs_per_cell;++i)
solution_f(local_dof_indices[i]) = hyperbolic_solution(i);
}
}
template <int dim>
void Advection<dim>::solve()
{
SolverControl solver_control(30000,1e-12,false,false);
SolverRichardson<> solver(solver_control);
solver.solve(hyperbolic_matrix,hyperbolic_solution,
hyperbolic_rhs,PreconditionIdentity());
}
template <int dim>
void Advection<dim>::output_results()
{
if(timestep_number%100 !=0)
return;
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(old_solution_f,"u");
data_out.build_patches ();
const std::string filename = "solution-"+
Utilities::int_to_string(timestep_number,6)+
".gnuplot";
std::ofstream output(filename.c_str());
data_out.write_gnuplot(output);
}
template <int dim>
void Advection<dim>::run()
{
make_grid_and_dofs();
ConstraintMatrix constraints;
constraints.close();
VectorTools::project(dof_handler,constraints,QGaussLobatto<dim>(degree+1),
InitialValues(),old_solution_f);
timestep_number = 0;
time = 0.0;
time_step = 0.001;
mid_solution =0.0;
int n_points = dof_handler.n_dofs();
do
{
output_results();
// Runge Kutta Time Stepping
for(int i=0;i<5;i++)
{
assemble_system_f();
for(int j=0;j<n_points;j++)
{
mid_solution(j) = a[i]*mid_solution(j) + time_step*solution_f(j);
solution_f(j) = old_solution_f(j) + b[i]*mid_solution(j);
}
old_solution_f = solution_f;
}
time +=time_step;
++timestep_number;
}
while(time<=0.41);
}
int main()
{
Advection<1> advection_1d(3);
advection_1d.run();
return 0;
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
