Finite Volume Method - Matrix/RHS assembly

[Nik Markus Leuenberger]

Dear all,

I'm new to deal.ii and I would like to use it to gradually build a simulator for charge/mass transport inside batteries.

Since the equations I'm trying to solve are conservations laws, I would like to use a finite volume type of implementation. I know that deal.ii is a finite element library. However, since I would like to use a fully implicit approach and later hopefully octree refined grids, I think deal.ii offers the things needed for this task as well as many functionalities and interfaces for future tasks.

To start, I've summarized a simplified problem I would like to solve in the appended .pdf file. It's a 1D Laplace equation using finite volumes.

I would like to build on the step-12 or step-12b tutorials because they use discontinuous elements too. I don't call the beta function of the original tutorials such that a 1D version should be possible.

The appended nik-step-12.cc file contains a hard-coded version of the three assembly functions that results in the correct matrix and solution for a very specific case described in the .pdf (section 3).

My questions would therefore be the following:

• How should I best write the three functions (cell, boundary and interior face worker) for the matrix/RHS assembly with as less hard-coding as possible for a finite volume scheme like equations (7-9) in the .pdf.
  
• Related - How can I access the elements of the solution vector (phi_i) of the current cell and neighbor cell that are used for example in equation (7) in the .pdf.

Thank you very much in advance for your help or inputs.

Best,

Nik

[Abbas Ballout]

Nick, 
You can look at step 74. It has EXACTLY what you need. 

[Wolfgang Bangerth]

> Since the equations I'm trying to solve are conservations laws, I would like
> to use a finite volume type of implementation. I know that deal.ii is a finite
> element library. However, since I would like to use a fully implicit approach
> and later hopefully octree refined grids, I think deal.ii offers the things
> needed for this task as well as many functionalities and interfaces for future
> tasks.

This is true, but before you get fully settled on FV methods, may I suggest
that you also think about the accuracy you get for the simplicity? If you use
piecewise constants (=finite volumes), the best you can hope for is O(h)
accuracy. That's quite inaccurate. By using finite element methods, you will
generally be able to obtain much better accuracy (not just in terms of the
convergence *rate*, but in *absolute errors*).

> To start, I've summarized a simplified problem I would like to solve in the
> appended .pdf file. It's a 1D Laplace equation using finite volumes.

I think I understand what you want to do, but you may want to see if you can
make your writing more precise. For example, in (4), you say you integrate
over S, but you're missing the normal vector that in your case is either +1 or
-1, but is necessary to ensure that the formula is correct. You can see that
in the transition between (5) and (6) where a minus sign magically appears in
front of the first term. Of course, (5) and (6) are also not exactly equal,
but only approximately.

Second, in (6) there magically appears a factor of h^2 for which I don't have
an explanation.

(No need to write back to these comments -- I just wanted to point it out, the
teacher in me can't just let it go ;-) )

> I would like to build on the step-12 or step-12b tutorials because they use
> discontinuous elements too. I don't call the beta function of the original
> tutorials such that a 1D version should be possible.
>
> The appended nik-step-12.cc file contains a hard-coded version of the three
> assembly functions that results in the correct matrix and solution for a very
> specific case described in the .pdf (section 3).
>
> My questions would therefore be the following:
>

> * How should I best write the three functions (cell, boundary and interior

> face worker) for the matrix/RHS assembly with as less hard-coding as
> possible for a finite volume scheme like equations (7-9) in the .pdf.

> * Related - How can I access the elements of the solution vector (phi_i) of

> the current cell and neighbor cell that are used for example in equation
> (7) in the .pdf.

As already mentioned by Abbas, step-74 gives a good overview.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Nik Markus Leuenberger]

Dear Abbas, Dear Prof. Bangerth, 

 

Thank you very much for your fast replies to my post. I'm sorry for not having found step-74 myself before.

step-74 is I think indeed what I am looking for: However, when I played around with step-74, I ran into the following two issues:

 

(1) trying to solve in 1D:

• I get errors of the following type:
  error: ‘class dealii::DoFAccessor<0, 1, 1, false>’ has no member named ‘measure’
    411 | re() / cell->face(face_no)->measure();  I can make it run by removing the calls to these functions and replace the penalty factor which would involve them by a constant number, which according to the output seems to work o.k. but might need further tuning of the penalty constant (here just chosen as a constant of 100).

degree = 1
| cycle | cells | dofs | L2        | L2...red.rate.log2 | H1        | H1...red.rate.log2 | Energy    |
| 0     | 4     | 8    | 5.426e-03 | -                  | 7.226e-02 | -                  | 7.458e-02 |
| 1     | 8     | 16   | 1.360e-03 | 2.00               | 3.619e-02 | 1.00               | 3.745e-02 |
| 2     | 16    | 32   | 3.428e-04 | 1.99               | 1.818e-02 | 0.99               | 1.894e-02 |
| 3     | 32    | 64   | 8.852e-05 | 1.95               | 9.267e-03 | 0.97               | 9.894e-03 |
| 4     | 64    | 128  | 4.032e-05 | 1.13               | 7.273e-03 | 0.35               | 9.568e-03 |
| 5     | 128   | 256  | 1.461e-05 | 1.46               | 4.961e-03 | 0.55               | 6.636e-03 |
| 6     | 256   | 512  | 4.876e-06 | 1.58               | 2.073e-03 | 1.26               | 2.706e-03 |
| 7     | 512   | 1024 | 1.837e-05 | -1.91              | 8.917e-03 | -2.10              | 1.260e-02 |
| 8     | 1024  | 2048 | 4.070e-06 | 2.17               | 2.706e-03 | 1.72               | 3.817e-03 |

| 9     | 2048  | 4096 | 2.449e-07 | 4.06               | 2.652e-04 | 3.35               | 3.474e-04 | 

• This first point however is less important than (2) below, because my desired usecase is anyways either 2D or 3D, so I could adjust the boundary conditions to make the 2D essentially 1D. However, it would of course be nice to be able to solve a simple 1D problem as the one described in the .pdf first.

 (2) trying to solve for polynomial degree 0 (i.e. piecewise constant): 

• When I try to solve the default convergence_rate test case in step-74 with degree p = 0, I don't see a decrease in the L2 norm of the solution. I tested different values for the penalty_factor for which the given function would evaluate to 0 because of the p*(p+1) term. I tested just removing the p*(p+1) term and also tested constants like 1000 or 0.001, but in all cases, I got results similar to the table below, whre the L2 norm does not decrease. The table below is for the case with penalty term  = 0.5 * (1. / cell_extent_left + 1. / cell_extent_right);

degree = 0
| cycle | cells | dofs  | L2        | L2...red.rate.log2 | H1        | H1...red.rate.log2 | Energy    |
| 0     | 16    | 16    | 3.016e-01 | -                  | 4.443e+00 | -                  | 6.045e+00 |
| 1     | 64    | 64    | 2.273e-01 | 0.41               | 4.443e+00 | 0.00               | 5.680e+00 |
| 2     | 256   | 256   | 2.284e-01 | -0.01              | 4.443e+00 | 0.00               | 5.554e+00 |
| 3     | 1024  | 1024  | 2.368e-01 | -0.05              | 4.443e+00 | 0.00               | 5.497e+00 |
| 4     | 4096  | 4096  | 2.428e-01 | -0.04              | 4.443e+00 | 0.00               | 5.470e+00 |
| 5     | 16384 | 16384 | 2.462e-01 | -0.02              | 4.443e+00 | 0.00               | 5.456e+00 |
| 6     | 65536 | 65536 | 2.481e-01 | -0.01              | 4.443e+00 | 0.00               | 5.448e+00 |

 

Maybe this is well-known and it shouldn't work. However, it would be nice to be able to solve with degree 0 because that would allow me to mimic a finite volume type implemenation. 

Thank you very much again on any hints how I could resolve either of the two issues.

Best regards,

Nik

On Friday, June 21, 2024 at 12:11:32 AM UTC+2 Wolfgang Bangerth wrote:

> Since the equations I'm trying to solve are conservations laws, I would like
> to use a finite volume type of implementation. I know that deal.ii is a finite
> element library. However, since I would like to use a fully implicit approach
> and later hopefully octree refined grids, I think deal.ii offers the things
> needed for this task as well as many functionalities and interfaces for future
> tasks.

This is true, but before you get fully settled on FV methods, may I suggest
that you also think about the accuracy you get for the simplicity? If you use
piecewise constants (=finite volumes), the best you can hope for is O(h)
accuracy. That's quite inaccurate. By using finite element methods, you will
generally be able to obtain much better accuracy (not just in terms of the
convergence *rate*, but in *absolute errors*).

I agree, the discontinuous elements are very promising for this because ideally, I could (1) recover the desired finite

volume formulation for degree 0, (2) keep the conservativeness and (3) easily get better accuracy by using higher order elements.

 

> To start, I've summarized a simplified problem I would like to solve in the
> appended .pdf file. It's a 1D Laplace equation using finite volumes.

I think I understand what you want to do, but you may want to see if you can
make your writing more precise. For example, in (4), you say you integrate
over S, but you're missing the normal vector that in your case is either +1 or
-1, but is necessary to ensure that the formula is correct. You can see that
in the transition between (5) and (6) where a minus sign magically appears in
front of the first term. Of course, (5) and (6) are also not exactly equal,
but only approximately.

Second, in (6) there magically appears a factor of h^2 for which I don't have
an explanation.

(No need to write back to these comments -- I just wanted to point it out, the
teacher in me can't just let it go ;-) )

Thank you very much for the comments - helps me to be more precise.

[Wolfgang Bangerth]

Nik:

> (1) trying to solve in 1D:
>

> * I get errors of the following type:
> *error: ‘class dealii::DoFAccessor<0, 1, 1, false>’ has no member
> named ‘measure’
>   411 | re() / cell->face(face_no)->measure();* / I can make it run

> by removing the calls to these functions and replace the penalty
> factor which would involve them by a constant number, which according
> to the output seems to work o.k. but might need further tuning of the

> penalty constant (here just chosen as a constant of 100)./

Yes, that's because the measure of a face in 1d is hard to define -- what is
the measure of a point in a zero-dimensional space?

If you have a good idea for how the 'extent1/2' variables should be defined in
1d, let us know and we can patch the program so that it will also work in 1d.

>  (2) trying to solve for polynomial degree 0 (i.e. piecewise constant):
>

> * /When I try to solve the default convergence_rate test case in step-74

> with degree p = 0, I don't see a decrease in the L2 norm of the solution.
> I tested different values for the penalty_factor for which the given
> function would evaluate to 0 because of the p*(p+1) term. I tested just
> removing the p*(p+1) term and also tested constants like 1000 or 0.001,
> but in all cases, I got results similar to the table below, whre the L2
> norm does not decrease. The table below is for the case with penalty term

> /= 0.5 *(1./cell_extent_left+1./cell_extent_right);

>
> degree = 0
> | cycle | cells | dofs  | L2        | L2...red.rate.log2 | H1        |
> H1...red.rate.log2 | Energy    |
> | 0     | 16    | 16    | 3.016e-01 | -                  | 4.443e+00 | -
>            | 6.045e+00 |
> | 1     | 64    | 64    | 2.273e-01 | 0.41               | 4.443e+00 | 0.00
>             | 5.680e+00 |
> | 2     | 256   | 256   | 2.284e-01 | -0.01              | 4.443e+00 | 0.00
>             | 5.554e+00 |
> | 3     | 1024  | 1024  | 2.368e-01 | -0.05              | 4.443e+00 | 0.00
>             | 5.497e+00 |
> | 4     | 4096  | 4096  | 2.428e-01 | -0.04              | 4.443e+00 | 0.00
>             | 5.470e+00 |
> | 5     | 16384 | 16384 | 2.462e-01 | -0.02              | 4.443e+00 | 0.00
>             | 5.456e+00 |
> | 6     | 65536 | 65536 | 2.481e-01 | -0.01              | 4.443e+00 | 0.00
>             | 5.448e+00 |
>
> Maybe this is well-known and it shouldn't work. However, it would be nice to
> be able to solve with degree 0 because that would allow me to mimic a finite
> volume type implemenation.

I don't know, and I don't know whether the people who wrote step-74 have
thought about the case. One would *hope* that the method also works for the
case p=0. A good first step would be to look at how the solution looks if you
visualize it -- that is often a good start towards debugging what is going wrong.

[Abbas Ballout]

Nick, 

2)
When you use DG0, all of the terms in the weak from cancel to zero except for the face term penallty*jump(u)*jump(v).  (du/dx and dv/dy are zero inside an element)
That penalty term in the DG community is seen as something that adds stability and isn't at the core of the scheme.  

In FV after doing your Gauss integral you end up asking yourself "what is the value of the gradient at the face" and set that to be mu*(u2-u1)/h. You get to construct one based on the values on an element and it's neighbour
That's a shortcoming of FV in a way because if you want to construct anything of higher order you'll have to querry values form the neighbours of neighbours of neighbours... In 2d/3d FV and with unstructured grids, constructing this isn't easy. 

The second term in the DG weak form: the jump(u) *mu * average(gradv), is supposed to be exactly that, but it gets cancelled out when you use DG0. In DG grad(v) on a face is calculated based on data of the local cell, in FV grad(v) gets calculated based on data from the neighbouring cells. BTW, that term is what you get from integrating by parts. 
Back to your question and if you insist on using DG0:
So to match the penalty*jump(u)*jump(v) with what you get from an FV  scheme (mu*(u2-u1)/h) maybe you might try to set the penalty to be something like (mu/h) where h is the distance between the element centers. (That's not it tho)

 

1) 

That penalty term is there to enhance stability but I don't know if there is more to it. It penalizes the jump across the elements. See if going higher with a constant penalty helps with the accuracy.
Or maybe scale it with the element size. 
There are other DG tutorials for the Poisson equation maybe they'll have something there for you too. 
(Never liked that penalty term tbh either. All my homies hate the penalty parameter. JK). 

Abbas

[Nik Markus Leuenberger]

Dear Prof. Bangerth, Dear Abbas,

Thank you very much again for your fast responses and apologies for my late reply. Further below are some small comments on your answers.
However, rather than trying to make the DG version work in 1D and using 0-order elements, I think I would like to proceed as follows:

• try to implement a discontinuous Galerkin-based version of the equations I'm trying to solve as described in https://www.sciencedirect.com/science/article/pii/S0377042718302875 in sections 1 (equations) and 2.1 (dg discretization of the equations). Only the fine scale equations, not the model reduction equations in section 2.2.
  The goal here would not be to recover piecewise constant finite volume formulations but actually to take advantage of the possibility to use higher order elements. Therefore, I could try to do this in d>=2 space dimensions using elements of order >= 1 to avoid the limitations discussed below. This way, I can take full advantage of the deal.ii capabilities for matrix assembly, mesh refinement etc.
  
  If anyone has the time to take a look at the equations (17)--(20) in the paper and sees something that doesn't seem possible with deal.ii, I would appreciate it. However, the formulation looks pretty close to the one in step-74.

•     if I want a version that really focuses and sticks to the finite volume, piecewise constant discretization, I would try to see, how a framework like OpenFoam could be used for this purpose. Because if I brute force a finite volume discretization in deal.ii, I think I might lose some of the nice features that it offers with FEM approaches like discontinuous Galerkin FEM methods.
  
  However, if you think I'm wrong and a finite volume implementation in deal.ii could be well integrated with most of its tools, I'm very eager to get feedback or hints on how I should approach this.

I know it's a lot of text and some open-ended questions, so please don't feel obligated to answer/comment on all of them.

Thank you very much in advance and best regards,
Nik

On Sunday, June 23, 2024 at 2:10:16 PM UTC+2 abbas.b...@gmail.com wrote:

Nick, 

2)
When you use DG0, all of the terms in the weak from cancel to zero except for the face term penallty*jump(u)*jump(v).  (du/dx and dv/dy are zero inside an element)
That penalty term in the DG community is seen as something that adds stability and isn't at the core of the scheme.  

In FV after doing your Gauss integral you end up asking yourself "what is the value of the gradient at the face" and set that to be mu*(u2-u1)/h. You get to construct one based on the values on an element and it's neighbour
That's a shortcoming of FV in a way because if you want to construct anything of higher order you'll have to querry values form the neighbours of neighbours of neighbours... In 2d/3d FV and with unstructured grids, constructing this isn't easy. 

The second term in the DG weak form: the jump(u) *mu * average(gradv), is supposed to be exactly that, but it gets cancelled out when you use DG0. In DG grad(v) on a face is calculated based on data of the local cell, in FV grad(v) gets calculated based on data from the neighbouring cells. BTW, that term is what you get from integrating by parts. 
Back to your question and if you insist on using DG0:
So to match the penalty*jump(u)*jump(v) with what you get from an FV  scheme (mu*(u2-u1)/h) maybe you might try to set the penalty to be something like (mu/h) where h is the distance between the element centers. (That's not it tho)

Thank you very much Abbas for the elaboration on the comparison between the finite volume and DG approach. This explanation is very helpful in understanding what is going on.

 

1) 

That penalty term is there to enhance stability but I don't know if there is more to it. It penalizes the jump across the elements. See if going higher with a constant penalty helps with the accuracy.
Or maybe scale it with the element size. 
There are other DG tutorials for the Poisson equation maybe they'll have something there for you too. 
(Never liked that penalty term tbh either. All my homies hate the penalty parameter. JK). 

Thank you very much for your inputs. Yes, these are definitely things that I could try if I really wanted to make it work for 1D and order 0.

Yes, it's the same for me, the penalty parameter makes it a bit confusing to understand.

Regarding other tutorials: exactly, 

• step-61 for example uses the so-called Weak-Galerkin method where results for piecewise constant elements are shown.
• step-21 also uses "zero-th degree elements" in the example.

Maybe I can go back to them if I really wanted to focus on getting the piecewise consant, 1D examples to work.

Abbas

On Sunday, June 23, 2024 at 1:24:38 AM UTC+2 Wolfgang Bangerth wrote:

Nik:

> (1) trying to solve in 1D:
>
> * I get errors of the following type:
> *error: ‘class dealii::DoFAccessor<0, 1, 1, false>’ has no member
> named ‘measure’
>   411 | re() / cell->face(face_no)->measure();* / I can make it run
> by removing the calls to these functions and replace the penalty
> factor which would involve them by a constant number, which according
> to the output seems to work o.k. but might need further tuning of the
> penalty constant (here just chosen as a constant of 100)./

Yes, that's because the measure of a face in 1d is hard to define -- what is
the measure of a point in a zero-dimensional space?

If you have a good idea for how the 'extent1/2' variables should be defined in
1d, let us know and we can patch the program so that it will also work in 1d.

Exactly. And no, I don't have a good idea how to define these extend variables to make the program work in 1D in all cases. 

For a fixed grid size, I think a constant penalty factor that would need to be found should work. If I go back and try it,

I will let you know if I have an answer.

I found this paper: https://www.sciencedirect.com/science/article/pii/S0377042714003057, that compares different DG approaches

for the same elliptic equation as in step-74. For some of the approaches, like weak Galerkin finite element methods (WGFEMs) in section 3.1 or mixed finite element methods (MFEMs) in section 3.2, they can use function spaces with order >=0 whereas for the *classical*  discontinuous Galerkin finite element methods (DGFEMs) in section 3.3, they define them using >=1 for the polynomial order. I don't know if that's a proven bound or if it's just a choice to not use the lowest order elements.

[Wolfgang Bangerth]

> * try to implement a discontinuous Galerkin-based version of the equations

> I'm trying to solve as described in
> https://www.sciencedirect.com/science/article/pii/S0377042718302875
> <https://www.sciencedirect.com/science/article/pii/S0377042718302875> in sections 1 (equations) and 2.1 (dg discretization of the equations). Only the fine scale equations, not the model reduction equations in section 2.2.
> The goal here would not be to recover piecewise constant finite volume
> formulations but actually to take advantage of the possibility to use
> higher order elements. Therefore, I could try to do this in d>=2 space
> dimensions using elements of order >= 1 to avoid the limitations discussed
> below. This way, I can take full advantage of the deal.ii capabilities for
> matrix assembly, mesh refinement etc.
>
> If anyone has the time to take a look at the equations (17)--(20) in the
> paper and sees something that doesn't seem possible with deal.ii, I would
> appreciate it. However, the formulation looks pretty close to the one in
> step-74.

If you start from (7), it's really just a heat equation coupled to a Laplace
equation. Looking at Fig. 5, the solution is also smooth. I don't think you
need anything like FV or DG discretizations -- they are good for problems with
non-smooth solutions. If you can separate the two parts of the domain into
individual cells, DG is fine but you might also get away with a continuous
element approach if you put the different solution variables that are
discontinuous along the interface into different vector conditions (-> step-46).


> *     if I want a version that really focuses and sticks to the finite

> volume, piecewise constant discretization, I would try to see, how a
> framework like OpenFoam could be used for this purpose. Because if I brute
> force a finite volume discretization in deal.ii, I think I might lose some
> of the nice features that it offers with FEM approaches like discontinuous
> Galerkin FEM methods.
>
> However, if you think I'm wrong and a finite volume implementation in
> deal.ii could be well integrated with most of its tools, I'm very eager to
> get feedback or hints on how I should approach this.

It's a framing problem. Many people think that FVM is simply a special case of
FEM where you choose discontinuous shape functions and define nonlinear fluxes
as appropriate. I don't think that FVM is conceptually different than FEM.
Whether you call it FVM or DG-FEM, it's two sides of the same coin.

Best
W>