From 6896f21083b2393ad957dd392232e06947840ca2 Mon Sep 17 00:00:00 2001 From: Bruno Blais Date: Mon, 9 Sep 2024 15:25:18 -0400 Subject: [PATCH] Some play around with DG --- include/core/parameters.h | 3 + include/solvers/copy_data.h | 9 ++ include/solvers/tracer.h | 14 ++- include/solvers/tracer_scratch_data.h | 24 +++- source/core/parameters.cc | 6 + source/solvers/tracer.cc | 162 ++++++++++++++++++++++---- 6 files changed, 192 insertions(+), 26 deletions(-) diff --git a/include/core/parameters.h b/include/core/parameters.h index 9a44fa9535..9709a62e3c 100644 --- a/include/core/parameters.h +++ b/include/core/parameters.h @@ -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; diff --git a/include/solvers/copy_data.h b/include/solvers/copy_data.h index bbffcab603..0f2b2209a1 100644 --- a/include/solvers/copy_data.h +++ b/include/solvers/copy_data.h @@ -116,11 +116,20 @@ class StabilizedMethodsCopyData } } + + struct CopyDataFace + { + FullMatrix cell_matrix; + std::vector joint_dof_indices; + }; + FullMatrix local_matrix; Vector local_rhs; std::vector local_dof_indices; Vector strong_residual; std::vector> strong_jacobian; + std::vector 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 diff --git a/include/solvers/tracer.h b/include/solvers/tracer.h index 0d4460d357..6c6de18945 100644 --- a/include/solvers/tracer.h +++ b/include/solvers/tracer.h @@ -83,8 +83,12 @@ class Tracer : public AuxiliaryPhysics else { // Usual case, for quad/hex meshes - fe = std::make_shared>( - simulation_parameters.fem_parameters.tracer_order); + if (simulation_parameters.fem_parameters.tracer_uses_dg) + fe = std::make_shared>( + simulation_parameters.fem_parameters.tracer_order); + else + fe = std::make_shared>( + simulation_parameters.fem_parameters.tracer_order); mapping = std::make_shared>(fe->degree); cell_quadrature = std::make_shared>(fe->degree + 1); face_quadrature = std::make_shared>(fe->degree + 1); @@ -293,6 +297,12 @@ class Tracer : public AuxiliaryPhysics 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 */ diff --git a/include/solvers/tracer_scratch_data.h b/include/solvers/tracer_scratch_data.h index c5dbfa5743..6b13dc0853 100644 --- a/include/solvers/tracer_scratch_data.h +++ b/include/solvers/tracer_scratch_data.h @@ -28,6 +28,7 @@ #include #include +#include #include #include #include @@ -79,6 +80,7 @@ class TracerScratchData TracerScratchData(const PhysicalPropertiesManager &properties_manager, const FiniteElement &fe_tracer, const Quadrature &quadrature, + const Quadrature &face_quadrature, const Mapping &mapping, const FiniteElement &fe_fd) : properties_manager(properties_manager) @@ -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(); @@ -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(), @@ -289,10 +303,12 @@ class TracerScratchData // FEValues for the Tracer problem - FEValues fe_values_tracer; - unsigned int n_dofs; - unsigned int n_q_points; - double cell_size; + FEValues fe_values_tracer; + FEInterfaceValues fe_interface_values_tracer; + + unsigned int n_dofs; + unsigned int n_q_points; + double cell_size; // Quadrature std::vector JxW; diff --git a/source/core/parameters.cc b/source/core/parameters.cc index 75dc1f3b41..419404975b 100644 --- a/source/core/parameters.cc +++ b/source/core/parameters.cc @@ -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(); } @@ -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 = diff --git a/source/solvers/tracer.cc b/source/solvers/tracer.cc index 9b6282268c..272a8bd4c6 100644 --- a/source/solvers/tracer.cc +++ b/source/solvers/tracer.cc @@ -23,6 +23,8 @@ #include #include +#include + #include template @@ -55,28 +57,38 @@ Tracer::assemble_system_matrix() simulation_parameters.source_term.tracer_source->set_time( simulation_control->get_current_time()); - const DoFHandler *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( - this->simulation_parameters.physical_properties_manager, - *this->fe, - *this->cell_quadrature, - *this->mapping, - dof_handler_fluid->get_fe()); + else + { + const DoFHandler *dof_handler_fluid = + multiphysics->get_dof_handler(PhysicsID::fluid_dynamics); + + auto scratch_data = TracerScratchData( + 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 void Tracer::assemble_local_system_matrix( @@ -180,6 +192,101 @@ Tracer::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 +void +Tracer::assemble_system_matrix_dg() +{ + const DoFHandler *dof_handler_fluid = + multiphysics->get_dof_handler(PhysicsID::fluid_dynamics); + + auto scratch_data = TracerScratchData( + 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::active_cell_iterator &cell, + TracerScratchData &scratch_data, + StabilizedMethodsCopyData ©_data) { + this->assemble_local_system_matrix(cell, scratch_data, copy_data); + }; + + const auto boundary_worker = + [&](const typename DoFHandler::active_cell_iterator &cell, + const unsigned int &face_no, + TracerScratchData &scratch_data, + StabilizedMethodsCopyData ©_data) {}; + + const auto face_worker = + [&](const typename DoFHandler::active_cell_iterator &cell, + const unsigned int &f, + const unsigned int &sf, + const typename DoFHandler::active_cell_iterator &ncell, + const unsigned int &nf, + const unsigned int &nsf, + TracerScratchData &scratch_data, + StabilizedMethodsCopyData ©_data) { + FEInterfaceValues &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 ©_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 &JxW = fe_iv.get_JxW_values(); + const std::vector> &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 void Tracer::assemble_system_rhs() @@ -201,6 +308,7 @@ Tracer::assemble_system_rhs() this->simulation_parameters.physical_properties_manager, *this->fe, *this->cell_quadrature, + *this->face_quadrature, *this->mapping, dof_handler_fluid->get_fe()); @@ -967,10 +1075,24 @@ Tracer::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,