An example of wave equation with time-dependent inhomogeneous Dirichlet boundary condition #4182

Wayne901 · 2024-03-09T17:01:20Z

Hi everyone, since there's no example for problems with time-dependent (especially nonzero) Dirichlet boundary conditions, I try to write one based on ex23. A special nature of ex23 is that wave equation has 2nd-order time derivatives. The code is as follow and I cannot guarantee it's correctness yet. Hope this can offer some help to people with similar needs.

ex23 applys a zero dirichlet boundary condtion, which doesn't need eliminating the rhs at all. Here I try a testcase here. In the $[0,1]^2$ region, the boundary conditions are
$u(x,t)=sin(4\pi t), \frac{du}{dt}(x,t)=4\pi cos(4\pi t), \frac{d^2u}{dt^2}(x,t)=-16\pi^2 sin(4\pi t)$, when $x=0, 1/6 < y < 5/6$, and 0 otherwise. In this case I need to set dirichlet BC for both $u, \frac{du}{dt}, \frac{d^2u}{dt^2}$.

For Newmark method, the wave equation is discretized as $[M+\beta(\Delta t)^2K]\frac{d^2u(t+\Delta t)}{dt^2} = -K[u(t)+\Delta t\frac{du(t)}{dt}+(\frac{1}{2}-\beta)(\Delta t)^2\frac{d^2u(t)}{dt^2}]$. Following @v-dobrev's suggestion #1720 , I derive the solution d2udt2 at t+dt after elimination as
Snipaste_2024-03-10_01-45-49
The u or dudt in photo represents u(t), du(t)/dt. I and B represents DoFs on interior region and boundary, respectively.

I have some quesitons about improving it

Is the ImplicitSolve() correct about eliminating the Dirichlet BC contributions?
I cannot access dt and beta inside ImplicitSolve, but I need dt and 0.5-beta to eliminate rhs. The current code calculates dt first, and then compute beta based on the input variable fac0=beta*dt*dt. Clearly it's only limit to the Newmark ODESolver. How to improve it for more general solvers?
Is there a better method to set boudnary value for solution vector, when the boundary value dependent on coordinate? Current method seems not robust, only limit to serial codes:

u[ess_tdof_list_main[dofi]] = u_gf[ess_tdof_list_main[dofi]];
dudt[ess_tdof_list_main[dofi]] = dudt_gf[ess_tdof_list_main[dofi]];

The whole code.

#include "mfem.hpp"
#include <fstream>
#include <iostream>

using namespace std;
using namespace mfem;

double InitialSolution(const Vector &x)
{
   return 0.0;
}


double InitialRate(const Vector &x)
{
   return 0.0;
}

double BoundaryValue_u(const Vector &x, double time)
{
    if (time <= 0.5 && x(0) <= 1e-7 && x(1) < 5.0/6.0 && x(1) > 1.0/6.0)
      return sin(4*M_PI*time);
    else
        return 0.;
}

double BoundaryValue_dudt(const Vector &x, double time)
{
    if (time <= 0.5 && x(0) <= 1e-7 && x(1) < 5.0/6.0 && x(1) > 1.0/6.0)
      return 4 * M_PI * cos(4*M_PI*time);
    else
        return 0.;
}

double BoundaryValue_d2udt2(const Vector &x, double time)
{
    if (time <= 0.5 && x(0) <= 1e-7 && x(1) < 5.0/6.0 && x(1) > 1.0/6.0)
      return -16 * M_PI * M_PI * sin(4*M_PI*time);
    else
        return 0.;
}

/** After spatial discretization, the conduction model can be written as:
 *
 *     d^2u/dt^2 = M^{-1}(-Ku)
 *
 *  where u is the vector representing the temperature, M is the mass matrix,
 *  and K is the diffusion operator with diffusivity depending on u:
 *  (\kappa + \alpha u).
 *
 *  Class WaveOperator represents the right-hand side of the above ODE.
 */
class WaveOperator : public SecondOrderTimeDependentOperator
{
protected:
   FiniteElementSpace &fespace;
   Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
   Array<int> ess_bdr;

   BilinearForm *M;
   BilinearForm *K;

   SparseMatrix Mmat, Kmat, Kmat0;
   SparseMatrix *T; // T = M + dt K
   double previous_time;
   double current_time;
   double current_dt;

   UMFPackSolver M_solver;   

   CGSolver T_solver; // Implicit solver for T = M + fac0*K
   DSmoother T_prec;  // Preconditioner for the implicit solver

   Coefficient *c2;

   FunctionCoefficient *u_bdr;
   FunctionCoefficient *dudt_bdr;
   FunctionCoefficient *d2udt2_bdr;

   mutable Vector z; // auxiliary vector

public:
   WaveOperator(FiniteElementSpace &f, Array<int> &ess_bdr,
                double speed, double init_time);

   using SecondOrderTimeDependentOperator::Mult;
   virtual void Mult(const Vector &u, const Vector &du_dt,
                     Vector &d2udt2) const;

   /** Solve the Backward-Euler equation:
       d2udt2 = f(u + fac0*d2udt2,dudt + fac1*d2udt2, t),
       for the unknown d2udt2. */
   using SecondOrderTimeDependentOperator::ImplicitSolve;
   virtual void ImplicitSolve(const double fac0, const double fac1,
                              const Vector &u, const Vector &dudt, Vector &d2udt2);

   ///
   void SetParameters(const Vector &u);

   virtual ~WaveOperator();
};


WaveOperator::WaveOperator(FiniteElementSpace &f,
                           Array<int> &ess_bdr, double speed, double init_time)
   : SecondOrderTimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL),
     K(NULL),
     T(NULL), current_dt(1.0e-2), z(height), previous_time(init_time), current_time(0.)
{
   const double rel_tol = 1e-8;

   this->ess_bdr = ess_bdr;
   fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);

   c2 = new ConstantCoefficient(speed*speed);

   K = new BilinearForm(&fespace);
   K->AddDomainIntegrator(new DiffusionIntegrator(*c2));
   K->Assemble();

   Array<int> dummy;
   K->FormSystemMatrix(dummy, Kmat0);
   K->FormSystemMatrix(ess_tdof_list, Kmat);

   M = new BilinearForm(&fespace);
   M->AddDomainIntegrator(new MassIntegrator());
   M->Assemble();
   M->FormSystemMatrix(ess_tdof_list, Mmat);

   M_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
   M_solver.SetOperator(Mmat);

   T_solver.iterative_mode = false;
   T_solver.SetRelTol(rel_tol);
   T_solver.SetAbsTol(0.0);
   T_solver.SetMaxIter(100);
   T_solver.SetPrintLevel(0);
   T_solver.SetPreconditioner(T_prec);
   T = NULL;

   u_bdr = new FunctionCoefficient(BoundaryValue_u);
   dudt_bdr = new FunctionCoefficient(BoundaryValue_dudt);
   d2udt2_bdr = new FunctionCoefficient(BoundaryValue_d2udt2);

   cout << "WaveOperator initialized" << endl;
}

void WaveOperator::Mult(const Vector &u, const Vector &du_dt,
                        Vector &d2udt2)  const
{
   // Compute:
   //    d2udt2 = M^{-1}*-K(u)
   // for d2udt2

//   Kmat.Mult(u, z);
//   z.Neg(); // z = -z
//   M_solver.Mult(z, d2udt2);

  GridFunction d2udt2_gf(&fespace);
  d2udt2_gf = 0.;
  d2udt2_bdr->SetTime(previous_time);
  Array<int> ess_bdr_ = this->ess_bdr;
  d2udt2_gf.ProjectBdrCoefficient(*d2udt2_bdr, ess_bdr_);
  for (int i = 0; i < ess_tdof_list.Size(); i++)
  {
     d2udt2[ess_tdof_list[i]] = d2udt2_gf[ess_tdof_list[i]]; // is this correct ?
  }

}

void WaveOperator::ImplicitSolve(const double fac0, const double fac1,
                                 const Vector &u, const Vector &dudt, Vector &d2udt2)
{
   if (!T)
   {
      T = Add(1.0, Mmat, fac0, Kmat);
      T_solver.SetOperator(*T);
   }

   // Get current_dt
   current_time = this->GetTime();
   current_dt = current_time - previous_time;

//   Kmat0.Mult(z, )
   Kmat.Mult(u, z);
   z.Neg();

   // Eliminate  -K_{IB}(u_B + dt dudt_B + (0.5-beta)(dt)^2 d2udt2_B) at t
   GridFunction u_gf(&fespace);
   u_bdr->SetTime(previous_time);
   u_gf = 0.;
   u_gf.ProjectBdrCoefficient(*u_bdr, this->ess_bdr);

   GridFunction dudt_gf(&fespace);
   dudt_bdr->SetTime(previous_time);
   dudt_gf = 0.;
   dudt_gf.ProjectBdrCoefficient(*dudt_bdr, this->ess_bdr);

   GridFunction d2udt2_gf(&fespace);
   d2udt2_bdr->SetTime(previous_time);
   d2udt2_gf = 0.;
   d2udt2_gf.ProjectBdrCoefficient(*d2udt2_bdr, this->ess_bdr);

   // u_B
   K->EliminateVDofsInRHS(ess_tdof_list, u_gf, z);
   // dt*dudt_B
   K->SpMatElim().AddMult(dudt_gf, z, -current_dt);
   K->SpMat().PartMult(ess_tdof_list, dudt_gf, z);
   // (0.5 - beta) * (dt)^2
   double beta = fac0 / (current_dt * current_dt);
   double fac_d2udt2 = (0.5 - beta) *(current_dt * current_dt);
   K->SpMatElim().AddMult(d2udt2_gf, z, -fac_d2udt2);
   K->SpMat().PartMult(ess_tdof_list, d2udt2_gf, z);

   // Eliminate -T_{IB}*d2udt2_B at t+dt
   // T = M + beta*(dt)^2*K
   d2udt2_bdr->SetTime(current_time);
   d2udt2_gf = 0.;
   d2udt2_gf.ProjectBdrCoefficient(*d2udt2_bdr, this->ess_bdr);
   M->EliminateVDofsInRHS(ess_tdof_list, d2udt2_gf, z);
   K->SpMatElim().AddMult(d2udt2_gf, z, -fac0);
   K->SpMat().PartMult(ess_tdof_list, d2udt2_gf, z);

   T_solver.Mult(z, d2udt2);
   
   // set Dirichlet BC for u, dudt, d2udt2 here
   // d2udt2:
   for (int i = 0; i < ess_tdof_list.Size(); i++)
   {
      d2udt2[ess_tdof_list[i]] = d2udt2_gf[ess_tdof_list[i]]; // is this correct ?
   }
   previous_time = current_time;
}

void WaveOperator::SetParameters(const Vector &u)
{
    if (0) // disable since dt is constant, so T never change
    {
      delete T;
      T = NULL; // re-compute T on the next ImplicitSolve
    }
}

WaveOperator::~WaveOperator()
{
   delete T;
   delete M;
   delete K;
   delete c2;
}

int main(int argc, char *argv[])
{
   // 1. Parse command-line options.
   const char *mesh_file = "../data/star.mesh";
   const char *ref_dir  = "";
   int ref_levels = 2;
   int order = 2;
   int ode_solver_type = 10;
   double t_final = 0.5;
   double dt = 1.0e-2;
   double speed = 1.0;
   bool visualization = true;
   bool visit = true;
   bool dirichlet = true;
   int vis_steps = 1;

   int precision = 8;
   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, "-o", "--order",
                  "Order (degree) of the finite elements.");
   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
                  "ODE solver: [0--10] - GeneralizedAlpha(0.1 * s),\n\t"
                  "\t   11 - Average Acceleration, 12 - Linear Acceleration\n"
                  "\t   13 - CentralDifference, 14 - FoxGoodwin");
   args.AddOption(&t_final, "-tf", "--t-final",
                  "Final time; start time is 0.");
   args.AddOption(&dt, "-dt", "--time-step",
                  "Time step.");
   args.AddOption(&speed, "-c", "--speed",
                  "Wave speed.");
   args.AddOption(&dirichlet, "-dir", "--dirichlet", "-neu",
                  "--neumann",
                  "BC switch.");
   args.AddOption(&ref_dir, "-r", "--ref",
                  "Reference directory for checking final solution.");
   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(&vis_steps, "-vs", "--visualization-steps",
                  "Visualize every n-th timestep.");
   args.Parse();
   if (!args.Good())
   {
      args.PrintUsage(cout);
      return 1;
   }
   args.PrintOptions(cout);

   // 2. Read the mesh from the given mesh file. We can handle triangular,
   //    quadrilateral, tetrahedral and hexahedral meshes with the same code.
   Mesh *mesh = new Mesh(mesh_file, 1, 1);
   int dim = mesh->Dimension();

   // 3. Define the ODE solver used for time integration. Several second order
   //    time integrators are available.
   SecondOrderODESolver *ode_solver;
   switch (ode_solver_type)
   {
      // Implicit methods
      case 0: ode_solver = new GeneralizedAlpha2Solver(0.0); break;
      case 1: ode_solver = new GeneralizedAlpha2Solver(0.1); break;
      case 2: ode_solver = new GeneralizedAlpha2Solver(0.2); break;
      case 3: ode_solver = new GeneralizedAlpha2Solver(0.3); break;
      case 4: ode_solver = new GeneralizedAlpha2Solver(0.4); break;
      case 5: ode_solver = new GeneralizedAlpha2Solver(0.5); break;
      case 6: ode_solver = new GeneralizedAlpha2Solver(0.6); break;
      case 7: ode_solver = new GeneralizedAlpha2Solver(0.7); break;
      case 8: ode_solver = new GeneralizedAlpha2Solver(0.8); break;
      case 9: ode_solver = new GeneralizedAlpha2Solver(0.9); break;
      case 10: ode_solver = new GeneralizedAlpha2Solver(1.0); break;

      case 11: ode_solver = new AverageAccelerationSolver(); break;
      case 12: ode_solver = new LinearAccelerationSolver(); break;
      case 13: ode_solver = new CentralDifferenceSolver(); break;
      case 14: ode_solver = new FoxGoodwinSolver(); break;
      case 15: ode_solver = new NewmarkSolver(); break;

      default:
         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
         delete mesh;
         return 3;
   }

   // 4. Refine the mesh to increase the resolution. In this example we do
   //    'ref_levels' of uniform refinement, where 'ref_levels' is a
   //    command-line parameter.
   for (int lev = 0; lev < ref_levels; lev++)
   {
      mesh->UniformRefinement();
   }
   cout << "Number of elements = " << mesh->GetNE() << endl;

   // 5. Define the vector finite element space representing the current and the
   //    initial temperature, u_ref.
   H1_FECollection fe_coll(order, dim);
   FiniteElementSpace fespace(mesh, &fe_coll);

   int fe_size = fespace.GetTrueVSize();
   cout << "Number of temperature unknowns: " << fe_size << endl;

   GridFunction u_gf(&fespace);
   GridFunction dudt_gf(&fespace);

   // 6. Set the initial conditions for u. All boundaries are considered
   //    natural.
   FunctionCoefficient u_0(InitialSolution);
   u_gf.ProjectCoefficient(u_0);
   Vector u;
   u_gf.GetTrueDofs(u);

   FunctionCoefficient dudt_0(InitialRate);
   dudt_gf.ProjectCoefficient(dudt_0);
   Vector dudt;
   dudt_gf.GetTrueDofs(dudt);

   // 7. Initialize the conduction operator and the visualization.
   Array<int> ess_bdr;
   if (mesh->bdr_attributes.Size())
   {
      ess_bdr.SetSize(mesh->bdr_attributes.Max());

      if (dirichlet)
      {
         ess_bdr = 1;
          cout << "Size of ess_bdr = " << ess_bdr.Size() << endl;
      }
      else
      {
         ess_bdr = 0;
      }
   }


   u_gf.SetFromTrueDofs(u);
   {
      ofstream omesh("ex23.mesh");
      omesh.precision(precision);
      mesh->Print(omesh);
      ofstream osol("ex23-init.gf");
      osol.precision(precision);
      u_gf.Save(osol);
      dudt_gf.Save(osol);
   }

   VisItDataCollection visit_dc("Example23", mesh);
   visit_dc.RegisterField("solution", &u_gf);
   visit_dc.RegisterField("rate", &dudt_gf);
   if (visit)
   {
      visit_dc.SetCycle(0);
      visit_dc.SetTime(0.0);
      visit_dc.Save();
   }

   socketstream sout;
   if (visualization)
   {
      char vishost[] = "localhost";
      int  visport   = 19916;
      sout.open(vishost, visport);
      if (!sout)
      {
         cout << "Unable to connect to GLVis server at "
              << vishost << ':' << visport << endl;
         visualization = false;
         cout << "GLVis visualization disabled.\n";
      }
      else
      {
         sout.precision(precision);
//         sout << "solution\n" << *mesh << dudt_gf;
         sout << "solution\n" << *mesh << u_gf;
         sout << "pause\n";
         sout << flush;
         cout << "GLVis visualization paused."
              << " Press space (in the GLVis window) to resume it.\n";
      }
   }

   // 8. Perform time-integration (looping over the time iterations, ti, with a
   //    time-step dt).
   double t = 0.0;
   WaveOperator oper(fespace, ess_bdr, speed, /*init_time*/t);
   ode_solver->Init(oper);
   bool last_step = false;

   Array<int> ess_tdof_list_main; // For DBC of u and dudt
   fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list_main);
   FunctionCoefficient u_bdr(BoundaryValue_u);
   FunctionCoefficient dudt_bdr(BoundaryValue_dudt);

   for (int ti = 1; !last_step; ti++)
   {

      if (t + dt >= t_final - dt/2)
      {
         last_step = true;
      }

      // Assign DBC for u and dudt
      u_bdr.SetTime(t);
      u_gf.ProjectBdrCoefficient(u_bdr, ess_bdr);
      dudt_bdr.SetTime(t);
      dudt_gf.ProjectBdrCoefficient(dudt_bdr, ess_bdr);
      for (int dofi = 0; dofi < ess_tdof_list_main.Size(); dofi++)
      {
         // I think directly set u is not legal
         u[ess_tdof_list_main[dofi]] = u_gf[ess_tdof_list_main[dofi]];
         dudt[ess_tdof_list_main[dofi]] = dudt_gf[ess_tdof_list_main[dofi]];
      }

      // Initial DBC for d2udt2 should set in f->Mult(x, dxdt, d2xdt2)
      ode_solver->Step(u, dudt, t, dt);

      if (last_step || (ti % vis_steps) == 0)
      {
         cout << "step " << ti << ", t = " << t << endl;

         // Assign DBC at t+dt
         u_bdr.SetTime(t);
         u_gf.ProjectBdrCoefficient(u_bdr, ess_bdr);
         dudt_bdr.SetTime(t);
         dudt_gf.ProjectBdrCoefficient(dudt_bdr, ess_bdr);
         for (int dofi = 0; dofi < ess_tdof_list_main.Size(); dofi++)
         {
            // I think directly set u is not legal
            u[ess_tdof_list_main[dofi]] = u_gf[ess_tdof_list_main[dofi]];
            dudt[ess_tdof_list_main[dofi]] = dudt_gf[ess_tdof_list_main[dofi]];
         }

         u_gf.SetFromTrueDofs(u);
         dudt_gf.SetFromTrueDofs(dudt);
         if (visualization)
         {
            sout << "solution\n" << *mesh << u_gf << flush;
         }

         if (visit)
         {
            visit_dc.SetCycle(ti);
            visit_dc.SetTime(t);
            visit_dc.Save();
         }
      }
      oper.SetParameters(u);
   }

   // 9. Save the final solution. This output can be viewed later using GLVis:
   //    "glvis -m ex23.mesh -g ex23-final.gf".
   {
      ofstream osol("ex23-final.gf");
      osol.precision(precision);
      u_gf.Save(osol);
      dudt_gf.Save(osol);
   }

   // 10. Free the used memory.
   delete ode_solver;
   delete mesh;

   return 0;
}

I use mesh uniform refined 7 times just like the example, but result seems differ.

The text was updated successfully, but these errors were encountered:

Wayne901 changed the title ~~A example of wave equation with time-dependent inhomogeneous Dirichlet boundary condition~~ An example of wave equation with time-dependent inhomogeneous Dirichlet boundary condition Mar 9, 2024

IdoAkkerman linked a pull request Apr 6, 2024 that will close this issue

Small correction to example 23 [mod-ex23] #4233

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

An example of wave equation with time-dependent inhomogeneous Dirichlet boundary condition #4182

An example of wave equation with time-dependent inhomogeneous Dirichlet boundary condition #4182

Wayne901 commented Mar 9, 2024 •

edited

An example of wave equation with time-dependent inhomogeneous Dirichlet boundary condition #4182

An example of wave equation with time-dependent inhomogeneous Dirichlet boundary condition #4182

Comments

Wayne901 commented Mar 9, 2024 • edited

Wayne901 commented Mar 9, 2024 •

edited