How to set non-zero BCs to TimeDependentOperator

[colmenaresj]

I am trying to take something like the ConductionOperator in Example 16 and apply non-zero BCs (both Neumann and Dirichlet), but so far I haven't had any success. None of the examples show how to do this for Time-dependent problems. Could someone please show an example of how to do this?

I just want to do a simple heat conduction test based on Ex 16 (with alpha=0), using the beam-tet.mesh geometry, where I fix the temperature on one end and apply a heat flux (neumann) on the other end. I tried simply adding a LinearForm object and adding it to the du_dt but it doesn't work. I tried doing something like K.AddBoundaryIntegrator, but I have no idea how to use that. Any guidance is welcome.

Thanks, -JD

[colmenaresj]

I've figured out how to do Dirichlet. I will post a code snippet soon. However, I am still struggling with Neumann bcs.

[colmenaresj]

Me being the dummy that I am, I had forgotten to do b->Assemble(), so I wasn't getting anything for the Neumann bc. I think I managed to do it now.

[colmenaresj]

Here is my what my ImplicitSolve and Mult operators looke like:

void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
{
   // Compute:
   //    du_dt = M^{-1}*-Ku
   // for du_dt, where K is linearized by using u from the previous timestep
   Kmat0.Mult(u, z);
   z.Neg(); // z = -z
   z += (*b); // b is a pointer, so need to de-reference using *
   M_solver.Mult(z, du_dt);
   du_dt.SetSubVector(ess_tdof_list,0.0);
}

void ConductionOperator::ImplicitSolve(const real_t dt,
                                       const Vector &u, Vector &du_dt)
{
   // Solve the equation:
   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
   // for du_dt, where K is linearized by using u from the previous timestep
   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
   Kmat0.Mult(u, z);
   z.Neg();
   z += (*b); // b is a pointer, so need to de-reference using *
   T_solver.Mult(z, du_dt);

   du_dt.SetSubVector(ess_tdof_list,0.0);
}

Note that in some parts, I am using Kmat0 which created using an empty ess_tdof_list. On the other hand, Kmat was created with the actual ess_tdof_list. I don't know if this is the right approach, but I took it from Example 23.

For Neumman, I just set the ParLinearForm* b object as you normally would, and added it to the z temporary array.

For enforcing Dirichlet BC, I force du_dt=0 at the dofs given by ess_tdof_list using du_dt.SetSubVector(ess_tdof_list,0.0);. Again, this worked for me but I don't know if there is a better/more appropriate approach.

[tangqi]

If your Kmat and Mmat are created with essential dofs eliminated, this looks correct to me.

I have not done non-zero Neumann. Sometimes, I found it is easy to solve u^{n+1} first in the ImplicitSolve and then construct du_dt.

[colmenaresj]

Thanks for the reply. I passed an empty dofs list to Kmat0, but the populated tdof list to Kmat and Mmat. Was that wrong?

   fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
   M = new ParBilinearForm(&fespace);
   M->AddDomainIntegrator(new MassIntegrator());
   M->Assemble(0); // keep sparsity pattern of M and K the same
   M->FormSystemMatrix(ess_tdof_list, Mmat);
   K = new ParBilinearForm(&fespace);
...
   GridFunctionCoefficient u_coeff(&u_alpha_gf);
   K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
   K->Assemble(0); // keep sparsity pattern of M and K the same
   Array<int> dummy;
   K->FormSystemMatrix(dummy, Kmat0);
   K->FormSystemMatrix(ess_tdof_list, Kmat);

[colmenaresj]

I saw this approach in Example 23, so it seemed correct.

[v-dobrev]

Hi @colmenaresj,

There was a discussion on this topic some time ago here: #1720; my comments can be found here #1720 (comment), and here #1720 (comment).

Another more recent related issues is #4182 which is about non-homogeneous b.c. for the wave equation.

It seems that you already got things working. Is that correct, or do you still need help resolving some issue?

[colmenaresj]

Thank you! I had looked for previous posts in the Discussion forum but did not think of checking the issues. Best,-JD

[colmenaresj]

@v-dobrev It is working, I just wanted to double-check if the approach of using an empty list for Kmat0 only and a full list for Kmat and Mmat was appropriate. I will look at your comments in the issue you shared, you seem to explain that there. Thanks.