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. > (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. 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. Best W. -- ------------------------------------------------------------------------ Wolfgang Bangerth email: bang...@colostate.edu www: http://www.math.colostate.edu/~bangerth/ -- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/fa054081-8ea7-4bfa-884b-5a0865ecabc6n%40googlegroups.com.
