Updates to year on copyright statement on makefiles

examples/ex41p.cpp

@@ -59,7 +59,7 @@ public:
       : M(M_),
         K(K_),
         S(S_),
-	A(nullptr),
+        A(nullptr),
         linear_solver(M.GetComm()),
         dt(1.0)
    {
@@ -129,9 +129,9 @@ public:
    virtual
    ~IMEX_Evolution()
    {
-       delete dg_solver;
-       delete lor_solver;
-       delete M_prec;
+      delete dg_solver;
+      delete lor_solver;
+      delete M_prec;
    }
 
    virtual

linalg/ode.cpp

@@ -1491,4 +1491,4 @@ void IMEX_DIRK_RK3::Step(Vector &x, real_t &t, real_t &dt)
 
 }
 
-}
+}

linalg/ode.hpp

@@ -1025,4 +1025,4 @@ public:
 
 }
 
-#endif
+#endif

linalg/operator.cpp

@@ -934,4 +934,4 @@ real_t PowerMethod::EstimateLargestEigenvalue(Operator& opr, Vector& v0,
    return eigenvalue;
 }
 
-}
+}

linalg/operator.hpp

@@ -1207,4 +1207,4 @@ public:
 
 }
 
-#endif
+#endif

miniapps/CMakeLists.txt

@@ -36,6 +36,4 @@ add_subdirectory(parelag)
 add_subdirectory(tribol)
 add_subdirectory(hooke)
 add_subdirectory(dpg)
-add_subdirectory(hdiv-linear-solver)
-add_subdirectory(dfem)
-add_subdirectory(diag-smoothers)
+add_subdirectory(hdiv-linear-solver)

miniapps/mtop/CMakeLists.txt

@@ -10,22 +10,21 @@
 # CONTRIBUTING.md for details.
 
 list(APPEND SEQMTOP_COMMON_SOURCES
-	paramnonlinearform.cpp
-        mtop_integrators.cpp)
+	darcy_heat_transfer_ex.cpp)
 
-list(APPEND SEQMTOP_COMMON_HEADERS
-	paramnonlinearform.hpp
-        mtop_integrators.hpp)
+# list(APPEND SEQMTOP_COMMON_HEADERS
+# 	paramnonlinearform.hpp
+#         mtop_integrators.hpp)
 
 convert_filenames_to_full_paths(SEQMTOP_COMMON_SOURCES)
-convert_filenames_to_full_paths(SEQMTOP_COMMON_HEADERS)
+//convert_filenames_to_full_paths(SEQMTOP_COMMON_HEADERS)
 
 set(SEQMTOP_COMMON_FILES
 	EXTRA_SOURCES ${SEQMTOP_COMMON_SOURCES}
 	EXTRA_HEADERS ${SEQMTOP_COMMON_HEADERS})
 
 add_mfem_miniapp(seqheat
-	MAIN seqheat.cpp
+	MAIN darcy_heat_transfer_ex.cpp
 	${SEQMTOP_COMMON_FILES}
 	LIBRARIES mfem)
 
@@ -51,4 +50,4 @@ add_mfem_miniapp(parheat
      ${PARMTOP_COMMON_FILES}
      LIBRARIES mfem)
 
-endif ()
+endif ()

miniapps/mtop/darcy_heat_transfer_ex.cpp

@@ -0,0 +1,758 @@
+//                               MFEM Darcy Test Run
+//
+// Compile with: make darcy_heat_transfer_ex
+//
+//
+// Description:  This code performs the forward and backward adjoint solve for advection diffusion, where the velocity field is given by Darcy
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+
+using namespace std;
+using namespace mfem;
+
+
+// *****Funtion definitions for the Advection-Diffusion solve******
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double theta0_function(const Vector &x);
+
+// true solution
+real_t theta_exact(const Vector &x, real_t t);
+
+// rhs
+double forcing_function(const Vector &x, real_t t);
+
+// inflow
+double inflow_function(const Vector &x);
+
+real_t f_natural(const Vector & x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+class DG_Solver : public Solver
+{
+private:
+   SparseMatrix &M, &K, &S, A;
+   CGSolver linear_solver;
+   BlockILU prec;
+   real_t dt;
+public:
+   DG_Solver(SparseMatrix &M_, SparseMatrix &K_, SparseMatrix &S_,
+             const FiniteElementSpace &fes)
+      : M(M_),
+        K(K_),
+        S(S_),
+        prec(fes.GetTypicalFE()->GetDof(),
+             BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
+        dt(1.0)
+   {
+      linear_solver.iterative_mode = false;
+      linear_solver.SetRelTol(1e-9);
+      linear_solver.SetAbsTol(0.0);
+      linear_solver.SetMaxIter(100);
+      linear_solver.SetPrintLevel(0);
+      linear_solver.SetPreconditioner(prec);
+   }
+
+   void SetTimeStep(real_t dt_)
+   {
+      if (dt_ != dt)
+      {
+         dt = dt_;
+         // Form operator A = M + dt*S
+         A = S;
+         A *= dt;
+         A += M;
+
+         // this will also call SetOperator on the preconditioner
+         linear_solver.SetOperator(A);
+      }
+   }
+
+   void SetOperator(const Operator &op) override
+   {
+      linear_solver.SetOperator(op);
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      linear_solver.Mult(x, y);
+   }
+};
+
+
+/** A time-dependent operator for the right-hand side of the ODE. The DG weak
+    form of the advection-diffusion equation is (M + dt S) du/dt = Su - K u + b, where M and K are the mass
+    and advection matrices, and b describes the flow on the boundary. In the case of IMEX evolution, the diffusion term is treated
+    implicitly, and the advection term is treated explicitly.  */
+class IMEX_Evolution : public SplitTimeDependentOperator
+{
+private:
+   BilinearForm &M, &K, &S;
+   const Vector &b;
+   unique_ptr<Solver> M_prec;
+   CGSolver M_solver;
+   unique_ptr<DG_Solver> dg_solver;
+
+   mutable Vector z;
+
+public:
+   IMEX_Evolution(BilinearForm &M_, BilinearForm &K_, BilinearForm &S_,
+                  const Vector &b_);
+
+   void Mult1(const Vector &x, Vector &y) const;
+   void ImplicitSolve2(const real_t dt, const Vector &x, Vector &k) override;
+};
+
+// *****Define the analytical solution and forcing terms / boundary conditions for Darcy*****
+void uFun_ex(const Vector & x, Vector & u);
+real_t pFun_ex(const Vector & x);
+void fFun(const Vector & x, Vector & f);
+real_t gFun(const Vector & x);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file =
+      "square-extended.mesh"; //reference square, but extended to be [-1, 1] x [-1, 1]
+   int order_darcy = 1;
+   int ref_levels = 2;
+   int order_ad = 3;
+   int ode_solver_type = 55;
+   double t_final = 10.0;
+   double d_coef = 0.01;
+   double dt = 0.01;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   bool visualization = true;
+   bool visit = false;
+   bool binary = false;
+   int vis_steps = 5;
+   bool paraview = false;
+   int precision = 16;
+   const char *device_config = "cpu";
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order_darcy, "-od", "--order_darcy",
+                  "Order (degree) of the finite elements for darcy solve.");
+   args.AddOption(&order_ad, "-oad", "--order_ad",
+                  "Order (degree) of the finite elements for advection diffusion.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "55 - Forward Backward Euler, 56 - IMEXRK2(2,2,2), 57 - IMEXRK2(2,3,2), 58 - IMEX_DIRK_RK3\n");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&d_coef, "-d", "--diff-coef",
+                  "Diffusion coefficient.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the two DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the two DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
+                  "--no-paraview-datafiles",
+                  "Save data files for ParaView (paraview.org) visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order_ad+1)*(order_ad+1);
+   }
+   args.PrintOptions(cout);
+   Device device(device_config);
+   device.Print();
+
+   // 2. Define the ODE solver used for time integration. Several explicit, implicit and IMEX
+   //    Runge-Kutta methods are available.
+   unique_ptr<SplitODESolver> ode_solver = SplitODESolver::Select(ode_solver_type);
+   unique_ptr<SplitODESolver> ode_solver_adj = SplitODESolver::Select(
+                                                  ode_solver_type);
+   // 3. Read the mesh from the given mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 4. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter.
+   for (int lev = 0; lev < ref_levels; lev++) {mesh.UniformRefinement();}
+   if (mesh.NURBSext) {mesh.SetCurvature(max(order_ad, 1));}
+   mesh.GetBoundingBox(bb_min, bb_max, max(order_ad, 1));
+
+   // ********DARCY SOLVE
+   // 5. Define a finite element space on the mesh. Here we use the
+   //    Raviart-Thomas finite elements of the specified order.
+   FiniteElementCollection *hdiv_coll(new RT_FECollection(order_darcy, dim));
+   FiniteElementCollection *l2_coll(new L2_FECollection(order_darcy, dim));
+
+   FiniteElementSpace *R_space = new FiniteElementSpace(&mesh, hdiv_coll);
+   FiniteElementSpace *W_space = new FiniteElementSpace(&mesh, l2_coll);
+
+   // 6. Define the BlockStructure of the problem, i.e. define the array of
+   //    offsets for each variable. The last component of the Array is the sum
+   //    of the dimensions of each block.
+   Array<int> block_offsets(3); // number of variables + 1
+   block_offsets[0] = 0;
+   block_offsets[1] = R_space->GetVSize();
+   block_offsets[2] = W_space->GetVSize();
+   block_offsets.PartialSum();
+
+   std::cout << "***********************************************************\n";
+   std::cout << "dim(R) = " << block_offsets[1] - block_offsets[0] << "\n";
+   std::cout << "dim(W) = " << block_offsets[2] - block_offsets[1] << "\n";
+   std::cout << "dim(R+W) = " << block_offsets.Last() << "\n";
+   std::cout << "***********************************************************\n";
+
+   // 7. Define the coefficients, analytical solution, and rhs of the Darcy PDE.
+   ConstantCoefficient one(1.0);
+
+   VectorFunctionCoefficient fcoeff(dim, fFun);
+   FunctionCoefficient fnatcoeff(f_natural);
+   FunctionCoefficient gcoeff(gFun);
+
+   VectorFunctionCoefficient ucoeff(dim, uFun_ex);
+   FunctionCoefficient pcoeff(pFun_ex);
+
+   // 8. Allocate memory for solution and rhs of Darcy
+   MemoryType mt = device.GetMemoryType();
+   BlockVector x(block_offsets, mt), rhs(block_offsets, mt);
+
+   LinearForm *fform(new LinearForm);
+   fform->Update(R_space, rhs.GetBlock(0), 0);
+   fform->AddDomainIntegrator(new VectorFEDomainLFIntegrator(fcoeff));
+   fform->AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(fnatcoeff));
+   fform->Assemble();
+   fform->SyncAliasMemory(rhs);
+
+   LinearForm *gform(new LinearForm);
+   gform->Update(W_space, rhs.GetBlock(1), 0);
+   gform->AddDomainIntegrator(new DomainLFIntegrator(gcoeff));
+   gform->Assemble();
+   gform->SyncAliasMemory(rhs);
+
+   // 9. Assemble the finite element matrices for the Darcy operator
+   //
+   //                            D = [ M  B^T ]
+   //                                [ B   0  ]
+   //     where:
+   //
+   //     M = \int_\Omega k u_h \cdot v_h d\Omega   u_h, v_h \in R_h
+   //     B   = -\int_\Omega \div u_h q_h d\Omega   u_h \in R_h, q_h \in W_h
+   BilinearForm *mVarf(new BilinearForm(R_space));
+   mVarf->AddDomainIntegrator(new VectorFEMassIntegrator(one));
+   mVarf->Assemble();
+   MixedBilinearForm *bVarf(new MixedBilinearForm(R_space, W_space));
+   bVarf->AddDomainIntegrator(new VectorFEDivergenceIntegrator);
+   bVarf->Assemble();
+   mVarf->Finalize();
+   bVarf->Finalize();
+   BlockOperator darcyOp(block_offsets);
+   TransposeOperator *Bt = NULL;
+   SparseMatrix &M(mVarf->SpMat());
+   SparseMatrix &B(bVarf->SpMat());
+   B *= -1.;
+   Bt = new TransposeOperator(&B);
+   darcyOp.SetBlock(0,0, &M);
+   darcyOp.SetBlock(0,1, Bt);
+   darcyOp.SetBlock(1,0, &B);
+
+   // 10. Construct the operators for preconditioner
+   //
+   //                 P = [ diag(M)         0         ]
+   //                     [  0       B diag(M)^-1 B^T ]
+   //
+   //     Here we use Symmetric Gauss-Seidel to approximate the inverse of the
+   //     pressure Schur Complement
+   SparseMatrix *MinvBt = NULL;
+   Vector Md(mVarf->Height());
+
+   BlockDiagonalPreconditioner darcyPrec(block_offsets);
+   Solver *invM, *invS;
+   SparseMatrix *S = NULL;
+   // SparseMatrix &M(mVarf->SpMat());
+   M.GetDiag(Md);
+   Md.HostReadWrite();
+   // SparseMatrix &B(bVarf->SpMat());
+   MinvBt = Transpose(B);
+   for (int i = 0; i < Md.Size(); i++)
+   {
+      MinvBt->ScaleRow(i, 1./Md(i));
+   }
+   S = Mult(B, *MinvBt);
+   invM = new DSmoother(M);
+#ifndef MFEM_USE_SUITESPARSE
+   invS = new GSSmoother(*S);
+#else
+   invS = new UMFPackSolver(*S);
+#endif
+   invM->iterative_mode = false;
+   invS->iterative_mode = false;
+
+   darcyPrec.SetDiagonalBlock(0, invM);
+   darcyPrec.SetDiagonalBlock(1, invS);
+
+   // 11. Solve the linear system with MINRES.
+   //     Check the norm of the unpreconditioned residual.
+   int maxIter(1000);
+   real_t rtol(1.e-6);
+   real_t atol(1.e-10);
+
+   MINRESSolver solver;
+   solver.SetAbsTol(atol);
+   solver.SetRelTol(rtol);
+   solver.SetMaxIter(maxIter);
+   solver.SetOperator(darcyOp);
+   solver.SetPreconditioner(darcyPrec);
+   solver.SetPrintLevel(1);
+   x = 0.0;
+   solver.Mult(rhs, x);
+
+   if (solver.GetConverged())
+   {
+      std::cout << "MINRES converged in " << solver.GetNumIterations()
+                << " iterations with a residual norm of "
+                << solver.GetFinalNorm() << ".\n";
+   }
+   else
+   {
+      std::cout << "MINRES did not converge in " << solver.GetNumIterations()
+                << " iterations. Residual norm is " << solver.GetFinalNorm()
+                << ".\n";
+   }
+
+   // 12. Create the grid functions u and p. Compute the L2 error norms.
+   GridFunction u, p;
+   u.MakeRef(R_space, x.GetBlock(0), 0);
+   p.MakeRef(W_space, x.GetBlock(1), 0);
+
+   int order_quad = max(2, 2*order_darcy+1);
+   const IntegrationRule *irs[Geometry::NumGeom];
+   for (int i=0; i < Geometry::NumGeom; ++i)
+   {
+      irs[i] = &(IntRules.Get(i, order_quad));
+   }
+
+   real_t err_u  = u.ComputeL2Error(ucoeff, irs);
+   real_t norm_u = ComputeLpNorm(2., ucoeff, mesh, irs);
+   real_t err_p  = p.ComputeL2Error(pcoeff, irs);
+   real_t norm_p = ComputeLpNorm(2., pcoeff, mesh, irs);
+
+   std::cout << "|| u_h - u_ex || / || u_ex || = " << err_u / norm_u << "\n";
+   std::cout << "|| p_h - p_ex || / || p_ex || = " << err_p / norm_p << "\n";
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream u_sock(vishost, visport);
+      u_sock.precision(8);
+      u_sock << "solution\n" << mesh << u << "window_title 'Velocity'" << endl;
+      socketstream p_sock(vishost, visport);
+      p_sock.precision(8);
+      p_sock << "solution\n" << mesh << p << "window_title 'Pressure'" << endl;
+   }
+
+
+   // ******Forward Advection-Diffusion solve
+   // 13. Define the DG finite element space on the
+   //    refined mesh of the given polynomial order.
+   DG_FECollection fec(order_ad, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec);
+   int num_dofs = fes.GetNDofs();
+
+   cout << "Number of unknowns (advection diffusion problem): " << fes.GetVSize()
+        << endl;
+
+   // 14. Set up and assemble the parallel bilinear and linear forms (and the
+   //    parallel hypre matrices) corresponding to the DG discretization. The
+   //    DGTraceIntegrator involves integrals over mesh interior faces.
+   const GridFunction* u_pointer = &u;
+   VectorGridFunctionCoefficient velocity(u_pointer);
+   FunctionCoefficient inflow(inflow_function);
+   ConstantCoefficient diff_coef(d_coef);
+
+   BilinearForm m(&fes);
+   m.AddDomainIntegrator(new MassIntegrator);
+
+   BilinearForm k(&fes);
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k.AddInteriorFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity,
+                                                                    -1.0));
+   k.AddBdrFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity, -1.0));
+
+   BilinearForm s(&fes);
+   s.AddDomainIntegrator(new DiffusionIntegrator(diff_coef));
+   s.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
+   s.AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(new BoundaryFlowIntegrator(inflow, velocity, -1.0));
+   //b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(U, diff_coef, sigma, kappa));
+
+   int skip_zeros = 0;
+   m.Assemble(skip_zeros);
+   k.Assemble(skip_zeros);
+   s.Assemble(skip_zeros);
+   b.Assemble();
+
+   m.Finalize(skip_zeros);
+   k.Finalize(skip_zeros);
+   s.Finalize(skip_zeros);
+
+   // 15. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   FunctionCoefficient theta0(theta0_function);
+   GridFunction theta(&fes);
+   theta.ProjectCoefficient(theta0);
+
+   // Set up visualization, if desired.
+   ParaViewDataCollection *pd_forward = NULL;
+   if (paraview)
+   {
+      pd_forward = new ParaViewDataCollection("darcy-adv-diff-forward", &mesh);
+      pd_forward->SetPrefixPath("ParaView");
+      pd_forward->RegisterField("solution_forward", &theta);
+      pd_forward->SetLevelsOfDetail(order_ad);
+      pd_forward->SetDataFormat(VTKFormat::BINARY);
+      pd_forward->SetHighOrderOutput(true);
+      pd_forward->SetCycle(0);
+      pd_forward->SetTime(0.0);
+      pd_forward->Save();
+   }
+
+   // 16. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   IMEX_Evolution adv(m, k, s, b);
+
+   real_t t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+   int n_steps = (int)ceil(t_final / dt);
+   double dt_real = t_final / n_steps;
+   // Vector err_vec(n_steps-1);
+
+   std::vector<GridFunction> theta_gf_vector;
+   theta_gf_vector.push_back(theta);
+
+   for (int ti = 0; ti < n_steps; ti++)
+   {
+      ode_solver->Step(theta, t, dt_real);
+      theta_gf_vector.push_back(theta);
+      if (ti % vis_steps == 0 || ti == n_steps -1)
+      {
+         cout << "time step: " << ti << ", time: " << t << endl;
+         if (paraview)
+         {
+            pd_forward->SetCycle(ti);
+            pd_forward->SetTime(t);
+            pd_forward->Save();
+         }
+      }
+   }
+
+   // ******Backward Advection-Diffusion solve
+   // 17. Define the DG finite element space on the
+   //    refined mesh of the given polynomial order.
+   DG_FECollection fec_adjoint(order_ad, dim);
+   FiniteElementSpace fes_adjoint(&mesh, &fec_adjoint);
+
+   // 18. Set up and assemble the parallel bilinear and linear forms (and the
+   //    parallel hypre matrices) corresponding to the DG discretization. The
+   //    DGTraceIntegrator involves integrals over mesh interior faces.
+   ConstantCoefficient zero(0.0);
+   GridFunctionCoefficient theta_coeff(&(theta_gf_vector[n_steps-1]));
+   FunctionCoefficient inflow_adj(inflow_function); //zero for now
+   ConstantCoefficient diff_coef_adj(-d_coef);
+
+   // FunctionCoefficient theta_exact_coeff(theta_exact);
+   BilinearForm m_adj(&fes_adjoint);
+   m_adj.AddDomainIntegrator(new MassIntegrator);
+
+   BilinearForm k_adj(&fes_adjoint);
+   k_adj.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k_adj.AddInteriorFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity,
+                                                                        -1.0));
+   k_adj.AddBdrFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity,
+                                                                   -1.0));
+
+   BilinearForm s_adj(&fes_adjoint);
+   s_adj.AddDomainIntegrator(new DiffusionIntegrator(diff_coef_adj));
+   s_adj.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef_adj, sigma,
+                                                             kappa));
+   s_adj.AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef_adj, sigma,
+                                                        kappa));
+
+   LinearForm b_adj(&fes_adjoint);
+   b_adj.AddDomainIntegrator(new DomainLFIntegrator(theta_coeff));
+   //b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(zero, diff_coef, sigma, kappa));
+
+   //int skip_zeros = 0;
+   m_adj.Assemble(skip_zeros);
+   m_adj.Finalize(skip_zeros);
+   k_adj.Assemble(skip_zeros);
+   k_adj.Finalize(skip_zeros);
+   s_adj.Assemble(skip_zeros);
+   s_adj.Finalize(skip_zeros);
+   b_adj.Assemble();
+
+   // 19. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction lam(&fes_adjoint);
+   lam.ProjectCoefficient(zero);
+   ParaViewDataCollection *pd_backward = NULL;
+   if (paraview)
+   {
+      pd_backward = new ParaViewDataCollection("darcy-adv-diff-backward", &mesh);
+      pd_backward->SetPrefixPath("ParaView");
+      pd_backward->RegisterField("solution-backward", &lam);
+      pd_backward->SetLevelsOfDetail(order_ad);
+      pd_backward->SetDataFormat(VTKFormat::BINARY);
+      pd_backward->SetHighOrderOutput(true);
+      pd_backward->SetCycle(0);
+      pd_backward->SetTime(t_final);
+      pd_backward->Save();
+   }
+
+   // 20. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   IMEX_Evolution adv_adj(m_adj, k_adj, s_adj, b_adj);
+
+   real_t t_adj = t_final;
+   adv_adj.SetTime(t_adj);
+   ode_solver_adj->Init(adv_adj);
+
+   // int n_steps = (int)ceil(t_final / dt);
+   double dt_real_adj = -dt;
+   std::cout << "dt back = " << dt_real_adj << std::endl;
+   //Vector err_vec(n_steps-1);
+
+   for (int ti = 0; ti < n_steps; ti++)
+   {
+      ode_solver_adj->Step(lam, t_adj, dt_real_adj);
+      Vector lam_vals(num_dofs);
+      Vector theta_values(num_dofs);
+      const GridFunction* theta_gf = theta_coeff.GetGridFunction();
+      theta_gf->GetTrueDofs(theta_values);
+      lam.GetTrueDofs(lam_vals);
+      theta_coeff = *(new GridFunctionCoefficient(&(theta_gf_vector[n_steps - ti -
+                                                                            1])));
+      b_adj = *(new LinearForm(&fes_adjoint));
+      b_adj.AddDomainIntegrator(new DomainLFIntegrator(theta_coeff));
+      b_adj.Assemble();
+      if (ti % vis_steps == 0 || ti == n_steps - 1)
+      {
+         cout << "time step: " << ti << ", time: " << t_adj << endl;
+         if (paraview)
+         {
+            pd_backward->SetCycle(ti);
+            pd_backward->SetTime(t_adj);
+            pd_backward->Save();
+         }
+      }
+   }
+
+   // 21. Free the used memory.
+   // delete &ode_solver;
+   // delete &adv;
+   // delete &adv_adj;
+   delete fform;
+   delete gform;
+   delete invM;
+   delete invS;
+   delete S;
+   delete Bt;
+   delete MinvBt;
+   delete mVarf;
+   delete bVarf;
+   delete W_space;
+   delete R_space;
+   delete l2_coll;
+   delete hdiv_coll;
+   // delete &b_adj;
+   // delete &theta_coeff;
+
+   return 0;
+}
+
+void uFun_ex(const Vector & x, Vector & u)
+{
+   real_t xi(x(0));
+   real_t yi(x(1));
+   real_t zi(0.0);
+   if (x.Size() == 3)
+   {
+      zi = x(2);
+   }
+
+   u(0) = - exp(xi)*sin(yi)*cos(zi);
+   u(1) = - exp(xi)*cos(yi)*cos(zi);
+
+   if (x.Size() == 3)
+   {
+      u(2) = exp(xi)*sin(yi)*sin(zi);
+   }
+}
+
+// Change if needed
+real_t pFun_ex(const Vector & x)
+{
+   real_t xi(x(0));
+   real_t yi(x(1));
+   real_t zi(0.0);
+
+   if (x.Size() == 3)
+   {
+      zi = x(2);
+   }
+
+   return exp(xi)*sin(yi)*cos(zi);
+}
+
+void fFun(const Vector & x, Vector & f)
+{
+   f = 0.0;
+}
+
+real_t gFun(const Vector & x)
+{
+   if (x.Size() == 3)
+   {
+      return -pFun_ex(x);
+   }
+   else
+   {
+      return 0;
+   }
+}
+
+real_t f_natural(const Vector & x)
+{
+   return (-pFun_ex(x));
+}
+
+// Implementation of class IMEX_Evolution
+IMEX_Evolution::IMEX_Evolution(BilinearForm &M_, BilinearForm &K_,
+                               BilinearForm &S_, const Vector &b_)
+   : SplitTimeDependentOperator(M_.FESpace()->GetTrueVSize()),
+     M(M_), K(K_), S(S_), b(b_), z(height)
+{
+   Array<int> ess_tdof_list;
+   if (M.GetAssemblyLevel() == AssemblyLevel::LEGACY)
+   {
+      M_prec = make_unique<DSmoother>(M.SpMat());
+      M_solver.SetOperator(M.SpMat());
+      dg_solver = make_unique<DG_Solver>(M.SpMat(), K.SpMat(), S.SpMat(),
+                                         *M.FESpace());
+   }
+   else
+   {
+      M_prec = make_unique<OperatorJacobiSmoother>(M, ess_tdof_list);
+      M_solver.SetOperator(M);
+      dg_solver = NULL;
+   }
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IMEX_Evolution::Mult1(const Vector &x, Vector &y) const
+{
+   // Perform the explicit step
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ImplicitSolve2(const real_t dt, const Vector &x, Vector &k)
+{
+   // Perform the implicit step
+   // solve for k, k = -(M+dt S)^{-1} S x
+   MFEM_VERIFY(dg_solver != NULL,
+               "Implicit time integration is not supported with partial assembly");
+   S.Mult(x, z);
+   z*= -1.0;
+   dg_solver->SetTimeStep(dt);
+   dg_solver->Mult(z, k);
+}
+
+
+// Initial condition
+double theta0_function(const Vector &x)
+{
+   int dim = x.Size();
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   // for (int i = 0; i < dim; i++)
+   // {
+   //    double center = (bb_min[i] + bb_max[i]) * 0.5;
+   //    X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   // }
+
+   double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+   if (dim == 3)
+   {
+      const double s = (1. + 0.25*cos(2*M_PI*x(2)));
+      rx *= s;
+      ry *= s;
+   }
+   return ( erfc(w*(x(0)-cx-rx))*erfc(-w*(x(0)-cx+rx))*erfc(w*(x(1)-cy-ry))*erfc(
+               -w*(x(1)-cy+ry)) )/16;
+}
+
+//forcing term
+real_t forcing_function(const Vector &x, real_t t)
+{
+   int dim = x.Size();
+   //map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+   return 0.0;
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x)
+{
+   return 0.0;
+}

miniapps/mtop/makefile

@@ -12,10 +12,11 @@
 # Use the MFEM build directory
 MFEM_DIR ?= ../..
 MFEM_BUILD_DIR ?= ../..
-MFEM_INSTALL_DIR ?= ../../mfem
 SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/mtop/,)
-CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
-   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
 
 # Include defaults.mk to get XLINKER
 DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
@@ -24,11 +25,10 @@ include $(DEFAULTS_MK)
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-MTOP_COMMON_SRC = mtop_integrators.cpp paramnonlinearform.cpp pparamnonlinearform.cpp
-
+MTOP_COMMON_SRC =  darcy_heat_transfer_ex.cpp
 MTOP_COMMON_OBJ = $(MTOP_COMMON_SRC:.cpp=.o)
 
-SEQ_MINIAPPS = seqheat
+SEQ_MINIAPPS = seqheat darcy_heat_transfer_ex
 PAR_MINIAPPS = parheat
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
@@ -76,4 +76,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -rf SeqHeat* ParHeat*
+	@rm -rf SeqHeat* ParHeat* for_adv_diff_solve* darcy_heat_transfer_ex*

miniapps/mtop/mtop_integrators.cpp

@@ -1,390 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "mtop_integrators.hpp"
-
-namespace mfem
-{
-
-real_t ParametricLinearDiffusion::GetElementEnergy(const
-                                                   Array<const FiniteElement *> &el,
-                                                   const Array<const FiniteElement *> &pel,
-                                                   ElementTransformation &Tr,
-                                                   const Array<const Vector *> &elfun,
-                                                   const Array<const Vector *> &pelfun)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dof_r0 = pel[0]->GetDof();
-
-   int dim = el[0]->GetDim();
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::GetElementEnergy"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   B=0.0;
-
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-
-   real_t energy =0.0;
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                 +pel[0]->GetOrder();
-      ir=&IntRules.Get(Tr.GetGeometryType(),order);
-   }
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w=Tr.Weight();
-      w = ip.weight * w;
-
-      el[0]->CalcPhysDShape(Tr,dsu0);
-      el[0]->CalcPhysShape(Tr,shu0);
-      pel[0]->CalcPhysShape(Tr,shr0);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // set the matrix B
-      for (int jj=0; jj<dim; jj++)
-      {
-         B.SetCol(jj,dsu0.GetColumn(jj));
-      }
-      B.SetCol(3,shu0);
-      B.MultTranspose(*elfun[0],uu);
-      energy=energy+w * qfun.QEnergy(Tr,ip,param,uu);
-   }
-   return energy;
-}
-
-
-void ParametricLinearDiffusion::AssembleElementVector(const
-                                                      Array<const FiniteElement *> &el,
-                                                      const Array<const FiniteElement *> &pel,
-                                                      ElementTransformation &Tr,
-                                                      const Array<const Vector *> &elfun,
-                                                      const Array<const Vector *> &pelfun,
-                                                      const Array<Vector *> &elvec)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dof_r0 = pel[0]->GetDof();
-
-   int dim = el[0]->GetDim();
-
-   elvec[0]->SetSize(dof_u0);
-   *elvec[0]=0.0;
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::AssembleElementVector"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   B=0.0;
-
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-   Vector rr(4);
-   Vector lvec; lvec.SetSize(dof_u0);
-
-   const IntegrationRule *ir = nullptr;
-   int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-              +pel[0]->GetOrder();
-   ir=&IntRules.Get(Tr.GetGeometryType(),order);
-
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w=Tr.Weight();
-      w = ip.weight * w;
-
-      el[0]->CalcPhysDShape(Tr,dsu0);
-      el[0]->CalcPhysShape(Tr,shu0);
-      pel[0]->CalcPhysShape(Tr,shr0);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // set the matrix B
-      for (int jj=0; jj<dim; jj++)
-      {
-         B.SetCol(jj,dsu0.GetColumn(jj));
-      }
-      B.SetCol(3,shu0);
-      B.MultTranspose(*elfun[0],uu);
-      qfun.QResidual(Tr,ip,param, uu, rr);
-
-      B.Mult(rr,lvec);
-      elvec[0]->Add(w,lvec);
-   }
-}
-
-void ParametricLinearDiffusion::AssembleElementGrad(const
-                                                    Array<const FiniteElement *> &el,
-                                                    const Array<const FiniteElement *> &pel,
-                                                    ElementTransformation &Tr,
-                                                    const Array<const Vector *> &elfun,
-                                                    const Array<const Vector *> &pelfun,
-                                                    const Array2D<DenseMatrix *> &elmats)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dof_r0 = pel[0]->GetDof();
-
-   int dim = el[0]->GetDim();
-
-   DenseMatrix* K=elmats(0,0);
-   K->SetSize(dof_u0,dof_u0);
-   (*K)=0.0;
-
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::AssembleElementGrad"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   DenseMatrix A(dof_u0, 4);
-   B=0.0;
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-   DenseMatrix hh(4,4);
-   Vector lvec; lvec.SetSize(dof_u0);
-
-   const IntegrationRule *ir = nullptr;
-   int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-              +pel[0]->GetOrder();
-   ir=&IntRules.Get(Tr.GetGeometryType(),order);
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w = Tr.Weight();
-      w = ip.weight * w;
-
-      el[0]->CalcPhysDShape(Tr,dsu0);
-      el[0]->CalcPhysShape(Tr,shu0);
-      pel[0]->CalcPhysShape(Tr,shr0);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // set the matrix B
-      for (int jj=0; jj<dim; jj++)
-      {
-         B.SetCol(jj,dsu0.GetColumn(jj));
-      }
-      B.SetCol(3,shu0);
-      B.MultTranspose(*elfun[0],uu);
-      qfun.QGradResidual(Tr,ip,param,uu,hh);
-      Mult(B,hh,A);
-      AddMult_a_ABt(w,A,B,*K);
-   }
-}
-
-void ParametricLinearDiffusion::AssemblePrmElementVector(
-   const Array<const FiniteElement *> &el,
-   const Array<const FiniteElement *> &pel,
-   ElementTransformation &Tr,
-   const Array<const Vector *> &elfun,
-   const Array<const Vector *> &alfun,
-   const Array<const Vector *> &pelfun,
-   const Array<Vector *> &elvec)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dof_r0 = pel[0]->GetDof();
-
-   int dim = el[0]->GetDim();
-   Vector& e0 = *(elvec[0]);
-
-   e0.SetSize(dof_r0);
-   e0=0.0;
-
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::AssemblePrmElementVector"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   B=0.0;
-
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-   Vector aa(4); aa=0.0;
-   Vector rr(1);
-   Vector lvec0; lvec0.SetSize(dof_r0);
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                 +pel[0]->GetOrder();
-      ir=&IntRules.Get(Tr.GetGeometryType(),order);
-   }
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w=Tr.Weight();
-      w = ip.weight * w;
-
-      el[0]->CalcPhysDShape(Tr,dsu0);
-      el[0]->CalcPhysShape(Tr,shu0);
-      pel[0]->CalcPhysShape(Tr,shr0);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // set the matrix B
-      for (int jj=0; jj<dim; jj++)
-      {
-         B.SetCol(jj,dsu0.GetColumn(jj));
-      }
-      B.SetCol(3,shu0);
-      B.MultTranspose(*elfun[0],uu);
-      B.MultTranspose(*alfun[0],aa);
-
-      qfun.AQResidual(Tr, ip, param, uu, aa, rr);
-
-      lvec0=shr0;
-      lvec0*=rr[0];
-
-      e0.Add(w,lvec0);
-   }
-}
-
-real_t DiffusionObjIntegrator::GetElementEnergy(const
-                                                Array<const FiniteElement *> &el,
-                                                ElementTransformation &Tr,
-                                                const Array<const Vector *> &elfun)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dim = el[0]->GetDim();
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("DiffusionObjIntegrator::GetElementEnergy"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-
-   real_t w;
-   real_t val;
-
-   real_t energy = 0.0;
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
-      ir=&IntRules.Get(Tr.GetGeometryType(),order);
-   }
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w=Tr.Weight();
-
-      w = ip.weight * w;
-
-      el[0]->CalcPhysShape(Tr,shu0);
-
-      val=shu0*(*elfun[0]);
-      energy=energy + w * val * val;
-   }
-   return 0.5*energy;
-}
-
-void DiffusionObjIntegrator::AssembleElementVector(const
-                                                   Array<const FiniteElement *> &el,
-                                                   ElementTransformation &Tr,
-                                                   const Array<const Vector *> &elfun,
-                                                   const Array<Vector *> &elvec)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dim = el[0]->GetDim();
-   int spaceDim = Tr.GetSpaceDim();
-
-   elvec[0]->SetSize(dof_u0);
-   *elvec[0]=0.0;
-
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("DiffusionObjIntegrator::GetElementEnergy"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-
-   real_t w;
-   real_t val;
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
-      ir=&IntRules.Get(Tr.GetGeometryType(),order);
-   }
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w=Tr.Weight();
-
-      w = ip.weight * w;
-
-      el[0]->CalcPhysShape(Tr,shu0);
-
-      val=shu0*(*elfun[0]);
-
-      elvec[0]->Add(w*val,shu0);
-   }
-}
-
-
-} // end mfem namespace

miniapps/mtop/mtop_integrators.hpp

@@ -1,233 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MTOPINTEGRATORS_HPP
-#define MTOPINTEGRATORS_HPP
-
-#include "mfem.hpp"
-#include "paramnonlinearform.hpp"
-
-#include <map>
-
-namespace mfem
-{
-
-/// Base class for representing function at integration points.
-class BaseQFunction
-{
-public:
-   virtual ~BaseQFunction() {}
-
-   /// Returns a user defined string identifying the function.
-   virtual std::string GetType()=0;
-
-   // Returns the energy at an integration point.
-   virtual
-   real_t QEnergy(ElementTransformation &T, const IntegrationPoint &ip,
-                  mfem::Vector &dd, mfem::Vector &uu)
-   {
-      return 0.0;
-   }
-
-   // Returns the residual at an integration point.
-   virtual
-   void QResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                  mfem::Vector &dd, mfem::Vector &uu, mfem::Vector &rr)=0;
-
-   /// Returns the gradient of the residual at a integration point.
-   virtual
-   void QGradResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                      mfem::Vector &dd, mfem::Vector &uu, mfem::DenseMatrix &hh)=0;
-
-   /// Returns the gradient of the residual with respect to the design
-   /// parameters, multiplied by the adjoint.
-   virtual
-   void AQResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                   mfem::Vector &dd, mfem::Vector &uu,
-                   mfem::Vector &aa, mfem::Vector &rr)=0;
-
-};
-
-/* QLinearDiffusion implements methods for computing the energy, the residual,
- * gradient of the residual and the product of the adjoint fields with the
- * derivative of the residual with respect to the parameters. All computations
- * are performed at a integration point. Therefore the vectors (vv,uu,aa,rr ..)
- * hold the fields' values and the fields' derivatives at the integration
- * point. For example for a single scalar parametric field representing the
- * density in topology optimization the vector dd will have size one and the
- * element will be the density at the integration point. The map between state
- * and parameter is not fixed and depends on the implementation of the QFunction
- * class. */
-class QLinearDiffusion:public BaseQFunction
-{
-public:
-   QLinearDiffusion(mfem::Coefficient& diffco, mfem::Coefficient& hsrco,
-                    real_t pp=1.0, real_t minrho=1e-7, real_t betac=4.0, real_t etac=0.5):
-      diff(diffco),load(hsrco), powerc(pp), rhomin(minrho), beta(betac), eta(etac)
-   {
-
-   }
-
-   std::string GetType() override
-   {
-      return "QLinearDiffusion";
-   }
-
-   real_t QEnergy(ElementTransformation &T, const IntegrationPoint &ip,
-                  Vector &dd, Vector &uu) override
-   {
-      // dd[0] - density
-      // uu[0] - grad_x
-      // uu[1] - grad_y
-      // uu[2] - grad_z
-      // uu[3] - temperature/scalar field
-
-      real_t di=diff.Eval(T,ip);
-      real_t ll=load.Eval(T,ip);
-      // Computes the physical density using projection.
-      real_t rz=0.5+0.5*std::tanh(beta*(dd[0]-eta)); //projection
-      // Computes the diffusion coefficient at the integration point.
-      real_t fd=di*(std::pow(rz,powerc)+rhomin);
-      // Computes the sum of the energy and the product of the temperature and
-      // the external input at the integration point.
-      real_t rez = 0.5*(uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2])*fd-uu[3]*ll;
-      return rez;
-   }
-
-   /// Returns the derivative of QEnergy with respect to the state vector uu.
-   void QResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                  Vector &dd, Vector &uu, Vector &rr) override
-   {
-      real_t di=diff.Eval(T,ip);
-      real_t ll=load.Eval(T,ip);
-      real_t rz=0.5+0.5*std::tanh(beta*(dd[0]-eta));
-      real_t fd=di*(std::pow(rz,powerc)+rhomin);
-
-      rr[0]=uu[0]*fd;
-      rr[1]=uu[1]*fd;
-      rr[2]=uu[2]*fd;
-      rr[3]=-ll;
-   }
-
-
-   // Returns the derivative, with respect to the density, of the product of
-   // the adjoint field with the residual at the integration point ip.
-   void AQResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                   Vector &dd, Vector &uu, Vector &aa, Vector &rr) override
-   {
-      real_t di=diff.Eval(T,ip);
-      real_t tt=std::tanh(beta*(dd[0]-eta));
-      real_t rz=0.5+0.5*tt;
-      real_t fd=di*powerc*std::pow(rz,powerc-1.0)*0.5*(1.0-tt*tt)*beta;
-
-      rr[0] = -(aa[0]*uu[0]+aa[1]*uu[1]+aa[2]*uu[2])*fd;
-   }
-
-   // Returns the gradient of the residual with respect to the state vector at
-   // the integration point ip.
-   void QGradResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                      Vector &dd, Vector &uu, DenseMatrix &hh) override
-   {
-      real_t di=diff.Eval(T,ip);
-      real_t tt=std::tanh(beta*(dd[0]-eta));
-      real_t rz=0.5+0.5*tt;
-      real_t fd=di*(std::pow(rz,powerc)+rhomin);
-      hh=0.0;
-
-      hh(0,0)=fd;
-      hh(1,1)=fd;
-      hh(2,2)=fd;
-      hh(3,3)=0.0;
-   }
-
-private:
-   mfem::Coefficient& diff; //diffusion coefficient
-   mfem::Coefficient& load; //load coefficient
-   real_t powerc; //penalization coefficient
-   real_t rhomin; //lower bound for the density
-   real_t beta;   //controls the sharpness of the projection
-   real_t eta;    //projection threshold for tanh
-};
-
-/// Provides implementation of an integrator for linear diffusion with
-/// parametrization provided by a density field. The setup is standard for
-/// topology optimization problems.
-class ParametricLinearDiffusion: public ParametricBNLFormIntegrator
-{
-public:
-   ParametricLinearDiffusion(BaseQFunction& qfunm): qfun(qfunm)
-   {
-
-   }
-
-   /// Computes the local energy.
-   real_t GetElementEnergy(const Array<const FiniteElement *> &el,
-                           const Array<const FiniteElement *> &pel,
-                           ElementTransformation &Tr,
-                           const Array<const Vector *> &elfun,
-                           const Array<const Vector *> &pelfun) override;
-
-   /// Computes the element's residual.
-   void AssembleElementVector(const Array<const FiniteElement *> &el,
-                              const Array<const FiniteElement *> &pel,
-                              ElementTransformation &Tr,
-                              const Array<const Vector *> &elfun,
-                              const Array<const Vector *> &pelfun,
-                              const Array<Vector *> &elvec) override;
-
-   /// Computes the stiffness/tangent matrix.
-   void AssembleElementGrad(const Array<const FiniteElement *> &el,
-                            const Array<const FiniteElement *> &pel,
-                            ElementTransformation &Tr,
-                            const Array<const Vector *> &elfun,
-                            const Array<const Vector *> &pelfun,
-                            const Array2D<DenseMatrix *> &elmats) override;
-
-   /// Computes the product of the adjoint solution and the derivative of the
-   /// residual with respect to the parametric fields.
-   void AssemblePrmElementVector(const Array<const FiniteElement *> &el,
-                                 const Array<const FiniteElement *> &pel,
-                                 ElementTransformation &Tr,
-                                 const Array<const Vector *> &elfun,
-                                 const Array<const Vector *> &alfun,
-                                 const Array<const Vector *> &pelfun,
-                                 const Array<Vector *> &elvec) override;
-private:
-   BaseQFunction& qfun;
-};
-
-
-/// Computes an example of nonlinear objective
-/// $\int \rm{field}*\rm{field}*\rm{weight})\rm{d}\Omega_e$.
-class DiffusionObjIntegrator:public BlockNonlinearFormIntegrator
-{
-public:
-
-   DiffusionObjIntegrator()
-   {
-
-   }
-
-   /// Returns the objective contribution at element level.
-   real_t GetElementEnergy(const Array<const FiniteElement *> &el,
-                           ElementTransformation &Tr,
-                           const Array<const Vector *> &elfun) override;
-
-   /// Returns the gradient of the objective contribution at element level.
-   void AssembleElementVector(const Array<const FiniteElement *> &el,
-                              ElementTransformation &Tr,
-                              const Array<const Vector *> &elfun,
-                              const Array<Vector *> &elvec) override;
-};
-
-}
-
-#endif

miniapps/mtop/paramnonlinearform.hpp

@@ -1,300 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_PRMNONLINEARFORM
-#define MFEM_PRMNONLINEARFORM
-
-#include "mfem.hpp"
-
-namespace mfem
-{
-
-/** The abstract base class ParametricBNLFormIntegrator is a generalization of
-    the BlockNonlinearFormIntegrator class suitable for block state and
-    parameter vectors. */
-class ParametricBNLFormIntegrator
-{
-public:
-   /// Compute the local energy
-   virtual real_t GetElementEnergy(const Array<const FiniteElement *>&el,
-                                   const Array<const FiniteElement *>&pel,
-                                   ElementTransformation &Tr,
-                                   const Array<const Vector *>&elfun,
-                                   const Array<const Vector *>&pelfun);
-
-   /// Perform the local action of the BlockNonlinearFormIntegrator
-   virtual void AssembleElementVector(const Array<const FiniteElement *> &el,
-                                      const Array<const FiniteElement *>&pel,
-                                      ElementTransformation &Tr,
-                                      const Array<const Vector *> &elfun,
-                                      const Array<const Vector *>&pelfun,
-                                      const Array<Vector *> &elvec);
-
-   /// Perform the local action of the BlockNonlinearFormIntegrator on element
-   /// faces
-   virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1,
-                                   const Array<const FiniteElement *> &el2,
-                                   const Array<const FiniteElement *> &pel1,
-                                   const Array<const FiniteElement *> &pel2,
-                                   FaceElementTransformations &Tr,
-                                   const Array<const Vector *> &elfun,
-                                   const Array<const Vector *>&pelfun,
-                                   const Array<Vector *> &elvect);
-
-   /// Perform the local action on the parameters of the BNLFormIntegrator
-   virtual void AssemblePrmElementVector(const Array<const FiniteElement *> &el,
-                                         const Array<const FiniteElement *>&pel,
-                                         ElementTransformation &Tr,
-                                         const Array<const Vector *> &elfun,
-                                         const Array<const Vector *> &alfun,
-                                         const Array<const Vector *>&pelfun,
-                                         const Array<Vector *> &pelvec);
-
-   /// Perform the local action on the parameters of the BNLFormIntegrator on
-   /// faces
-   virtual void AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
-                                      const Array<const FiniteElement *> &el2,
-                                      const Array<const FiniteElement *> &pel1,
-                                      const Array<const FiniteElement *> &pel2,
-                                      FaceElementTransformations &Tr,
-                                      const Array<const Vector *> &elfun,
-                                      const Array<const Vector *> &alfun,
-                                      const Array<const Vector *>&pelfun,
-                                      const Array<Vector *> &pelvect);
-
-   /// Assemble the local gradient matrix
-   virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
-                                    const Array<const FiniteElement *>&pel,
-                                    ElementTransformation &Tr,
-                                    const Array<const Vector *> &elfun,
-                                    const Array<const Vector *>&pelfun,
-                                    const Array2D<DenseMatrix *> &elmats);
-
-   /// Assemble the local gradient matrix on faces of the elements
-   virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
-                                 const Array<const FiniteElement *>&el2,
-                                 const Array<const FiniteElement *> &pel1,
-                                 const Array<const FiniteElement *> &pel2,
-                                 FaceElementTransformations &Tr,
-                                 const Array<const Vector *> &elfun,
-                                 const Array<const Vector *>&pelfun,
-                                 const Array2D<DenseMatrix *> &elmats);
-
-
-   virtual ~ParametricBNLFormIntegrator() { }
-};
-
-
-/** @brief A class representing a general parametric block nonlinear operator
-    defined on the Cartesian product of multiple FiniteElementSpace%s. */
-class ParametricBNLForm : public Operator
-{
-protected:
-   /// FE spaces on which the form lives.
-   Array<FiniteElementSpace*> fes;
-
-   /// FE spaces for the parametric fields
-   Array<FiniteElementSpace*> paramfes;
-
-   int paramheight;
-   int paramwidth;
-
-   /// Set of Domain Integrators to be assembled (added).
-   Array<ParametricBNLFormIntegrator*> dnfi;
-
-   /// Set of interior face Integrators to be assembled (added).
-   Array<ParametricBNLFormIntegrator*> fnfi;
-
-   /// Set of Boundary Face Integrators to be assembled (added).
-   Array<ParametricBNLFormIntegrator*> bfnfi;
-   Array<Array<int>*>           bfnfi_marker;
-
-   /** Auxiliary block-vectors for wrapping input and output vectors or holding
-       GridFunction-like block-vector data (e.g. in parallel). */
-   mutable BlockVector xs, ys;
-   mutable BlockVector prmxs, prmys;
-
-   /** Auxiliary block-vectors for holding GridFunction-like block-vector data
-       (e.g. in parallel). */
-   mutable BlockVector xsv;
-
-   /** Auxiliary block-vectors for holding GridFunction-like block-vector data
-       for the parameter fields (e.g. in parallel). */
-   mutable BlockVector xdv;
-   /** Auxiliary block-vectors for holding GridFunction-like block-vector data
-       for the adjoint fields (e.g. in parallel). */
-   mutable BlockVector adv;
-
-   mutable Array2D<SparseMatrix*> Grads, cGrads;
-   mutable BlockOperator *BlockGrad;
-
-   // A list of the offsets
-   Array<int> block_offsets;
-   Array<int> block_trueOffsets;
-   // A list with the offsets for the parametric fields
-   Array<int> paramblock_offsets;
-   Array<int> paramblock_trueOffsets;
-
-   // Array of Arrays of tdofs for each space in 'fes'
-   Array<Array<int> *> ess_tdofs;
-
-   // Array of Arrays of tdofs for each space in 'paramfes'
-   Array<Array<int> *> paramess_tdofs;
-
-   /// Array of pointers to the prolongation matrix of fes, may be NULL
-   Array<const Operator *> P;
-
-   /// Array of pointers to the prolongation matrix of paramfes, may be NULL
-   Array<const Operator *> Pparam;
-
-   /// Array of results of dynamic-casting P to SparseMatrix pointer
-   Array<const SparseMatrix *> cP;
-
-   /// Array of results of dynamic-casting Pparam to SparseMatrix pointer
-   Array<const SparseMatrix *> cPparam;
-
-   /// Indicator if the Operator is part of a parallel run
-   bool is_serial = true;
-
-   /// Indicator if the Operator needs prolongation on assembly
-   bool needs_prolongation = false;
-
-   /// Indicator if the Operator needs prolongation on assembly
-   bool prmneeds_prolongation = false;
-
-   mutable BlockVector aux1, aux2;
-
-   mutable BlockVector prmaux1, prmaux2;
-
-   const BlockVector &Prolongate(const BlockVector &bx) const;
-
-   const BlockVector &ParamProlongate(const BlockVector &bx) const;
-
-   real_t GetEnergyBlocked(const BlockVector &bx, const BlockVector &dx) const;
-
-
-   /// Specialized version of Mult() for BlockVector%s
-   /// Block L-Vector to Block L-Vector
-   void MultBlocked(const BlockVector &bx, const BlockVector &dx,
-                    BlockVector &by) const;
-
-   /// Specialized version of Mult() for BlockVector%s
-   /// Block L-Vector to Block L-Vector
-   /// bx - state vector, ax - adjoint vector, dx - parametric fields
-   /// dy = ax' d(residual(bx))/d(dx)
-   void MultParamBlocked(const BlockVector &bx, const BlockVector & ax,
-                         const BlockVector &dx, BlockVector &dy) const;
-
-
-   /// Specialized version of GetGradient() for BlockVector
-   void ComputeGradientBlocked(const BlockVector &bx, const BlockVector &dx) const;
-
-public:
-   /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
-   ParametricBNLForm();
-
-   /// Construct a BlockNonlinearForm on the given set of FiniteElementSpace%s.
-   ParametricBNLForm(Array<FiniteElementSpace *> &statef,
-                     Array<FiniteElementSpace *> &paramf);
-
-   /// Return the @a k-th FE space of the ParametricBNLForm.
-   FiniteElementSpace *FESpace(int k) { return fes[k]; }
-
-   /// Return the @a k-th parametric FE space of the ParametricBNLForm.
-   FiniteElementSpace *ParamFESpace(int k) { return paramfes[k]; }
-
-
-   /// Return the @a k-th FE space of the BlockNonlinearForm (const version).
-   const FiniteElementSpace *FESpace(int k) const { return fes[k]; }
-
-   /// Return the @a k-th parametric FE space of the BlockNonlinearForm (const
-   /// version).
-   const FiniteElementSpace *ParamFESpace(int k) const { return paramfes[k]; }
-
-   /// Return the integrators
-   Array<ParametricBNLFormIntegrator*>& GetDNFI() { return dnfi;}
-
-
-   /// (Re)initialize the ParametricBNLForm.
-   /** After a call to SetSpaces(), the essential b.c. must be set again. */
-   void SetSpaces(Array<FiniteElementSpace *> &statef,
-                  Array<FiniteElementSpace *> &paramf);
-
-   /// Return the regular dof offsets.
-   const Array<int> &GetBlockOffsets() const { return block_offsets; }
-
-   /// Return the true-dof offsets.
-   const Array<int> &GetBlockTrueOffsets() const { return block_trueOffsets; }
-
-   /// Return the regular dof offsets for the parameters.
-   const Array<int> &ParamGetBlockOffsets() const { return paramblock_offsets; }
-
-   /// Return the true-dof offsets for the parameters.
-   const Array<int> &ParamGetBlockTrueOffsets() const { return paramblock_trueOffsets; }
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(ParametricBNLFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
-
-   /// Adds new Interior Face Integrator.
-   void AddInteriorFaceIntegrator(ParametricBNLFormIntegrator *nlfi)
-   { fnfi.Append(nlfi); }
-
-   /// Adds new Boundary Face Integrator.
-   void AddBdrFaceIntegrator(ParametricBNLFormIntegrator *nlfi)
-   { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
-
-   /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
-       attributes. */
-   void AddBdrFaceIntegrator(ParametricBNLFormIntegrator *nlfi,
-                             Array<int> &bdr_marker);
-
-   /// Set the essential boundary conditions.
-   virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                               Array<Vector *> &rhs);
-
-   /// Set the essential boundary conditions on the parametric fields.
-   virtual void SetParamEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                                    Array<Vector *> &rhs);
-
-
-   /// Computes the energy for a state vector x.
-   virtual real_t GetEnergy(const Vector &x) const;
-
-   /// Method is only called in serial, the parallel version calls MultBlocked
-   /// directly.
-   void Mult(const Vector &x, Vector &y) const override;
-
-   /// Method is only called in serial, the parallel version calls MultBlocked
-   /// directly.
-   virtual void ParamMult(const Vector &x, Vector &y) const;
-
-   /// Method is only called in serial, the parallel version calls
-   /// GetGradientBlocked directly.
-   BlockOperator &GetGradient(const Vector &x) const override;
-
-   /// Set the state fields
-   virtual void SetStateFields(const Vector &xv) const;
-
-   /// Set the adjoint fields
-   virtual void SetAdjointFields(const Vector &av) const;
-
-   /// Set the parameters/design fields
-   virtual void SetParamFields(const Vector &dv) const;
-
-   /// Destructor.
-   virtual ~ParametricBNLForm();
-
-};
-
-}
-
-#endif

miniapps/mtop/parheat.cpp

@@ -1,354 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-//
-//     ----------------------------------------------------------------
-//     ParHeat Miniapp: Gradients of PDE constrained objective function
-//     ----------------------------------------------------------------
-//                           (Parallel Version)
-//
-// The following example computes the gradients of a specified objective
-// function with respect to parametric fields. The objective function is having
-// the following form f(u(\rho)) where u(\rho) is a solution of a specific state
-// problem (in the example that is the diffusion equation), and \rho is a
-// parametric field discretized by finite elements. The parametric field (also
-// called density in topology optimization) controls the coefficients of the
-// state equation. For the considered case, the density controls the diffusion
-// coefficient within the computational domain.
-//
-// For more information, the users are referred to:
-//
-//    Hinze, M.; Pinnau, R.; Ulbrich, M. & Ulbrich, S.
-//    Optimization with PDE Constraints
-//    Springer Netherlands, 2009
-//
-//    Bendsøe, M. P. & Sigmund, O.
-//    Topology Optimization - Theory, Methods and Applications
-//    Springer Verlag, Berlin Heidelberg, 2003
-//
-// Compile with: make parheat
-//
-// Sample runs:
-//
-//    mpirun -np 4 parheat --visualization
-//    mpirun -np 4 parheat --visualization -m ../../data/beam-quad.mesh
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-#include "pparamnonlinearform.hpp"
-#include "mtop_integrators.hpp"
-
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI and HYPRE.
-   mfem::Mpi::Init(argc, argv);
-   int myrank = mfem::Mpi::WorldRank();
-   mfem::Hypre::Init();
-
-   // Parse command-line options.
-   const char *mesh_file = "../../data/star.mesh";
-   int order = 1;
-   bool static_cond = false;
-   int ser_ref_levels = 1;
-   int par_ref_levels = 1;
-   real_t newton_rel_tol = 1e-7;
-   real_t newton_abs_tol = 1e-12;
-   int newton_iter = 10;
-   int print_level = 1;
-   bool visualization = false;
-
-   mfem::OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ser_ref_levels,
-                  "-rs",
-                  "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels,
-                  "-rp",
-                  "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree) or -1 for"
-                  " isoparametric space.");
-   args.AddOption(&visualization,
-                  "-vis",
-                  "--visualization",
-                  "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&newton_rel_tol,
-                  "-rel",
-                  "--relative-tolerance",
-                  "Relative tolerance for the Newton solve.");
-   args.AddOption(&newton_abs_tol,
-                  "-abs",
-                  "--absolute-tolerance",
-                  "Absolute tolerance for the Newton solve.");
-   args.AddOption(&newton_iter,
-                  "-it",
-                  "--newton-iterations",
-                  "Maximum iterations for the Newton solve.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myrank == 0)
-      {
-         args.PrintUsage(std::cout);
-      }
-      return 1;
-   }
-
-   if (myrank == 0)
-   {
-      args.PrintOptions(std::cout);
-   }
-
-   // Read the (serial) mesh from the given mesh file on all processors.  We
-   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   // and volume meshes with the same code.
-   mfem::Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-
-   // Refine the serial mesh on all processors to increase the resolution. In
-   // this example we do 'ref_levels' of uniform refinement. We choose
-   // 'ref_levels' to be the largest number that gives a final mesh with no
-   // more than 10,000 elements.
-   {
-      int ref_levels =
-         (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   // Define a parallel mesh by a partitioning of the serial mesh. Refine
-   // this mesh further in parallel to increase the resolution. Once the
-   // parallel mesh is defined, the serial mesh can be deleted.
-   mfem::ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   mesh.Clear();
-   {
-      for (int l = 0; l < par_ref_levels; l++)
-      {
-         pmesh.UniformRefinement();
-      }
-   }
-
-   // Define the Diffusion coefficient.
-   mfem::ConstantCoefficient* diffco=new mfem::ConstantCoefficient(1.0);
-   // Define the Heat source.
-   mfem::ConstantCoefficient* loadco=new mfem::ConstantCoefficient(1.0);
-   // Define the q-function.
-   mfem::QLinearDiffusion* qfun=new mfem::QLinearDiffusion(*diffco,*loadco,1.0,
-                                                           1e-7,4.0,0.5);
-
-   // Define FE collection and space for the state solution.
-   mfem::H1_FECollection sfec(order, dim);
-   mfem::ParFiniteElementSpace* sfes=new mfem::ParFiniteElementSpace(&pmesh,&sfec,
-                                                                     1);
-   // Define FE collection and space for the density field.
-   mfem::L2_FECollection pfec(order, dim);
-   mfem::ParFiniteElementSpace* pfes=new mfem::ParFiniteElementSpace(&pmesh,&pfec,
-                                                                     1);
-
-   // Define the arrays for the nonlinear form.
-   mfem::Array<mfem::ParFiniteElementSpace*> asfes;
-   mfem::Array<mfem::ParFiniteElementSpace*> apfes;
-
-   asfes.Append(sfes);
-   apfes.Append(pfes);
-
-   // Define parametric block nonlinear form using single scalar H1 field
-   // and L2 scalar density field.
-   mfem::ParParametricBNLForm* nf=new mfem::ParParametricBNLForm(asfes,apfes);
-   // Add a parametric integrator.
-   nf->AddDomainIntegrator(new mfem::ParametricLinearDiffusion(*qfun));
-
-   // Define true block vectors for state, adjoint, resudual.
-   mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets());    solbv=0.0;
-   mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets());    adjbv=0.0;
-   mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets());    resbv=0.0;
-   // Define true block vectors for parametric field and gradients.
-   mfem::BlockVector prmbv; prmbv.Update(nf->ParamGetBlockTrueOffsets());
-   prmbv=0.0;
-   mfem::BlockVector grdbv; grdbv.Update(nf->ParamGetBlockTrueOffsets());
-   grdbv=0.0;
-
-   // Set the BCs for the physics.
-   mfem::Array<mfem::Array<int> *> ess_bdr;
-   mfem::Array<mfem::Vector*>      ess_rhs;
-   ess_bdr.Append(new mfem::Array<int>(pmesh.bdr_attributes.Max()));
-   ess_rhs.Append(nullptr);
-   (*ess_bdr[0]) = 1;
-   nf->SetEssentialBC(ess_bdr,ess_rhs);
-   delete ess_bdr[0];
-
-   // Set the density field to 0.5.
-   prmbv=0.5;
-   // Set the density as parametric field in the parametric BNLForm.
-   nf->SetParamFields(prmbv); //set the density
-
-   // Compute the stiffness/tangent matrix for density prmbv=0.5.
-   mfem::BlockOperator *A = &nf->GetGradient(solbv);
-   mfem::HypreBoomerAMG* prec=new mfem::HypreBoomerAMG();
-   prec->SetPrintLevel(print_level);
-   // Use only block (0,0) as in this case we have a single field.
-   prec->SetOperator(A->GetBlock(0,0));
-
-   // Construct block preconditioner for the BNLForm.
-   mfem::BlockDiagonalPreconditioner *blpr = new mfem::BlockDiagonalPreconditioner(
-      nf->GetBlockTrueOffsets());
-   blpr->SetDiagonalBlock(0,prec);
-
-   // Define the solvers.
-   mfem::GMRESSolver *gmres;
-   gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
-   gmres->SetAbsTol(newton_abs_tol/10);
-   gmres->SetRelTol(newton_rel_tol/10);
-   gmres->SetMaxIter(100);
-   gmres->SetPrintLevel(print_level);
-   gmres->SetPreconditioner(*blpr);
-   gmres->SetOperator(*A);
-
-
-   // Solve the problem.
-   solbv=0.0;
-   nf->Mult(solbv,resbv); resbv.Neg(); //compute RHS
-   gmres->Mult(resbv, solbv);
-
-   // Compute the energy of the state system.
-   real_t energy = nf->GetEnergy(solbv);
-   if (myrank==0)
-   {
-      std::cout << "energy =" << energy << std::endl;
-   }
-
-   // Define the block nonlinear form utilized for representing the objective -
-   // use the state array from the BNLForm.
-   mfem::ParBlockNonlinearForm* ob=new mfem::ParBlockNonlinearForm(asfes);
-   // Add the integrator for the objective.
-   ob->AddDomainIntegrator(new mfem::DiffusionObjIntegrator());
-
-   // Compute the objective.
-   real_t obj=ob->GetEnergy(solbv);
-   if (myrank==0)
-   {
-      std::cout << "Objective =" << obj << std::endl;
-   }
-
-   // Solve the adjoint.
-   {
-      mfem::BlockVector adjrhs; adjrhs.Update(nf->GetBlockTrueOffsets());  adjrhs=0.0;
-      // Compute the RHS for the adjoint, i.e., the gradients with respect to
-      // the parametric fields.
-      ob->Mult(solbv, adjrhs);
-      // Get the tangent matrix from the state problem. We do not need to
-      // transpose the operator for diffusion.  Compute the adjoint solution.
-      gmres->Mult(adjrhs, adjbv);
-   }
-
-   // Compute gradients.
-   // First set the adjoint field.
-   nf->SetAdjointFields(adjbv);
-   // Set the state field.
-   nf->SetStateFields(solbv);
-   // Call the parametric Mult.
-   nf->ParamMult(prmbv, grdbv);
-
-   // Dump out the data.
-   if (visualization)
-   {
-      mfem::ParaViewDataCollection *dacol=new mfem::ParaViewDataCollection("ParHeat",
-                                                                           &pmesh);
-      mfem::ParGridFunction gfgrd(pfes); gfgrd.SetFromTrueDofs(grdbv.GetBlock(0));
-      mfem::ParGridFunction gfdns(pfes); gfdns.SetFromTrueDofs(prmbv.GetBlock(0));
-      // Define state grid function.
-      mfem::ParGridFunction gfsol(sfes); gfsol.SetFromTrueDofs(solbv.GetBlock(0));
-      mfem::ParGridFunction gfadj(sfes); gfadj.SetFromTrueDofs(adjbv.GetBlock(0));
-
-      dacol->SetLevelsOfDetail(order);
-      dacol->RegisterField("sol", &gfsol);
-      dacol->RegisterField("adj", &gfadj);
-      dacol->RegisterField("dns", &gfdns);
-      dacol->RegisterField("grd", &gfgrd);
-
-      dacol->SetTime(1.0);
-      dacol->SetCycle(1);
-      dacol->Save();
-
-      delete dacol;
-   }
-
-   // FD check
-   {
-      mfem::BlockVector prtbv;
-      mfem::BlockVector tmpbv;
-      prtbv.Update(nf->ParamGetBlockTrueOffsets());
-      tmpbv.Update(nf->ParamGetBlockTrueOffsets());
-      prtbv.GetBlock(0).Randomize();
-      prtbv*=1.0;
-      real_t lsc=1.0;
-
-      real_t gQoI=ob->GetEnergy(solbv);
-      real_t lQoI;
-
-      real_t nd=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,prtbv);
-      real_t td=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,grdbv);
-      td=td/nd;
-
-      for (int l = 0; l < 10; l++)
-      {
-         lsc/=10.0;
-         prtbv/=10.0;
-         add(prmbv,prtbv,tmpbv);
-         nf->SetParamFields(tmpbv);
-         // Solve the physics.
-         solbv=0.0;
-         nf->Mult(solbv,resbv); resbv.Neg(); //compute RHS
-         A = &nf->GetGradient(solbv);
-         prec->SetPrintLevel(0);
-         prec->SetOperator(A->GetBlock(0,0));
-         gmres->SetOperator(*A);
-         gmres->SetPrintLevel(0);
-         gmres->Mult(resbv,solbv);
-         // Compute the objective.
-         lQoI=ob->GetEnergy(solbv);
-         real_t ld=(lQoI-gQoI)/lsc;
-         if (myrank==0)
-         {
-            std::cout << "dx=" << lsc <<" FD approximation=" << ld/nd
-                      << " adjoint gradient=" << td
-                      << " err=" << std::fabs(ld/nd-td) << std::endl;
-         }
-      }
-   }
-
-   delete ob;
-   delete gmres;
-   delete blpr;
-   delete prec;
-
-   delete nf;
-   delete pfes;
-   delete sfes;
-
-   delete qfun;
-   delete loadco;
-   delete diffco;
-
-   return 0;
-}

miniapps/mtop/pparamnonlinearform.cpp

@@ -1,362 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "mfem.hpp"
-#include "pparamnonlinearform.hpp"
-
-#ifdef MFEM_USE_MPI
-
-namespace mfem
-{
-
-ParParametricBNLForm::ParParametricBNLForm(Array<ParFiniteElementSpace *>
-                                           &statef,
-                                           Array<ParFiniteElementSpace *> &paramf)
-   :ParametricBNLForm()
-{
-   pBlockGrad = nullptr;
-   SetParSpaces(statef,paramf);
-}
-
-void ParParametricBNLForm::SetParSpaces(Array<ParFiniteElementSpace *> &statef,
-                                        Array<ParFiniteElementSpace *> &paramf)
-{
-   delete pBlockGrad;
-   pBlockGrad = nullptr;
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete phBlockGrad(s1,s2);
-      }
-   }
-
-   Array<FiniteElementSpace *> serialSpaces(statef.Size());
-   Array<FiniteElementSpace *> prmserialSpaces(paramf.Size());
-   for (int s=0; s<statef.Size(); s++)
-   {
-      serialSpaces[s] = (FiniteElementSpace *) statef[s];
-   }
-   for (int s=0; s<paramf.Size(); s++)
-   {
-      prmserialSpaces[s] = (FiniteElementSpace *) paramf[s];
-   }
-
-   SetSpaces(serialSpaces,prmserialSpaces);
-
-   phBlockGrad.SetSize(fes.Size(), fes.Size());
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2) = new OperatorHandle(Operator::Hypre_ParCSR);
-      }
-   }
-}
-
-ParFiniteElementSpace * ParParametricBNLForm::ParFESpace(int k)
-{
-   return (ParFiniteElementSpace *)fes[k];
-}
-
-const ParFiniteElementSpace *ParParametricBNLForm::ParFESpace(int k) const
-{
-   return (const ParFiniteElementSpace *)fes[k];
-}
-
-
-ParFiniteElementSpace * ParParametricBNLForm::ParParamFESpace(int k)
-{
-   return (ParFiniteElementSpace *)paramfes[k];
-}
-
-const ParFiniteElementSpace *ParParametricBNLForm::ParParamFESpace(int k) const
-{
-   return (const ParFiniteElementSpace *)paramfes[k];
-}
-
-// Here, rhs is a true dof vector
-void ParParametricBNLForm::SetEssentialBC(const
-                                          Array<Array<int> *>&bdr_attr_is_ess,
-                                          Array<Vector *> &rhs)
-{
-   Array<Vector *> nullarray(fes.Size());
-   nullarray = NULL;
-
-   ParametricBNLForm::SetEssentialBC(bdr_attr_is_ess, nullarray);
-
-   for (int s = 0; s < fes.Size(); ++s)
-   {
-      if (rhs[s])
-      {
-         rhs[s]->SetSubVector(*ess_tdofs[s], 0.0);
-      }
-   }
-}
-
-void ParParametricBNLForm::SetParamEssentialBC(const
-                                               Array<Array<int> *>&bdr_attr_is_ess,
-                                               Array<Vector *> &rhs)
-{
-   Array<Vector *> nullarray(fes.Size());
-   nullarray = NULL;
-
-   ParametricBNLForm::SetParamEssentialBC(bdr_attr_is_ess, nullarray);
-
-   for (int s = 0; s < paramfes.Size(); ++s)
-   {
-      if (rhs[s])
-      {
-         rhs[s]->SetSubVector(*paramess_tdofs[s], 0.0);
-      }
-   }
-}
-
-real_t ParParametricBNLForm::GetEnergy(const Vector &x) const
-{
-   xs_true.Update(const_cast<Vector&>(x), block_trueOffsets);
-   xs.Update(block_offsets);
-
-   for (int s = 0; s < fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   real_t enloc = ParametricBNLForm::GetEnergyBlocked(xs,xdv);
-   real_t englo = 0.0;
-
-   MPI_Allreduce(&enloc, &englo, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
-                 ParFESpace(0)->GetComm());
-
-   return englo;
-}
-
-void ParParametricBNLForm::Mult(const Vector &x, Vector &y) const
-{
-   xs_true.Update(const_cast<Vector&>(x), block_trueOffsets);
-   ys_true.Update(y, block_trueOffsets);
-   xs.Update(block_offsets);
-   ys.Update(block_offsets);
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   ParametricBNLForm::MultBlocked(xs, xdv, ys);
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->MultTranspose(
-         ys.GetBlock(s), ys_true.GetBlock(s));
-
-      ys_true.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
-   }
-}
-
-/// Block T-Vector to Block T-Vector
-void ParParametricBNLForm::ParamMult(const Vector &x, Vector &y) const
-{
-   xs_true.Update(const_cast<Vector&>(x), paramblock_trueOffsets);
-   ys_true.Update(y, paramblock_trueOffsets);
-   prmxs.Update(paramblock_offsets);
-   prmys.Update(paramblock_offsets);
-
-   for (int s=0; s<paramfes.Size(); ++s)
-   {
-      paramfes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), prmxs.GetBlock(s));
-   }
-
-   ParametricBNLForm::MultParamBlocked(xsv,adv,xdv,prmys);
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s=0; s<paramfes.Size(); ++s)
-   {
-      paramfes[s]->GetProlongationMatrix()->MultTranspose(
-         prmys.GetBlock(s), ys_true.GetBlock(s));
-
-      ys_true.GetBlock(s).SetSubVector(*paramess_tdofs[s], 0.0);
-   }
-
-}
-
-/// Return the local gradient matrix for the given true-dof vector x
-const BlockOperator & ParParametricBNLForm::GetLocalGradient(
-   const Vector &x) const
-{
-   xs_true.Update(const_cast<Vector&>(x), block_trueOffsets);
-   xs.Update(block_offsets);
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   ParametricBNLForm::ComputeGradientBlocked(xs,
-                                             xdv); // (re)assemble Grad with b.c.
-
-   delete BlockGrad;
-   BlockGrad = new BlockOperator(block_offsets);
-
-   for (int i = 0; i < fes.Size(); ++i)
-   {
-      for (int j = 0; j < fes.Size(); ++j)
-      {
-         BlockGrad->SetBlock(i, j, Grads(i, j));
-      }
-   }
-   return *BlockGrad;
-}
-
-// Set the operator type id for the parallel gradient matrix/operator.
-void ParParametricBNLForm::SetGradientType(Operator::Type tid)
-{
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2)->SetType(tid);
-      }
-   }
-}
-
-BlockOperator & ParParametricBNLForm::GetGradient(const Vector &x) const
-{
-   if (pBlockGrad == NULL)
-   {
-      pBlockGrad = new BlockOperator(block_trueOffsets);
-   }
-
-   Array<const ParFiniteElementSpace *> pfes(fes.Size());
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      pfes[s1] = ParFESpace(s1);
-
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2)->Clear();
-      }
-   }
-
-   GetLocalGradient(x); // gradients are stored in 'Grads'
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                        Ph(phBlockGrad(s1,s2)->Type()),
-                        Rh(phBlockGrad(s1,s2)->Type());
-
-         if (s1 == s2)
-         {
-            dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
-                                   pfes[s1]->GetDofOffsets(), Grads(s1,s1));
-            Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-            OperatorHandle Ae;
-            Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
-         }
-         else
-         {
-            dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
-                                        pfes[s1]->GlobalVSize(),
-                                        pfes[s2]->GlobalVSize(),
-                                        pfes[s1]->GetDofOffsets(),
-                                        pfes[s2]->GetDofOffsets(),
-                                        Grads(s1,s2));
-            Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-            phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-            phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-            phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
-         }
-
-         pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
-      }
-   }
-
-   return *pBlockGrad;
-}
-
-ParParametricBNLForm::~ParParametricBNLForm()
-{
-   delete pBlockGrad;
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete phBlockGrad(s1,s2);
-      }
-   }
-}
-
-
-void ParParametricBNLForm::SetStateFields(const Vector &xv) const
-{
-   xs_true.Update(const_cast<Vector&>(xv), block_trueOffsets);
-   xsv.Update(block_offsets);
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xsv.GetBlock(s));
-   }
-}
-
-
-void ParParametricBNLForm::SetAdjointFields(const Vector &av) const
-{
-   xs_true.Update(const_cast<Vector&>(av), block_trueOffsets);
-   adv.Update(block_offsets);
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), adv.GetBlock(s));
-   }
-}
-
-void ParParametricBNLForm::SetParamFields(const Vector &dv) const
-{
-   xs_true.Update(const_cast<Vector&>(dv),paramblock_trueOffsets);
-   xdv.Update(paramblock_offsets);
-   for (int s=0; s<paramfes.Size(); ++s)
-   {
-      paramfes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xdv.GetBlock(s));
-   }
-}
-
-}
-
-#endif

miniapps/mtop/pparamnonlinearform.hpp

@@ -1,114 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_PPRMNONLINEARFORM
-#define MFEM_PPRMNONLINEARFORM
-
-
-#ifdef MFEM_USE_MPI
-
-#include "mfem.hpp"
-#include "paramnonlinearform.hpp"
-
-namespace mfem
-{
-
-/** @brief A class representing a general parametric parallel block nonlinear
-    operator defined on the Cartesian product of multiple
-    ParFiniteElementSpace%s. */
-/** The ParParametricBNLForm takes as input, and returns as output, vectors on
-    the true dofs. */
-class ParParametricBNLForm : public ParametricBNLForm
-{
-protected:
-   mutable BlockVector xs_true, ys_true;
-   mutable Array2D<OperatorHandle *> phBlockGrad;
-   mutable BlockOperator *pBlockGrad;
-
-public:
-   /// Computes the energy of the system
-   real_t GetEnergy(const Vector &x) const override;
-
-   /// Construct an empty ParParametricBNLForm. Initialize with SetParSpaces().
-   ParParametricBNLForm() : pBlockGrad(nullptr) { }
-
-   /** @brief Construct a ParParametricBNLForm on the given set of
-       parametric and state ParFiniteElementSpace%s. */
-   ParParametricBNLForm(Array<ParFiniteElementSpace *> &statef,
-                        Array<ParFiniteElementSpace *> &paramf);
-
-   /// Return the @a k-th parallel FE state space of the ParParametricBNLForm.
-   ParFiniteElementSpace *ParFESpace(int k);
-   /** @brief Return the @a k-th parallel FE state space of the
-       ParParametricBNLForm (const version). */
-   const ParFiniteElementSpace *ParFESpace(int k) const;
-
-   /// Return the @a k-th parallel FE parameters space of the
-   /// ParParametricBNLForm.
-   ParFiniteElementSpace *ParParamFESpace(int k);
-   /** @brief Return the @a k-th parallel FE parameters space of the
-       ParParametricBNLForm (const version). */
-   const ParFiniteElementSpace *ParParamFESpace(int k) const;
-
-   /** @brief Set the parallel FE spaces for the state and the parametric
-    *  fields. After a call to SetParSpaces(), the essential b.c. and the
-    *  gradient-type (if different from the default) must be set again. */
-   void SetParSpaces(Array<ParFiniteElementSpace *> &statef,
-                     Array<ParFiniteElementSpace *> &paramf);
-
-   /// Set the state essential BCs. Here, rhs is a true dof vector!
-   void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                       Array<Vector *> &rhs) override;
-
-   // Set the essential BCs for the parametric fields. Here, rhs is a true dof
-   // vector!
-   void SetParamEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                            Array<Vector *> &rhs) override;
-
-
-   /** @brief Calculates the residual for a state input given by block T-Vector.
-    *  The result is Block T-Vector! The parametric fields should be set in
-    *  advance by calling SetParamFields(). */
-   void Mult(const Vector &x, Vector &y) const override;
-
-   /** @brief Calculates the product of the adjoint field and the derivative of
-    *  the state residual with respect to the parametric fields. The adjoint and
-    *  the state fields should be set in advance by calling SetAdjointFields()
-    *  and SetStateFields(). The input and the result are block T-Vectors!*/
-   void ParamMult(const Vector &x, Vector &y) const override;
-
-   /// Return the local block gradient matrix for the given true-dof vector x
-   const BlockOperator &GetLocalGradient(const Vector &x) const;
-
-   /// Return the block gradient matrix for the given true-dof vector x
-   BlockOperator &GetGradient(const Vector &x) const override;
-
-   /** @brief Set the operator type id for the blocks of the parallel gradient
-       matrix/operator. The default type is Operator::Hypre_ParCSR. */
-   void SetGradientType(Operator::Type tid);
-
-   /// Destructor.
-   virtual ~ParParametricBNLForm();
-
-   /// Set the state fields
-   void SetStateFields(const Vector &xv) const override;
-
-   /// Set the adjoint fields
-   void SetAdjointFields(const Vector &av) const override;
-
-   /// Set the parameters/design fields
-   void SetParamFields(const Vector &dv) const override;
-};
-
-}
-
-#endif
-#endif

miniapps/mtop/seqheat.cpp

@@ -1,308 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-//
-//     ----------------------------------------------------------------
-//     SeqHeat Miniapp: Gradients of PDE constrained objective function
-//     ----------------------------------------------------------------
-//                          (Sequential Version)
-//
-// The following example computes the gradients of a specified objective
-// function with respect to parametric fields. The objective function is having
-// the following form f(u(\rho)) where u(\rho) is a solution of a specific state
-// problem (in the example that is the diffusion equation), and \rho is a
-// parametric field discretized by finite elements. The parametric field (also
-// called density in topology optimization) controls the coefficients of the
-// state equation. For the considered case, the density controls the diffusion
-// coefficient within the computational domain.
-//
-// For more information, the users are referred to:
-//
-//    Hinze, M.; Pinnau, R.; Ulbrich, M. & Ulbrich, S.
-//    Optimization with PDE Constraints
-//    Springer Netherlands, 2009
-//
-//    Bendsøe, M. P. & Sigmund, O.
-//    Topology Optimization - Theory, Methods and Applications
-//    Springer Verlag, Berlin Heidelberg, 2003
-//
-// Compile with: make seqheat
-//
-// Sample runs:
-//
-//    seqheat -m ../../data/star-mixed.mesh
-//    seqheat --visualization
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-#include "mtop_integrators.hpp"
-
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   const char *mesh_file = "../../data/star.vtk";
-   int ser_ref_levels = 1;
-   int order = 2;
-   bool visualization = false;
-   real_t newton_rel_tol = 1e-4;
-   real_t newton_abs_tol = 1e-6;
-   int newton_iter = 10;
-   int print_level = 0;
-
-   mfem::OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
-   args.AddOption(&ser_ref_levels,
-                  "-rs",
-                  "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&order,
-                  "-o",
-                  "--order",
-                  "Order (degree) of the finite elements.");
-   args.AddOption(&visualization,
-                  "-vis",
-                  "--visualization",
-                  "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&newton_rel_tol,
-                  "-rel",
-                  "--relative-tolerance",
-                  "Relative tolerance for the Newton solve.");
-   args.AddOption(&newton_abs_tol,
-                  "-abs",
-                  "--absolute-tolerance",
-                  "Absolute tolerance for the Newton solve.");
-   args.AddOption(&newton_iter,
-                  "-it",
-                  "--newton-iterations",
-                  "Maximum iterations for the Newton solve.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(std::cout);
-      return 1;
-   }
-   args.PrintOptions(std::cout);
-
-   // Read the (serial) mesh from the given mesh file on all processors. We
-   // can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
-   // with the same code.
-   mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
-
-   // Refine the mesh in serial to increase the resolution. In this example
-   // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
-   // a command-line parameter.
-   for (int lev = 0; lev < ser_ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // Diffusion coefficient
-   mfem::ConstantCoefficient* diffco=new mfem::ConstantCoefficient(1.0);
-   // Heat source
-   mfem::ConstantCoefficient* loadco=new mfem::ConstantCoefficient(1.0);
-   // Define the q-function
-   mfem::QLinearDiffusion* qfun=new mfem::QLinearDiffusion(*diffco,*loadco,1.0,
-                                                           1e-7,4.0,0.5);
-
-   // Define FE collection and space for the state solution
-   mfem::H1_FECollection sfec(order, dim);
-   mfem::FiniteElementSpace* sfes=new mfem::FiniteElementSpace(mesh,&sfec,1);
-   // Define FE collection and space for the density field
-   mfem::L2_FECollection pfec(order, dim);
-   mfem::FiniteElementSpace* pfes=new mfem::FiniteElementSpace(mesh,&pfec,1);
-
-   // Define the arrays for the nonlinear form
-   mfem::Array<mfem::FiniteElementSpace*> asfes;
-   mfem::Array<mfem::FiniteElementSpace*> apfes;
-
-   asfes.Append(sfes);
-   apfes.Append(pfes);
-   // Define parametric block nonlinear form using single scalar H1 field
-   // and L2 scalar density field
-   mfem::ParametricBNLForm* nf=new mfem::ParametricBNLForm(asfes,apfes);
-   // Add the parametric integrator
-   nf->AddDomainIntegrator(new mfem::ParametricLinearDiffusion(*qfun));
-
-   // Define true block vectors for state, adjoint, residual
-   mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets());    solbv=0.0;
-   mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets());    adjbv=0.0;
-   mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets());    resbv=0.0;
-   // Define true block vectors for parametric field and gradients
-   mfem::BlockVector prmbv; prmbv.Update(nf->ParamGetBlockTrueOffsets());
-   prmbv=0.0;
-   mfem::BlockVector grdbv; grdbv.Update(nf->ParamGetBlockTrueOffsets());
-   grdbv=0.0;
-
-   // Set the BC for the physics
-   mfem::Array<mfem::Array<int> *> ess_bdr;
-   mfem::Array<mfem::Vector*>      ess_rhs;
-   ess_bdr.Append(new mfem::Array<int>(mesh->bdr_attributes.Max()));
-   ess_rhs.Append(nullptr);
-   (*ess_bdr[0]) = 1;
-   nf->SetEssentialBC(ess_bdr,ess_rhs);
-   delete ess_bdr[0];
-
-   // Define the linear solvers
-   mfem::GMRESSolver *gmres;
-   gmres = new mfem::GMRESSolver();
-   gmres->SetAbsTol(newton_abs_tol/10);
-   gmres->SetRelTol(newton_rel_tol/10);
-   gmres->SetMaxIter(300);
-   gmres->SetPrintLevel(print_level);
-
-   // Define the Newton solver
-   mfem::NewtonSolver *ns;
-   ns = new mfem::NewtonSolver();
-   ns->iterative_mode = true;
-   ns->SetSolver(*gmres);
-   ns->SetOperator(*nf);
-   ns->SetPrintLevel(print_level);
-   ns->SetRelTol(newton_rel_tol);
-   ns->SetAbsTol(newton_abs_tol);
-   ns->SetMaxIter(newton_iter);
-
-   // Solve the problem
-   // Set the density to 0.5
-   prmbv=0.5;
-   nf->SetParamFields(prmbv); // Set the density
-   // Define the RHS
-   mfem::Vector b;
-   solbv=0.0;
-   // Newton solve
-   ns->Mult(b, solbv);
-
-   // Compute the residual
-   nf->Mult(solbv,resbv);
-   std::cout<<"Norm residual="<<resbv.Norml2()<<std::endl;
-
-   // Compute the energy of the state system
-   real_t energy = nf->GetEnergy(solbv);
-   std::cout<<"energy ="<< energy<<std::endl;
-
-   // Define the block nonlinear form utilized for representing the
-   // objective. The input is the state array asfes defined earlier.
-   mfem::BlockNonlinearForm* ob=new mfem::BlockNonlinearForm(asfes);
-
-   // Add the integrator for the objective
-   ob->AddDomainIntegrator(new mfem::DiffusionObjIntegrator());
-
-   // Compute the objective
-   real_t obj=ob->GetEnergy(solbv);
-   std::cout<<"Objective ="<<obj<<std::endl;
-
-   // Solve the adjoint
-   {
-      mfem::BlockVector adjrhs; adjrhs.Update(nf->GetBlockTrueOffsets());  adjrhs=0.0;
-      // Compute the RHS for the adjoint
-      ob->Mult(solbv, adjrhs);
-      // Get the tangent matrix from the state problem
-      mfem::BlockOperator& A=nf->GetGradient(solbv);
-      // We do not need to transpose the operator for diffusion
-      gmres->SetOperator(A.GetBlock(0,0));
-      // Compute the adjoint solution
-      gmres->Mult(adjrhs.GetBlock(0), adjbv.GetBlock(0));
-   }
-
-   // Compute gradients
-   nf->SetAdjointFields(adjbv);
-   nf->SetStateFields(solbv);
-   nf->ParamMult(prmbv, grdbv);
-
-   // Dump out the data
-   if (visualization)
-   {
-      mfem::ParaViewDataCollection *dacol=new mfem::ParaViewDataCollection("SeqHeat",
-                                                                           mesh);
-      mfem::GridFunction gfgrd(pfes); gfgrd.SetFromTrueDofs(grdbv.GetBlock(0));
-      mfem::GridFunction gfdns(pfes); gfdns.SetFromTrueDofs(prmbv.GetBlock(0));
-      // Define state grid function
-      mfem::GridFunction gfsol(sfes); gfsol.SetFromTrueDofs(solbv.GetBlock(0));
-      mfem::GridFunction gfadj(sfes); gfadj.SetFromTrueDofs(adjbv.GetBlock(0));
-
-      dacol->SetLevelsOfDetail(order);
-      dacol->RegisterField("sol", &gfsol);
-      dacol->RegisterField("adj", &gfadj);
-      dacol->RegisterField("dns", &gfdns);
-      dacol->RegisterField("grd", &gfgrd);
-
-      dacol->SetTime(1.0);
-      dacol->SetCycle(1);
-      dacol->Save();
-
-      delete dacol;
-   }
-
-   // FD check
-   {
-      // Perturbation vector
-      mfem::BlockVector prtbv;
-      mfem::BlockVector tmpbv;
-      prtbv.Update(nf->ParamGetBlockTrueOffsets());
-      tmpbv.Update(nf->ParamGetBlockTrueOffsets());
-      // Generate the perturbation
-      prtbv.GetBlock(0).Randomize();
-      prtbv*=1.0;
-      // Scaling parameter
-      real_t lsc=1.0;
-
-      // Compute initial objective
-      real_t gQoI=ob->GetEnergy(solbv);
-      real_t lQoI;
-
-      // Norm of the perturbation
-      real_t nd=mfem::InnerProduct(prtbv,prtbv);
-      // Projection of the adjoint gradient on the perturbation
-      real_t td=mfem::InnerProduct(prtbv,grdbv);
-      // Normalize the directional derivative
-      td=td/nd;
-
-      for (int l = 0; l < 10; l++)
-      {
-         lsc/=10.0;
-         // Scale the perturbation
-         prtbv/=10.0;
-         // Add the perturbation to the original density
-         add(prmbv,prtbv,tmpbv);
-         nf->SetParamFields(tmpbv);
-         // Solve the physics
-         ns->Mult(b,solbv);
-         // Compute the objective
-         lQoI=ob->GetEnergy(solbv);
-         // FD approximation
-         real_t ld=(lQoI-gQoI)/lsc;
-         std::cout << "dx=" << lsc << " FD gradient=" << ld/nd
-                   << " adjoint gradient=" << td
-                   << " err=" << std::fabs(ld/nd-td) << std::endl;
-      }
-   }
-
-   delete ob;
-
-   delete ns;
-   delete gmres;
-
-   delete nf;
-   delete pfes;
-   delete sfes;
-
-   delete qfun;
-   delete loadco;
-   delete diffco;
-
-   delete mesh;
-
-   return 0;
-}

miniapps/mtop/square-extended.mesh

@@ -0,0 +1,35 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+1
+1 3 0 1 2 3
+
+boundary
+4
+1 1 0 1
+2 1 1 2
+3 1 2 3
+4 1 3 0
+
+vertices
+4
+2
+-1 -1
+1 -1
+1 1
+-1 1