Add transient_ex3 example demonstrating explicit DG/FV formulation of 2D advection equation

configure.ac

@@ -471,6 +471,7 @@ AS_IF([test "x$enableexamples" = "xyes"],
                          examples/reduced_basis/reduced_basis_ex7/Makefile
                          examples/transient/transient_ex1/Makefile
                          examples/transient/transient_ex2/Makefile
+                         examples/transient/transient_ex3/Makefile
                          examples/vector_fe/vector_fe_ex1/Makefile
                          examples/vector_fe/vector_fe_ex2/Makefile
                          examples/vector_fe/vector_fe_ex3/Makefile

examples/Makefile.am

@@ -66,6 +66,7 @@ SUBDIRS = \
 	reduced_basis/reduced_basis_ex7 \
 	transient/transient_ex1 \
 	transient/transient_ex2 \
+	transient/transient_ex3 \
 	vector_fe/vector_fe_ex1 \
 	vector_fe/vector_fe_ex2 \
 	vector_fe/vector_fe_ex3 \

examples/transient/transient_ex3/Makefile.am

@@ -0,0 +1,20 @@
+example_name  = transient_ex3
+check_SCRIPTS = run.sh
+install_dir   = $(examples_install_path)/transient/ex3
+data          = advection_2D.in advection_system.C advection_system.h claw_system.C claw_system.h transient_ex3.C run.sh
+sources       = $(data) run.sh
+
+CLEANFILES = claw_solution.*.e
+
+if LIBMESH_VPATH_BUILD
+  BUILT_SOURCES = .linkstamp
+.linkstamp:
+	-rm -f advection_2D.in && $(LN_S) -f $(srcdir)/advection_2D.in .
+	$(AM_V_GEN)touch .linkstamp
+
+  CLEANFILES += advection_2D.in .linkstamp
+endif
+
+##############################################
+# include common example environment
+include $(top_srcdir)/examples/Make.common

examples/transient/transient_ex3/advection_2D.in

@@ -0,0 +1,30 @@
+# The number of elements in each dimension in the mesh
+n_elem = 30
+
+# The advection velocity vector
+u1 = 1
+u2 = 1
+
+# The Lax-Friedrichs flux constant
+# We guarantee stability for the advection problem if C >= |u|
+LxF_constant = 1.414
+
+# The time-step size
+delta_t = 0.001
+
+# The number of time-steps
+n_time_steps = 500
+
+# The time-step interval at which we write out solutions
+write_interval = 50
+
+# The temporal discretization type
+temporal_discretization_type = RK4
+
+# Equivalent to a Finite Volume (FV) formulation
+fe_order = CONSTANT
+fe_family = MONOMIAL
+
+# Discontinuous Galerkin (DG) formulation
+# fe_order = SECOND
+# fe_family = L2_LAGRANGE

examples/transient/transient_ex3/advection_system.C

@@ -0,0 +1,137 @@
+// local includes
+#include "advection_system.h"
+
+// LibMesh includes
+#include "libmesh/equation_systems.h" // EquationSystems::comm()
+#include "libmesh/getpot.h" // GetPot input file parsing
+#include "libmesh/libmesh_logging.h" // LOG_SCOPE
+#include "libmesh/numeric_vector.h"
+#include "libmesh/sparse_matrix.h"
+#include "libmesh/string_to_enum.h"
+
+namespace libMesh
+{
+
+AdvectionSystem::AdvectionSystem (EquationSystems& es,
+		                  const std::string& name,
+		                  const unsigned int number)
+  : Parent(es, name, number),
+    _q1_var(0),
+    _fe_order(CONSTANT),
+    _fe_family(MONOMIAL)
+{
+  // Allocate NumericVectors in _Fh. I could not figure out
+  // how to do this in the initialization list, I don't think it's
+  // possible.
+  // TODO: the number of problem dimensions is hard-coded here, we should
+  // make this depend on the input file parameters instead.
+  for (unsigned int i=0; i<2; ++i)
+    _Fh.push_back(NumericVector<Number>::build(es.comm()));
+}
+
+
+AdvectionSystem::~AdvectionSystem() = default;
+
+void AdvectionSystem::assemble_claw_rhs (NumericVector<Number> & q)
+{
+  LOG_SCOPE("assemble_claw_rhs()", "AdvectionSystem");
+
+  // The input to this function is the solution vector, "q", from
+  // either the initial condition or the previous timestep.
+  this->update_Fh(q);
+
+  // Allocate storage for temporary vector used in the computations below
+  auto temp = NumericVector<Number>::build(this->comm());
+  temp->init(this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
+
+  rhs->zero();
+
+  // rhs += (Ai - AVG_i)*Fi
+  for (auto i : index_range(_Fh))
+  {
+    this->get_advection_matrix(i).vector_mult(*temp, *_Fh[i]); // temp = Ai*Fi
+    rhs->add(+1., *temp);
+    this->get_avg_matrix(i).vector_mult(*temp, *_Fh[i]); // temp = AVG_i*Fi
+    rhs->add(-1., *temp);
+  }
+
+  // Jump term
+  // rhs -= LxF * (J*q)
+  this->get_jump_matrix().vector_mult(*temp, q);
+  rhs->add(-get_LxF_constant(), *temp);
+
+  // Apply boundary conditions. This is also expressed as a matrix-vector product
+  // rhs -= BC_i * Fh[i]
+  for (auto i : index_range(_Fh))
+  {
+    this->get_boundary_condition_matrix(i).vector_mult(*temp, *_Fh[i]); // temp = BC_i * Fi
+    rhs->add(-1., *temp);
+  }
+}
+
+void AdvectionSystem::update_Fh (NumericVector<Number> & q)
+{
+  LOG_SCOPE("update_Fh()", "AdvectionSystem");
+
+  // Fi = u(i)*q
+  for (auto i : index_range(_Fh))
+  {
+    _Fh[i]->zero();
+    _Fh[i]->add(_u(i), q);
+    _Fh[i]->close();
+  }
+}
+
+void AdvectionSystem::init_data ()
+{
+  // Print info about the type of discretization being used.  Note: we
+  // can use the same exact assembly code while changing only the
+  // approximation order to test both "DG" and "FV" discretizations.
+  // In the FV case, the "interior" integral contributions are zero
+  // since they depend on "dphi".
+  libMesh::out << "Adding q1 variable using ("
+               << Utility::enum_to_string(_fe_order)
+               << ", "
+               << Utility::enum_to_string(_fe_family)
+               << ") approximation to the system."
+               << std::endl;
+
+  _q1_var = this->add_variable ("q1", _fe_order, _fe_family);
+
+  Parent::init_data();
+
+  for (auto & vec : _Fh)
+    vec->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
+}
+
+void AdvectionSystem::process_parameters_file (const std::string& parameters_filename)
+{
+  Parent::process_parameters_file(parameters_filename);
+
+  // First read in data from parameters_filename
+  GetPot infile(parameters_filename);
+  const Real u1_in = infile("u1", 1.);
+  const Real u2_in = infile("u2", 1.);
+  _u = Point(u1_in, u2_in, 0.);
+
+  std::string fe_order_str = infile("fe_order", std::string("CONSTANT"));
+  std::string fe_family_str = infile("fe_family", std::string("MONOMIAL"));
+
+  // Convert input strings to enumerations
+  _fe_order = Utility::string_to_enum<Order>(fe_order_str);
+  _fe_family = Utility::string_to_enum<FEFamily>(fe_family_str);
+}
+
+void AdvectionSystem::print_info ()
+{
+  Parent::print_info();
+
+  libMesh::out << std::endl << "==== AdvectionSystem ====" << std::endl;
+  libMesh::out << "u1 = " << _u(0) << ", u2 = " << _u(1) << std::endl;
+  libMesh::out << std::endl;
+}
+
+} // namespace libMesh
+
+
+

examples/transient/transient_ex3/advection_system.h

@@ -0,0 +1,113 @@
+#ifndef ADVECTION_SYSTEM_H
+#define ADVECTION_SYSTEM_H
+
+// Application includes
+#include "claw_system.h" // base class
+
+namespace libMesh
+{
+
+/**
+ * This class extends ClawSystem to implement pure advection
+ * in 2D.
+ *
+ * @author David J. Knezevic, 2012
+ * @author John W. Peterson, 2025 (modernization and libmesh example)
+ */
+class AdvectionSystem : public ClawSystem
+{
+public:
+
+  /**
+   * Constructor.  Optionally initializes required
+   * data structures.
+   */
+  AdvectionSystem (EquationSystems& es,
+                   const std::string& name,
+                   const unsigned int number);
+
+  /**
+   * Destructor. Defaulted out-of-line.
+   */
+  virtual ~AdvectionSystem ();
+
+  /**
+   * The type of system.
+   */
+  typedef AdvectionSystem sys_type;
+
+  /**
+   * @return a reference to ourself... I'm not sure what the point of
+   * this is. The LinearImplicitSystem base class has a similar API,
+   * but it's not virtual, so we aren't overriding it by doing
+   * this. This could probably be removed.
+   */
+  sys_type & system () { return *this; }
+
+  /**
+   * The type of the parent.
+   */
+  typedef ClawSystem Parent;
+
+  /**
+   * Right-hand side assembly. This function is called from the
+   * base class function, ClawSystem::solve_conservation_law().
+   * The input vector is either the initial condition (at the
+   * first timestep) or the solution vector from the previous
+   * timestep.
+   */
+  virtual void assemble_claw_rhs(NumericVector<Number> & q) override;
+
+  /**
+   * Initialize the system (e.g. add variables)
+   */
+  virtual void init_data() override;
+
+  /**
+   * Read in data from input file.
+   * Calls the base class version of this function and then reads in
+   * some additional parameters.
+   */
+  virtual void process_parameters_file (const std::string & parameters_filename) override;
+
+  /**
+   * Print out some info about the system's configuration.
+   */
+  virtual void print_info() override;
+
+protected:
+
+  /**
+   * Computes the vectors "uq" and "vq" where (u,v) are the (given,
+   * constant) advective velocity coefficients, and q is either the
+   * initial condition or the previous solution vector.
+   */
+  void update_Fh(NumericVector<Number> & q);
+
+  /**
+   * The advective flux vectors.  The outer std::vector length
+   * corresponds to the number of space dimensions in the problem.
+   */
+  std::vector<std::unique_ptr<NumericVector<Number>>> _Fh;
+
+  /**
+   * Variable number for q1
+   */
+  unsigned int _q1_var;
+
+  /**
+   * The (assumed constant) advection velocity
+   */
+  Point _u;
+
+  /**
+   * Store the FE Order and family specified by the input file.
+   * Defaults to (CONSTANT, MONOMIAL).
+   */
+  Order _fe_order;
+  FEFamily _fe_family;
+};
+
+} // namespace libMesh
+
+#endif