Hi
I am trying a different time integration approach using the generalized
alpha method rather then reposing with a velocity variable. I am also
switching to the BC handled by affine constraints.
I am not sure why the wave wont move through the mesh. I just see the BC at
the edge which quickly fades out in the tutorial.
Any hints or nudges would be appreciated code attached.
Matt
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/279806f2-2cf3-4c0d-84c6-ff44b8e4d7ebn%40googlegroups.com.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/base/utilities.h>
namespace Final_Proj
{
using namespace dealii;
template <int dim>
class WaveEquation
{
public:
WaveEquation();
void run();
private:
void setup_system();
void solve();
void output_results() const;
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> mass_matrix;
SparseMatrix<double> laplace_matrix;
SparseMatrix<double> system_matrix;
SparseMatrix<double> matrix_temp;
Vector<double> solution_u, solution_udot, solution_udotdot;
Vector<double> old_solution_u, old_solution_udot,
old_solution_udotdot;
Vector<double> system_rhs;
double time_step;
double time;
unsigned int timestep_number;
const double beta;
const double gamma;
const double alpha_f;
const double alpha_m;
const double wave_speed;
};
template <int dim>
class ExactSolution : public Function<dim>
{
public:
virtual double value(const Point<dim> &p, unsigned int component =
0) const override
{
(void)component;
double t = this-> get_time();
double exact = -1*std::exp(-1*dim*std::pow(M_PI,2)*t);
for (auto i = 0; i < dim; i++)
exact *= std::sin(M_PI*p[i]); // fix me!!!
return exact;
}
};
template <int dim>
class InitialValuesU : public Function<dim>
{
public:
virtual double value(const Point<dim> & /*p*/, const unsigned int
component = 0) const override
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0;
}
};
template <int dim>
class InitialValuesUdot : public Function<dim>
{
public:
virtual double value(const Point<dim> & /*p*/, const unsigned int
component = 0) const override
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0;
}
};
template <int dim>
class RightHandSide : public Function<dim>
{
public:
virtual double value(const Point<dim> & /*p*/, const unsigned int
component = 0) const override
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0;
}
};
template <int dim>
class BoundaryValuesU : public Function<dim>
{
public:
virtual double value(const Point<dim> & p, const unsigned int
component = 0) const override
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
if ((this->get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3)
&& (p[1] > -1. / 3))
return -std::sin(this->get_time() * 4 * numbers::PI);
else
return 0;
}
};
template <int dim>
WaveEquation<dim>::WaveEquation()
: fe(1)
, dof_handler(triangulation)
, time_step(1. / 64)
, time(time_step)
, timestep_number(1)
, beta(0.25) // unconditionally stable
, gamma(0.5) // unconditionally stable
, alpha_f(0.25)
, alpha_m(0.250)
, wave_speed(1)
{}
template <int dim>
void WaveEquation<dim>::setup_system()
{
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(7);
std::cout << "Number of active cells: " <<
triangulation.n_active_cells()
<< std::endl;
dof_handler.distribute_dofs(fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl
<< std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
mass_matrix.reinit(sparsity_pattern);
laplace_matrix.reinit(sparsity_pattern);
system_matrix.reinit(sparsity_pattern);
MatrixCreator::create_mass_matrix(dof_handler,
QGauss<dim>(fe.degree + 1),
mass_matrix);
MatrixCreator::create_laplace_matrix(dof_handler,
QGauss<dim>(fe.degree + 1),
laplace_matrix);
solution_u.reinit(dof_handler.n_dofs());
solution_udot.reinit(dof_handler.n_dofs());
solution_udotdot.reinit(dof_handler.n_dofs());
old_solution_u.reinit(dof_handler.n_dofs());
old_solution_udot.reinit(dof_handler.n_dofs());
old_solution_udotdot.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
void WaveEquation<dim>::solve()
{
SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
SolverCG<Vector<double>> cg(solver_control);
cg.solve(system_matrix, solution_u, system_rhs, PreconditionIdentity());
constraints.distribute(solution_u);
std::cout << " u-equation: " << solver_control.last_step() << " CG
iterations." << std::endl;
}
template <int dim>
void WaveEquation<dim>::output_results() const
{
// Compute the pointwise maximum error:
Vector<double> max_error_per_cell(triangulation.n_active_cells());
{
ExactSolution<dim> exact_solution;
exact_solution.set_time(time);
MappingQGeneric<dim> mapping(1);
VectorTools::integrate_difference(mapping,
dof_handler,
solution_u,
exact_solution,
max_error_per_cell,
QIterated<dim>(QGauss<1>(2), 2),
VectorTools::NormType::Linfty_norm);
std::cout << "maximum error = " <<
*std::max_element(max_error_per_cell.begin(),max_error_per_cell.end()) <<
std::endl;
}
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution_u, "U");
data_out.add_data_vector(solution_udot, "V");
data_out.add_data_vector(solution_udotdot, "A");
data_out.add_data_vector(max_error_per_cell, "max_error_per_cell");
data_out.build_patches();
const std::string filename = "solution-" +
Utilities::int_to_string(timestep_number, 3) + ".vtu";
DataOutBase::VtkFlags vtk_flags;
vtk_flags.compression_level =
DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
data_out.set_flags(vtk_flags);
std::ofstream output(filename);
data_out.write_vtu(output);
}
}
template <int dim>
void WaveEquation<dim>::run()
{
setup_system();
Vector<double> tmp(solution_u.size());
Vector<double> system_rhs(solution_u.size());
Vector<double> forcing_terms(solution_u.size());
for (; time <= 5; time += time_step, ++timestep_number)
{
BoundaryValuesU<dim> boundary_values_u_function;
boundary_values_u_function.set_time(time);
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
boundary_values_u_function,
constraints
);
constraints.close();
std::cout << "Time step " << timestep_number << " at t=" << time <<
std::endl;
// Building RHS
mass_matrix.vmult(tmp,old_solution_u);
// M*u_n
system_rhs.add((1-alpha_m)/(beta*time_step*time_step),tmp);
// M*(1-alpha_m)/(beta*delt^2)*u_n
mass_matrix.vmult(tmp,old_solution_udot);
//M*udot_n
system_rhs.add((1-alpha_m)/(beta*time_step),tmp);
// M*(1-alpha_m)/(beta*delt^2)*u_n +
M*(1-alpha_m)/(beta*delt)*udot_n
mass_matrix.vmult(tmp,old_solution_udotdot);
//M*udotdot_n
system_rhs.add((1+alpha_m-2*beta)/(2*beta),tmp);
// M*(1-alpha_m)/(beta*delt^2)*u_n +
M*(1-alpha_m)/(beta*delt)*udot_n + M*(1- alpha_m - 2*beta)/(2*beta)*udotdot)
laplace_matrix.vmult(tmp, old_solution_u);
//K*u_n
system_rhs.add(-1*wave_speed*wave_speed*alpha_f,tmp);
// - K*c^2 *alpha_f *u_n + M*(1-alpha_m)/(beta*delt^2)*u_n +
M*(1-alpha_m)/(beta*delt)*udot_n + M*(1- alpha_m - 2*beta)/(2*beta)*udotdot
RightHandSide<dim> rhs_function;
rhs_function.set_time(time);
VectorTools::create_right_hand_side(dof_handler,
QGauss<dim>(fe.degree + 1),
rhs_function,
forcing_terms,
constraints);
//
system_rhs.add((1-alpha_f),forcing_terms);
system_rhs.add(alpha_f,forcing_terms);
// F - K*c^2 *alpha_f *u_n + M*(1-alpha_m)/(beta*delt^2)*u_n +
M*(1-alpha_m)/(beta*delt)*udot_n + M*(1- alpha_m - 2*beta)/(2*beta)*udotdot
// Buidling system matrix
system_matrix.copy_from(mass_matrix);
// M
system_matrix *= (1-alpha_m)/(beta*time_step*time_step);
// M*(1-alpha_m)/(beta*delt^2)
system_matrix.add(wave_speed*wave_speed*(1-
alpha_f),laplace_matrix);
// M*(1-alpha_m)/(beta*delt^2) + c^2 * (1-alpha_f)*K
constraints.condense(system_matrix,system_rhs);
solve();
//solution_udot = gamma/(beta*time_step)*(solution_u-
old_solution_u) -
// (gamma -beta)/beta*old_solution_udot -
//
(gamma-2*beta)/(2*beta)*time_step*old_solution_udotdot;
solution_udot = solution_u;
solution_udot.sadd(gamma/(beta*time_step),-1*gamma/(beta*time_step),old_solution_u);
solution_udot.add(-1*(gamma -beta)/beta,old_solution_udot);
solution_udot.add(-1*(gamma-2*beta)/(2*beta)*time_step,old_solution_udotdot);
//solution_udotdot = 1/(beta*time_step*time_step)*(solution_u-
old_solution_u) -
// 1/(beta*time_step)*old_solution_udot -
// (1-2*beta)/(2*beta)*old_solution_udotdot;
solution_udotdot = solution_u;
solution_udotdot.sadd(1/(beta*time_step*time_step),-1/(beta*time_step*time_step),
old_solution_u);
solution_udotdot.add(-1/(beta*time_step)/beta,old_solution_udot);
solution_udotdot.add((1-2*beta)/(2*beta),old_solution_udotdot);
std::cout << " Total energy: " <<
(mass_matrix.matrix_norm_square(solution_udot) +
laplace_matrix.matrix_norm_square(solution_u)) / 2 << std::endl;
output_results();
old_solution_u = solution_u;
old_solution_udot = solution_udot;
old_solution_udotdot = solution_udotdot;
}
}
}
int main()
{
try
{
using namespace Final_Proj;
WaveEquation<2> wave_equation_solver;
wave_equation_solver.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
