Update boundary conditions in a time dependent problem without reassembling matrix with multigrid solver

[Markus Stricker]

Hi,

i haven't found a solution to a problem I'm working on right now. Neither in the tutorials nor in the FAQ nor in this user group.

I am trying to solve a time dependent problem with the multigrid approach. Each new time step requires to update the boundary conditions (Dirichlet and Neumann BC). The nodes at which these BC are applied do not change over time, only their corresponding value is updated.

As far as my initial test go, each time I change the value I have to go through the process of reassembling the matrix, which makes my calculation unusably slow (I need to perform on the order of 25000 time steps or more).

Is there a way to update the BC within the multigrid environment without reassembling the matrix every time step?

Thanks in advance,
Markus

[Martin Kronbichler]

Dear Markus,

> I am trying to solve a time dependent problem with the multigrid
> approach. Each new time step requires to update the boundary
> conditions (Dirichlet and Neumann BC). The nodes at which these BC are
> applied do not change over time, only their corresponding value is
> updated.
>
> As far as my initial test go, each time I change the value I have to
> go through the process of reassembling the matrix, which makes my
> calculation unusably slow (I need to perform on the order of 25000
> time steps or more).

I guess your concern are the Dirichlet boundary conditions, right? There
are at least three ways to deal with this case:
a) You hold a two global matrices. You assemble into the first one
before the time loop. In the time loop, you copy the matrix into another
sparse matrix where you call "apply_boundary_values" (or the respective
ConstraintMatrix::condense call with matrix and vector if you have your
boundary conditions in a ConstraintMatrix).
b) You identify cells where inhomogeneous boundary conditions are,
assemble a matrix there also in the time loop, and pass this matrix to
the distribute_local_to_global call of the constraint matrix (where the
local matrix is the fourth argument). This provides the necessary
columns of the inhomogeneously constrained DoFs that enter the rhs
vector.
c) You reformulate your linear system so that you actually solve for
increments instead of the solution itself. This reflects to first
applying the current values to the solution components at the boundary,
adding terms similar to the matrix entries to the right hand side of
your weak form (the same way as you would compute a residual of your
differential operator e.g. with Newton's method), and solving a linear
system with homogeneous boundary conditions.

Which version to prefer depends on your application. I rarely advise to
use variant a) because it considerably increase matrix memory and
requires care when used in parallel with MPI. Variant b) is the simplest
way because all you need to add is a few checks if the current cell is
at a boundary of interest and then you go through the usual matrix
assembly code that you have already. Variant c) requires new code
(unless you use Newton's method) but is most flexible.

Have also a look at
http://www.dealii.org/developer/doxygen/deal.II/group__constraints.html
(section 'Treatment of inhomogeneous constraints').

Best,
Martin

[Praveen C]

I have seen the use of FilteredMatrixin such situations, see Session 7 here

http://numerik.uni-hd.de/~lehre/SS10/dealcompact/

How does this method compare with other approaches listed by Martin in terms of efficiency and parallelization ?

Thanks

praveen

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Markus Stricker]

Martin,

thanks for your quick response. I'll look into that and get back here with the solution I applied.

Best,
Markus

[Martin Kronbichler]

Hi Praveen,

> I have seen the use of FilteredMatrixin such situations

> How does this method compare with other approaches listed by Martin in
> terms of efficiency and parallelization ?

I have not used FilteredMatrix before (as I said, there are several
alternatives). I would not use it in parallel because it needs to set
elements of the vector in the methods pre_filter and post_filter (check
the implementation in include/deal.II/lac/filtered_matrix.h), but due to
the 'generic' vector type it cannot do the data exchange. It should be
fine anyway but I have not checked all details. Regarding efficiency, it
should be similarly fast as variant a) described by me above from a
matrix point of view, but without the memory overhead. The solution
phase will be a bit slower from what you can get from variants b) and c)
where you may drop the matrix entries over constrained entries, reducing
the matrix nonzeros by a few percent. (You set that option when creating
the sparsity pattern.) However, variants b) and c) involve additional
evaluations during assembly, which may increase total cost depending on
the application. I typically prefer faster solvers and somewhat more
work to do during assembly because solvers take often most of the time.

A more important restriction with FilteredMatrix is the choice of
preconditioners, which you should check carefully.

Best,
Martin

[Denis Davydov]

Approach (b) worked fine for me in a similar context where Dirichlet BC are changing 

and the matrix stays the same.

Regards,

Denis.