mfem / mfem

Lightweight, general, scalable C++ library for finite element methods
http://mfem.org
BSD 3-Clause "New" or "Revised" License
1.63k stars 472 forks source link

Best example to start with for time-dependent conduction problem? #1634

Closed ericman314 closed 3 years ago

ericman314 commented 3 years ago

I am new to MFEM. I'm trying to figure out how to create a time-dependent heat simulation with Neumann boundary conditions. My input is a 2D tri-mesh. The initial condition is constant temperature throughout, and then there should be a period of constant non-zero heat flux through the boundary, followed by a period of no heating (insulated boundary).

I've looked at ex16.cpp, since it seemed like it was pretty similar to my problem. I got it to compile and run, and even modified it to add Dirichlet boundary conditions following ex1.cpp. But I can't figure out how to add the Neumann boundary conditions. I also don't need the temperature-dependent conduction coefficient, and I think the problem should be much simpler than what's implemented in ex16.

I'm really impressed by what this library can do. I'm just frustrated because this is one of the simplest problems and I just can't figure out where to start. Thanks in advance.

v-dobrev commented 3 years ago

You can take a look at example 27 which illustrates different types of boundary conditions. Specifically, the Neumann b.c. is applied here: https://github.com/mfem/mfem/blob/0effa7abffb019f087ef4a84a9ae0c917c29e6da/examples/ex27.cpp#L243-L244

In a time-dependent setting, if the coefficient m_nbcCoef depends on time, the linear form will need to be re-assembled when time changes.

ericman314 commented 3 years ago

Oh thanks, I didn't see that example before.

How then do we use that in a time-dependent problem? So far, example 16 is the closest I found to what I would like to do. But there is no LinearForm to add the boundary integrator to--that's where I'm getting stuck.

vladotomov commented 3 years ago

Hi Eric,

For explicit time integration, the RHS in ex 16 is constructed here: https://github.com/mfem/mfem/blob/0effa7abffb019f087ef4a84a9ae0c917c29e6da/examples/ex16.cpp#L331-L332 If you create a LinearForm (as in ex 27), you can Assemble it and add its contribution to the RHS at that point.

For implicit time integration, if your BC doesn't depend on u, you can do the same, but here: https://github.com/mfem/mfem/blob/0effa7abffb019f087ef4a84a9ae0c917c29e6da/examples/ex16.cpp#L349-L350

Hope this helps, Vladimir

ericman314 commented 3 years ago

Thank you. This is what I have so far:

void ConductionOperator::ImplicitSolve(const double dt,
                                       const Vector &u, Vector &du_dt)
{
  // Solve the equation:
  //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
  // for du_dt
  if (!T)
  {
    T = Add(1.0, Mmat, dt, Kmat);
    current_dt = dt;
    T_solver.SetOperator(*T);
  }
  MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
  Kmat.Mult(u, z); // z = K * u
  z.Neg(); // z = -z

  LinearForm n(&fespace);
  ConstantCoefficient qCoef(1.0);
  n.AddBoundaryIntegrator(new BoundaryLFIntegrator(qCoef));
  n.Assemble();

  z.Add( ?? );

  T_solver.Mult(z, du_dt);
}

But I'm not sure what is meant by adding the contribution of the LinearForm to the RHS.

vladotomov commented 3 years ago

Try z.Add(1.0, n); or simply z += n;

The class LinearForm inherits from class Vector, thus it can be treated as a Vector.

ericman314 commented 3 years ago

Oh, gotcha! I got it working. Thank you so much!

I'm simulating local heating of a metal part from laser cutting. The results are looking brilliant. Thanks for this terrific software.

image