formatting and include desired qfunction

fem/dfem/doperator.hpp

@@ -968,12 +968,17 @@ void DifferentiableOperator::AddIntegrator(
          // Quadrature point local derivative cache for each element, with data
          // layout:
          // [test_vdim, test_op_dim, trial_vdim, trial_op_dim, qp, num_entities].
-         derivative_qp_caches[derivative_id] = Vector(test_vdim * test_op_dim *
-                                                      trial_vdim *
-                                                      total_trial_op_dim * num_qp * num_entities);
+         size_t derivative_qp_size = test_vdim * test_op_dim *
+                                     trial_vdim *
+                                     total_trial_op_dim * num_qp * num_entities;
+         if (derivative_qp_caches.count(derivative_qp_size) == 0)
+         {
+            derivative_qp_caches[derivative_qp_size] = Vector(derivative_qp_size);
+         }
+
          // Create local references for MSVC lambda capture compatibility
          auto& fields_ref = this->fields;
-         auto& derivative_qp_caches_ref = this->derivative_qp_caches[derivative_id];
+         auto& derivative_qp_caches_ref = this->derivative_qp_caches[derivative_qp_size];
 
          // In each of the callbacks we're saving the derivatives in the quadrature point
          // caches. This trades memory with computational effort but also minimizes

tests/unit/dfem/test_diffusion.cpp

@@ -69,9 +69,9 @@ template <int DIM> struct Diffusion
 };
 
 template <int DIM>
-void diffusion(const char *filename, int p)
+void diffusion(const char *filename, int p, bool include_mass)
 {
-   CAPTURE(filename, DIM, p);
+   CAPTURE(filename, DIM, p, include_mass);
 
    Mesh smesh(filename);
    ParMesh pmesh(MPI_COMM_WORLD, smesh);
@@ -116,6 +116,10 @@ void diffusion(const char *filename, int p)
 
    ParBilinearForm blf_fa(&pfes);
    blf_fa.AddDomainIntegrator(new DiffusionIntegrator(rho_coeff, ir));
+   if(include_mass)
+   {
+      blf_fa.AddDomainIntegrator(new MassIntegrator(ir));
+   }
    blf_fa.SetAssemblyLevel(AssemblyLevel::FULL);
    blf_fa.Assemble();
    blf_fa.Finalize();
@@ -130,6 +134,16 @@ void diffusion(const char *filename, int p)
    static constexpr int U = 0, Coords = 1, Rho = 3;
    const auto sol = std::vector{ FieldDescriptor{ U, &pfes } };
 
+   // Define mass qfunction in case it's used.
+      const auto mf_mass_qf =
+         [] MFEM_HOST_DEVICE(const dscalar_t u, const real_t &rho /*Needed for some reason*/,
+                             const tensor<real_t, DIM, DIM> &J, const real_t &w)
+      { return tuple{u * w * det(J)}; };
+      const auto desired_mass_qf =
+         [] MFEM_HOST_DEVICE(const dscalar_t u,
+                             const tensor<real_t, DIM, DIM> &J, const real_t &w)
+      { return tuple{u * w * det(J)}; };
+
    SECTION("action")
    {
       DOperator dop_mf(sol, {{Rho, &rho_ps}, {Coords, mfes}}, pmesh);
@@ -139,6 +153,21 @@ void diffusion(const char *filename, int p)
                                         Gradient<Coords>{}, Weight{} },
                                  tuple{ Gradient<U>{} }, *ir,
                                  all_domain_attr);
+      if(include_mass)
+      {
+         // This version works, but had to include unused Rho parameter in qfunction.
+         dop_mf.AddDomainIntegrator(mf_mass_qf,
+                                    tuple{ Value<U>{}, Identity<Rho>{},
+                                           Gradient<Coords>{}, Weight{} },
+                                    tuple{ Value<U>{} }, *ir,
+                                    all_domain_attr);
+         // Below version does not work.
+         // dop_mf.AddDomainIntegrator(desired_mass_qf,
+         //                            tuple{ Value<U>{},
+         //                                   Gradient<Coords>{}, Weight{} },
+         //                            tuple{ Value<U>{} }, *ir,
+         //                            all_domain_attr);
+      }
       dop_mf.SetParameters({ &rho_coeff_cv, nodes });
 
       pfes.GetRestrictionMatrix()->Mult(x, X);
@@ -159,6 +188,9 @@ void diffusion(const char *filename, int p)
 
    SECTION("action partial assembly")
    {
+      // is there a better way to skip a test than this?
+      if(!include_mass)
+      {
       static constexpr int QData = 2;
       UniformParameterSpace qd_ps(pmesh, *ir, DIM * DIM);
       ParameterFunction qdata(qd_ps);
@@ -196,6 +228,7 @@ void diffusion(const char *filename, int p)
 
       REQUIRE(norm_global == MFEM_Approx(0.0));
       MPI_Barrier(MPI_COMM_WORLD);
+      }
    }
 
    SECTION("action linearized")
@@ -208,6 +241,14 @@ void diffusion(const char *filename, int p)
                                         Gradient<Coords>{}, Weight{} },
                                  tuple{ Gradient<U>{} }, *ir,
                                  all_domain_attr, derivatives);
+      if(include_mass)
+      {
+         dop_mf.AddDomainIntegrator(mf_mass_qf,
+                                    tuple{ Value<U>{}, Identity<Rho>{} /*need this along with dummy argument to q function for some reason*/,
+                                           Gradient<Coords>{}, Weight{} },
+                                    tuple{ Value<U>{} }, *ir,
+                                    all_domain_attr, derivatives);
+      }
       dop_mf.SetParameters({ &rho_coeff_cv, nodes });
       auto dRdU = dop_mf.GetDerivative(U, {&x}, {&rho_coeff_cv, nodes});
 
@@ -258,6 +299,10 @@ void diffusion(const char *filename, int p)
          ConstantCoefficient one(1.0);
          ParBilinearForm vblf_fa(&vpfes);
          vblf_fa.AddDomainIntegrator(new VectorDiffusionIntegrator(one, ir));
+         // if(include_mass)
+         // {
+         //    vblf_fa.AddDomainIntegrator(new VectorMassIntegrator(one, ir));
+         // }
          vblf_fa.SetAssemblyLevel(AssemblyLevel::LEGACYFULL);
          vblf_fa.Assemble();
          vblf_fa.Finalize();
@@ -290,6 +335,14 @@ void diffusion(const char *filename, int p)
                                         Gradient<Coords>{}, Weight{} },
                                  tuple{ Gradient<U>{} }, *ir,
                                  all_domain_attr, derivatives);
+      if(include_mass)
+      {
+         dop_mf.AddDomainIntegrator(mf_mass_qf,
+                                    tuple{ Value<U>{}, Identity<Rho>{},
+                                           Gradient<Coords>{}, Weight{} },
+                                    tuple{ Value<U>{} }, *ir,
+                                    all_domain_attr, derivatives);
+      }
       dop_mf.SetParameters({ &rho_coeff_cv, nodes });
       auto dRdU = dop_mf.GetDerivative(U, {&x}, {&rho_coeff_cv, nodes});
 
@@ -306,6 +359,8 @@ TEST_CASE("dFEM Diffusion", "[Parallel][dFEM][GPU]")
 
    const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
 
+   const bool include_mass = GENERATE(false, true);
+
    SECTION("2d")
    {
       const auto filename2d =
@@ -316,7 +371,7 @@ TEST_CASE("dFEM Diffusion", "[Parallel][dFEM][GPU]")
             "../../data/inline-quad.mesh",
             "../../data/periodic-square.mesh"
          );
-      diffusion<2>(filename2d, p);
+      diffusion<2>(filename2d, p, include_mass);
    }
 
    SECTION("3d")
@@ -329,7 +384,7 @@ TEST_CASE("dFEM Diffusion", "[Parallel][dFEM][GPU]")
             "../../data/toroid-hex.mesh",
             "../../data/periodic-cube.mesh"
          );
-      diffusion<3>(filename3d, p);
+      diffusion<3>(filename3d, p, include_mass);
    }
 }