Nonhomogeneous Dirichlet Boundary in Time-dependent Solver

[Liu0813]

In most examples of MFEM, the homogeneous dirichlet boundary is used. But I want to set nonhomogeneous Dirichlet boundary in a time-dependent solver. I hope to know what I should do. And is there any similar examples? Thanks.

[psocratis]

Hi @Liu0813 ,

You can define a time-dependent function coefficient and project it to a GridFunction which would hold the non-homogeneous Dirichlet data. Take a look at coefficient.hpp

mfem/fem/coefficient.hpp

Lines 127 to 147 in b3beafe

	/// A general C-function coefficient
	class FunctionCoefficient : public Coefficient
	{
	protected:
	double (*Function)(const Vector &);
	double (*TDFunction)(const Vector &, double);
	public:
	/// Define a time-independent coefficient from a pointer to a C-function
	FunctionCoefficient(double (*f)(const Vector &))
	{
	Function = f;
	TDFunction = NULL;
	}
	/// Define a time-dependent coefficient from a pointer to a C-function
	FunctionCoefficient(double (*tdf)(const Vector &, double))
	{
	Function = NULL;
	TDFunction = tdf;
	}

Thanks,

Socratis

[Liu0813]

Thanks for your reply. I use the function ProjectBdrCoefficient to project the Coefficient on the GridFuntion to be solved. Is that right? By the way, do I need to add these boundary degrees of freedom to the ess_tdof_list?

[psocratis]

Hi @Liu0813,

Yes, for an H1 space you can use ProjectBdrCoefficient, and yes you have to construct the list of essential true dofs by fespace.GetEssentialTrueDofs(ess_bdr,ess_tdof_list) and then use it when you form the linear system.

Best,

Socratis

[Liu0813]

Hi @psocratis,
Thank you for your suggestion. I try to solve my problem using the method you mentioned, but it seems that there is a bug.
When I put the degrees of freedom on the non-homogeneous Dirichlet boundary in the list of essential true dofs ess_tdof_list and use it to form the linear system, the sparsity patterns of mass matrix and stiffness matrix is shown as the following figure.

mfem/fem/bilinearform.cpp

Lines 956 to 976 in 4f521a1

	void BilinearForm::EliminateVDofs(const Array<int> &vdofs,
	DiagonalPolicy dpolicy)
	{
	if (mat_e == NULL)
	{
	mat_e = new SparseMatrix(height);
	}
	for (int i = 0; i < vdofs.Size(); i++)
	{
	int vdof = vdofs[i];
	if ( vdof >= 0 )
	{
	mat -> EliminateRowCol (vdof, *mat_e, dpolicy);
	}
	else
	{
	mat -> EliminateRowCol (-1-vdof, *mat_e, dpolicy);
	}
	}
	}

That is, all the corresponding rows and columns are eliminated when the function FormSystemMatrix(ess_tdof_list, Mmat) is performed, as the following code in bilinearform.cpp (release version mfem-4.1)

There is no problem for homogeneous Dirichlet boundary, but for non-homogeneous Dirichlet boundary, I don't think this treatment is reasonable.

I am a beginner of MFEM, so I am not sure of this problem. Should I change the right-side item of the linear system non-homogeneous Dirichlet boundary?

Thanks,
Liu

[psocratis]

Hi @Liu0813,

Suppose you have BilinearForm * a. Then, in order to form the linear system, you need to call
a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);. This method modifies both the system matrix (by calling the method FormSystemMatrix) and the right-hand side. The essential dofs (truedof indices) are read from the Array<int> ess_tdof_list and the essential data from the GridFunction x.

Please let us know if this solves the issue.

Best,

Socratis

[Liu0813]

Hi @psocratis,

That is right. The method FormLinearSystem modifies both the system matrix and the right-hand side. However, my problem to solve is time-dependent. After spacial discretization, the partial different equation should be writen as an ODE: du/dt=M^{-1}(-Ku). For the implicit solver, the Backward-Euler equation k=f(u+dt*k,t) should be solved. Take a look at the Example 16. According to the example, the variable du/dt is firstly solved in time-dependent problem. So, I am not sure if the method FormLinearSystem is unsuitable for this case. And I hope to know how to use it in time-dependent problem.

Thanks,
Liu

[cjvogl]

@Liu0813, you can still use FormLinearSystem for time-dependent essential data. The steps to solve k=f(u+dt*k,t) could be something like:

1. use your time-dependent FunctionCoefficient to project the correct boundary data for time t onto x
2. call FormLinearSystem using the updated x containing the essential boundary data
3. solve AX=B

A disclaimer is that you'll need to be a little more careful if you want higher-than-first order accuracy in time... but I believe this will work fine for forward/backward Euler.

I can also elaborate on my suggested approach if needed.

[Liu0813]

Hi @cjvogl,
Thank you for your reply. I still have some questions about your suggested approach. Considering a parabolic model problem, such as the unsteady heat equation with time-dependent essential boundary conditions, the semi-discrete form can be expressed as

According to your suggested approach, I should obtain the full-discrete form. For backward Euler method, it is

and the linear system can be writen as

That means the sparse matrix is A=M+dt*K, and the right-hand side is M*u. I have to construct the mass matrix M and stiffness matrix K before forming the linear system. Thus, the questions are

1. When M and K are constructed, Do I need to introduce the ess_tdof_list?
2. In the standard ODE solver of MFEM, the Step method calls ImplicitSolve method to calculate du/dt. You suggested approach does not seem to match with this process. Do I need to override the ODESolver?

Thanks,
Liu

[v-dobrev]

Hi @Liu0813,

Here are some notes about applying b.c. that might be helpful.

If you split the dofs in to two groups: I (for interior) and B (for essential boundary) then the matrices and the solution vector can be written in the following block form:

M = | M_I  M_IB |   K = | K_I  K_IB |   u = | u_I |
    | M_BI M_B  |       | K_BI K_B  |       | u_B |

Since we want to solve just for the I block, the real system of ODEs has the form:

M_I du_I/dt = - K_I u_I - K_IB u_B

Since we do not want to change the sizes of the vectors and matrices depending on the essential b.c., we can augment the above system of ODEs by

    dx_B/dt = 0

and solve for the vector:

x = | u_I |
    | x_B |

e.g. with initial condition x(0) = u(0). Note that when we need to save or visualize the solution at time t we'll need to replace the x_B part with the real u_B(t). With this modification, the block system of ODEs can be written as:

dx/dt = | M_I 0 |^{-1} | - K_I u_I - K_IB u_B |
        |  0  I |      |   0                  |

Note that when FormLinearSystem is called the returned matrices have the form:

MM = | M_I 0 |   KM = | K_I 0 |
     | 0   I |        | 0   I |

and the complement matrices MC = M - MM and KC = K - KM are stored internally, so that they can be used to eliminate boundary conditions. The elimination modifies the r.h.s. vector b in the following way:

b = | b_I | --> b = | b_I - K_IB u_B |
    | b_B |         | u_B            |

where u_B is the essential b.c. taken from the input solution vector. In Example 16, this means that inside ConductionOperator::Mult(...), before the call to M_solver.Mult(...), one can call K->FormLinearSystem(...) to transform the vector z by subtracting K_IB u_B from z_I. Then after M_solver.Mult(...), one can simply reset the boundary component of the output vector, du_dt, by calling du_dt.SetSubVector(ess_tdof_list, 0.0).

To support SDIRK methods (and backward Euler, in particular) one has to update the implementation of ConductionOperator::ImplicitSolve too. In this case the system for k (i.e. du_dt) can be written as

MM | k_I | = - | K_I 0 | | u_I+dt*k_I | - | K_BI u_B(t) |
   | k_B |     | 0   I | | 0          |   | 0           |

or

| M_I+dt*K_I  0 | | k_I | = - | K_I 0 | | u_I | - | K_BI u_B(t) |
| 0           I | | k_B |     | 0   I | | 0   |   | 0           |

Note that the matrix of the l.h.s. can also be replaced by MM+dt*KM -- this is what is already done in Example 16. The modification for nonzero u_B here is very similar to the one in ConductionOperator::Mult.

[Liu0813]

Hi @v-dobrev,
Thanks for your guidance. Everything is clear.

[chang-yiyi]

Hi @v-dobrev ,

If the boundary condition is nonhomogeneous Dirichlet and time-dependent, i.e., dx_B/dt != 0. It seems not simple.
Before set essential boundary conditions, we get

| M_I  M_IB | | dxI_dt | = - |  K_I  u_I + K_IB u_B |
| M_BI  M_B | | dxB_dt |     |  K_BI u_I +  K_B u_B |

In terms of the above explanations, after set essential boundary conditions, we should get

| M_I  0 | | dxI_dt | =  | - K_I u_I - K_IB u_B - M_IB dxB_dt|
| 0    I | | dxB_dt |    |               dxB_dt              |

It's not clear to me how to obtain the rhs?

[stale]

⚠️ This issue or PR has been automatically marked as stale because it has not had any activity in the last month. If no activity occurs in the next week, it will be automatically closed. Thank you for your contributions.

[v-dobrev]

Hi @chang-yiyi,

Sorry for not replying earlier.

You are correct: when du_B/dt != 0 one has the additional term - M_IB du_B/dt in the r.h.s. One way to deal with this term (when your Dirichlet b.c. is given analytically) is to differentiate the Dirichlet b.c. function with respect to time, define a separate coefficient based on that derivative and then project it on the boundary to obtain the vector du_B/dt. The multiplication with M_IB can be done the same way as the multiplication with K_IB described above.

If your u_B vector is only given at discrete times, one can probably use a backward differentiation formula (BDF) like backward-Euler or higher-order approximations to compute an approximation of du_B/dt.

Another possibility (I'm not sure if this works well) is to use under-integration for the mass matrix, so that it becomes diagonal and then M_IB = 0.

[chang-yiyi]

Thanks @v-dobrev ,

Your answer is very clear. Thanks again.

Making a small suggestion that if replacing homogenous boundary conditions into inhomogeneous ones and replacing only Dirichlet boundary conditions into mixed Dirichlet and Neumann boundary conditions in examples in mfem/examples maybe could solve many issues. It's just my thinking. Sorry if it's rude to ask for a lot of demands.

[tzanio]

Hi @chang-yiyi,

Thank you for the feedback, that is useful and definitely appreciated 👍

Yes, we should have more boundary condition in the examples, that's why we added Example 27, and we should continue to add more. We always welcome contributions, so if you have any specific code suggestions, please feel free to open a PR.

Tzanio

[tangqi]

Hi @v-dobrev
I have a question related to time-dependent Dirichlet following your suggestion.

If I am solving the system fully implicitly, and use backward Euler, I can let dxB_dt=(u_B^{n+1}-u_B^n)/dt. I think it should be fine to get first order. Can I use such an approach to extend it to high-order using SDIRK, assuming u_B^{n+1} is always updated analytically even in the intermediate stages? Or should I always construct an analytical dxB_dt?

Due to my preconditioner, I prefer to think the problem in a fully discretized sense. So I always target at solving u^{n+1} in my reduced system, and then recover the backward Euler update through (u^{n+1}-u^n)/dt. So I do not fully understand the impact of dxB_dt in a high-order integrator

Qi

[mlstowell]

Regarding the need for more examples, I fully agree that more examples would be helpful. However, there are so many special cases that this is difficult. For one example you could look at the maxwell mini app in the electromagnetics subdirectory. This miniapp uses a variable order ODE integrator which requires implicit solves using intermediate time steps of differing lengths. Since these time steps repeat cyclically we build a separate linear system for each distinct time step and reuse them as needed. This approach, therefore, uses additional memory to avoid rebuilding the linear system. This may not be optimal but it demonstrates one approach which, I think, is well suited to this particular problem.

Another approach can be seen in the joule miniapp which uses a time dependent voltage source. In this case I don't think the linear system is implicit (I don't even think the voltage update is done as a time evolution) but it may be useful.

[stefanozampini]

@chang-yiyi the PETSc classes can handle time-dependent non-homogeneous boundary conditions (linearity in state vector is assumed), see e.g https://github.com/mfem/mfem/blob/master/linalg/petsc.cpp#L3659. This function computes the value of the residual of the ODE written in implicit form using exact BC. Then, it modifies the residual to properly compute the update when solving the (non)linear system of equations (take a look at how the Jacobian is modified here https://github.com/mfem/mfem/blob/master/linalg/petsc.cpp#L3760). The interface for boundary conditions is quite general https://github.com/mfem/mfem/blob/master/linalg/petsc.hpp#L452. You only need to override the evaluation method https://github.com/mfem/mfem/blob/master/linalg/petsc.hpp#L478.

[lmolin3]

Hi @chang-yiyi,

Sorry for not replying earlier.

You are correct: when du_B/dt != 0 one has the additional term - M_IB du_B/dt in the r.h.s. One way to deal with this term (when your Dirichlet b.c. is given analytically) is to differentiate the Dirichlet b.c. function with respect to time, define a separate coefficient based on that derivative and then project it on the boundary to obtain the vector du_B/dt. The multiplication with M_IB can be done the same way as the multiplication with K_IB described above.

If your u_B vector is only given at discrete times, one can probably use a backward differentiation formula (BDF) like backward-Euler or higher-order approximations to compute an approximation of du_B/dt.

Another possibility (I'm not sure if this works well) is to use under-integration for the mass matrix, so that it becomes diagonal and then M_IB = 0.

Hi @v-dobrev ,
Sorry to comment on this old issue.
I implemented a heat transfer solver with Neumann, Robin, and Dirichlet time-dependent conditions.
For the latter I followed the suggestions and computed a BDF approximation of du_B/dt with order increasing up to 6 (and consistent with the order of the time integration scheme).
In the solver I set du/dt (ess_dofs) = du_B/dt and u(ess_dofs) = u_n so that after the step i correctly get u(ess_dofs) = u_n + du_B/dt = u_B.
I tested the code running a convergence analysis and everything seems to work.

What I was wondering is: can we instead set du/dt (ess_dofs) = 0 and u(ess_dofs) = u_B. This would satisfy u(ess_dofs) = u_B, but I was wondering if it's doable or would compromise anything in the rest of the code.

Also, for higher order schemes, should the derivative approximation be computed for each intermediate step when calling the overridden SetTime() method?

Thanks!