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 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. > > > > (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. > > 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/6f34f4a1-c254-4f49-969c-08513c7b5d6cn%40googlegroups.com.
