Skip to content

Commit

Permalink
Some play around with DG
Browse files Browse the repository at this point in the history
  • Loading branch information
@blaisb @oguevremont
blaisb authored and oguevremont committed Sep 24, 2024
1 parent fb18bc6 commit 6896f21
Show file tree
Hide file tree
Showing 6 changed files with 192 additions and 26 deletions.
3 changes: 3 additions & 0 deletions include/core/parameters.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1060,6 +1060,9 @@ namespace Parameters
// Interpolation order tracer
unsigned int tracer_order;

// Switch tracer to DG formulation instead of CG
bool tracer_uses_dg;

// Interpolation order vof model
unsigned int VOF_order;

Expand Down
9 changes: 9 additions & 0 deletions include/solvers/copy_data.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,11 +116,20 @@ class StabilizedMethodsCopyData
}
}


struct CopyDataFace
{
FullMatrix<double> cell_matrix;
std::vector<types::global_dof_index> joint_dof_indices;
};

FullMatrix<double> local_matrix;
Vector<double> local_rhs;
std::vector<types::global_dof_index> local_dof_indices;
Vector<double> strong_residual;
std::vector<Vector<double>> strong_jacobian;
std::vector<CopyDataFace> face_data;


// Boolean used to indicate if the cell being assembled is local or not
// This information is used to indicate to the copy_local_to_global function
Expand Down
14 changes: 12 additions & 2 deletions include/solvers/tracer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,8 +83,12 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
else
{
// Usual case, for quad/hex meshes
fe = std::make_shared<FE_Q<dim>>(
simulation_parameters.fem_parameters.tracer_order);
if (simulation_parameters.fem_parameters.tracer_uses_dg)
fe = std::make_shared<FE_DGQ<dim>>(
simulation_parameters.fem_parameters.tracer_order);
else
fe = std::make_shared<FE_Q<dim>>(
simulation_parameters.fem_parameters.tracer_order);
mapping = std::make_shared<MappingQ<dim>>(fe->degree);
cell_quadrature = std::make_shared<QGauss<dim>>(fe->degree + 1);
face_quadrature = std::make_shared<QGauss<dim - 1>>(fe->degree + 1);
Expand Down Expand Up @@ -293,6 +297,12 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
void
assemble_system_matrix() override;

/**
* @brief Assembles the matrix associated with the solver
*/
void
assemble_system_matrix_dg();

/**
* @brief Assemble the rhs associated with the solver
*/
Expand Down
24 changes: 20 additions & 4 deletions include/solvers/tracer_scratch_data.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>

#include <deal.II/fe/fe_interface_values.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping.h>
Expand Down Expand Up @@ -79,6 +80,7 @@ class TracerScratchData
TracerScratchData(const PhysicalPropertiesManager &properties_manager,
const FiniteElement<dim> &fe_tracer,
const Quadrature<dim> &quadrature,
const Quadrature<dim - 1> &face_quadrature,
const Mapping<dim> &mapping,
const FiniteElement<dim> &fe_fd)
: properties_manager(properties_manager)
Expand All @@ -87,6 +89,12 @@ class TracerScratchData
quadrature,
update_values | update_quadrature_points |
update_JxW_values | update_gradients | update_hessians)
, fe_interface_values_tracer(mapping,
fe_tracer,
face_quadrature,
update_values | update_gradients |
update_quadrature_points |
update_JxW_values | update_normal_vectors)
, fe_values_fd(mapping, fe_fd, quadrature, update_values)
{
allocate();
Expand All @@ -112,6 +120,12 @@ class TracerScratchData
sd.fe_values_tracer.get_quadrature(),
update_values | update_quadrature_points |
update_JxW_values | update_gradients | update_hessians)
, fe_interface_values_tracer(sd.fe_interface_values_tracer.get_mapping(),
sd.fe_interface_values_tracer.get_fe(),
sd.fe_interface_values_tracer.get_quadrature(),
update_values | update_gradients |
update_quadrature_points |
update_JxW_values | update_normal_vectors)
, fe_values_fd(sd.fe_values_fd.get_mapping(),
sd.fe_values_fd.get_fe(),
sd.fe_values_fd.get_quadrature(),
Expand Down Expand Up @@ -289,10 +303,12 @@ class TracerScratchData


// FEValues for the Tracer problem
FEValues<dim> fe_values_tracer;
unsigned int n_dofs;
unsigned int n_q_points;
double cell_size;
FEValues<dim> fe_values_tracer;
FEInterfaceValues<dim> fe_interface_values_tracer;

unsigned int n_dofs;
unsigned int n_q_points;
double cell_size;

// Quadrature
std::vector<double> JxW;
Expand Down
6 changes: 6 additions & 0 deletions source/core/parameters.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1463,6 +1463,11 @@ namespace Parameters
"1",
Patterns::Integer(),
"interpolation order chemical potential in the Cahn-Hilliard equations");

prm.declare_entry("tracer uses dg",
"false",
Patterns::Bool(),
"Switch tracer to discontinuous galerkin formulation");
}
prm.leave_subsection();
}
Expand All @@ -1477,6 +1482,7 @@ namespace Parameters
void_fraction_order = prm.get_integer("void fraction order");
temperature_order = prm.get_integer("temperature order");
tracer_order = prm.get_integer("tracer order");
tracer_uses_dg = prm.get_bool("tracer uses dg");
VOF_order = prm.get_integer("VOF order");
phase_cahn_hilliard_order = prm.get_integer("phase cahn hilliard order");
potential_cahn_hilliard_order =
Expand Down
162 changes: 142 additions & 20 deletions source/solvers/tracer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>

#include <deal.II/meshworker/mesh_loop.h>

#include <deal.II/numerics/vector_tools.h>

template <int dim>
Expand Down Expand Up @@ -55,28 +57,38 @@ Tracer<dim>::assemble_system_matrix()
simulation_parameters.source_term.tracer_source->set_time(
simulation_control->get_current_time());

const DoFHandler<dim> *dof_handler_fluid =
multiphysics->get_dof_handler(PhysicsID::fluid_dynamics);
if (simulation_parameters.fem_parameters.tracer_uses_dg)
assemble_system_matrix_dg();

auto scratch_data = TracerScratchData<dim>(
this->simulation_parameters.physical_properties_manager,
*this->fe,
*this->cell_quadrature,
*this->mapping,
dof_handler_fluid->get_fe());
else
{
const DoFHandler<dim> *dof_handler_fluid =
multiphysics->get_dof_handler(PhysicsID::fluid_dynamics);

auto scratch_data = TracerScratchData<dim>(
this->simulation_parameters.physical_properties_manager,
*this->fe,
*this->cell_quadrature,
*this->face_quadrature,
*this->mapping,
dof_handler_fluid->get_fe());

WorkStream::run(this->dof_handler.begin_active(),
this->dof_handler.end(),
*this,
&Tracer::assemble_local_system_matrix,
&Tracer::copy_local_matrix_to_global_matrix,
scratch_data,
StabilizedMethodsCopyData(this->fe->n_dofs_per_cell(),
this->cell_quadrature->size()));
}

WorkStream::run(this->dof_handler.begin_active(),
this->dof_handler.end(),
*this,
&Tracer::assemble_local_system_matrix,
&Tracer::copy_local_matrix_to_global_matrix,
scratch_data,
StabilizedMethodsCopyData(this->fe->n_dofs_per_cell(),
this->cell_quadrature->size()));

system_matrix.compress(VectorOperation::add);
}



template <int dim>
void
Tracer<dim>::assemble_local_system_matrix(
Expand Down Expand Up @@ -180,6 +192,101 @@ Tracer<dim>::copy_local_matrix_to_global_matrix(
}



/// Assemble the system matrix when the system is DG
// This uses the MeshWorker paradigm instead of the WorkStream paradigm
// This is because the DG system matrix is assembled in a different way
template <int dim>
void
Tracer<dim>::assemble_system_matrix_dg()
{
const DoFHandler<dim> *dof_handler_fluid =
multiphysics->get_dof_handler(PhysicsID::fluid_dynamics);

auto scratch_data = TracerScratchData<dim>(
this->simulation_parameters.physical_properties_manager,
*this->fe,
*this->cell_quadrature,
*this->face_quadrature,
*this->mapping,
dof_handler_fluid->get_fe());

StabilizedMethodsCopyData copy_data(this->fe->n_dofs_per_cell(),
this->cell_quadrature->size());

const auto cell_worker =
[&](const typename DoFHandler<dim>::active_cell_iterator &cell,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {
this->assemble_local_system_matrix(cell, scratch_data, copy_data);
};

const auto boundary_worker =
[&](const typename DoFHandler<dim>::active_cell_iterator &cell,
const unsigned int &face_no,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {};

const auto face_worker =
[&](const typename DoFHandler<dim>::active_cell_iterator &cell,
const unsigned int &f,
const unsigned int &sf,
const typename DoFHandler<dim>::active_cell_iterator &ncell,
const unsigned int &nf,
const unsigned int &nsf,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {
FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values_tracer;
fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
const auto &q_points = fe_iv.get_quadrature_points();
copy_data.face_data.emplace_back();

auto &copy_data_face = copy_data.face_data.back();

const unsigned int n_dofs = fe_iv.n_current_interface_dofs();
copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);

const std::vector<double> &JxW = fe_iv.get_JxW_values();
const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
Tensor<1, dim> beta;
beta[0] = 1;
for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
{
const double beta_dot_n = beta * normals[qpoint];
for (unsigned int i = 0; i < n_dofs; ++i)
for (unsigned int j = 0; j < n_dofs; ++j)
copy_data_face.cell_matrix(i, j) +=
fe_iv.jump_in_shape_values(i, qpoint) // [\phi_i]
*
fe_iv.shape_value((beta_dot_n > 0), j, qpoint) // phi_j^{upwind}
* beta_dot_n // (\beta .n)
* JxW[qpoint]; // dx
}
};


const auto copier = [&](const StabilizedMethodsCopyData &c) {
this->copy_local_matrix_to_global_matrix(c);
};



MeshWorker::mesh_loop(this->dof_handler.begin_active(),
this->dof_handler.end(),
cell_worker,
copier,
scratch_data,
copy_data,
MeshWorker::assemble_own_cells |
MeshWorker::assemble_boundary_faces |
MeshWorker::assemble_own_interior_faces_once,
boundary_worker,
face_worker);
}



template <int dim>
void
Tracer<dim>::assemble_system_rhs()
Expand All @@ -201,6 +308,7 @@ Tracer<dim>::assemble_system_rhs()
this->simulation_parameters.physical_properties_manager,
*this->fe,
*this->cell_quadrature,
*this->face_quadrature,
*this->mapping,
dof_handler_fluid->get_fe());

Expand Down Expand Up @@ -967,10 +1075,24 @@ Tracer<dim>::setup_dofs()

// Sparse matrices initialization
DynamicSparsityPattern dsp(locally_relevant_dofs);
DoFTools::make_sparsity_pattern(this->dof_handler,
dsp,
nonzero_constraints,
/*keep_constrained_dofs = */ true);


if (simulation_parameters.fem_parameters.tracer_uses_dg)
{
DoFTools::make_flux_sparsity_pattern(this->dof_handler,
dsp,
nonzero_constraints,
/*keep_constrained_dofs = */ true);
}
else
{
DoFTools::make_sparsity_pattern(this->dof_handler,
dsp,
nonzero_constraints,
/*keep_constrained_dofs = */ true);
}



SparsityTools::distribute_sparsity_pattern(dsp,
locally_owned_dofs,
Expand Down

0 comments on commit 6896f21

Please sign in to comment.