dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jun 01, 2010 at 01:42:54PM +0100, Garth N. Wells wrote:
>
>
> On 01/06/10 13:39, Anders Logg wrote:
> >On Thu, May 27, 2010 at 04:37:42PM +0100, Garth N. Wells wrote:
> >>
> >>
> >>On 27/05/10 16:34, Anders Logg wrote:
> >>>On Thu, May 27, 2010 at 04:22:33PM +0100, Garth N. Wells wrote:
> >>>>
> >>>>
> >>>>On 27/05/10 15:58, Anders Logg wrote:
> >>>>>Wouldn't it work to just to this:
> >>>>>
> >>>>>v = v + (h/2.0*vnorm)*dot(velocity, grad(v))
> >>>>>U = 0.5*(u0 + u)
> >>>>>F = v*(u - u0)*dx + 0.5*k*(v*dot(velocity, grad(U))*dx + c*dot(grad(v), grad(U))*dx) - k*v*f*dx
> >>>>>a = lhs(F)
> >>>>>L = rhs(F)
> >>>>>
> >>>>>That's my understanding of SUPG, but perhaps I've missed something?
> >>>>>
> >>>>
> >>>>Looks like the same thing (I'd have to check). I prefer to
> >>>>expressing the method as weighting the residual, e.g. Galerkin +
> >>>>weighted residual.
> >>>
> >>>I also think of it as weighting the residual, but with a modified test
> >>>function v. The above looks like a weighted (weak) residual to me.
> >>>
> >>>The difference from your implementation is that you multiply the
> >>>stabilization part of the test function with the strong residual for
> >>>which the Laplacian (for P1 elements) is zero. In the above case, the
> >>>stabilization also hits the nonzero integrated by parts grad-grad term.
> >>>
> >>>I don't know whether or not that is important.
> >>>
> >>>>The code I added can be made more compact using lhs/rhs, but I was
> >>>>having some trouble and initially wanted to make a bare bones
> >>>>implementation.
> >>>
> >>>ok.
> >>>
> >>>I might try to make some changes. I'd like to get a feeling for how
> >>>SUPG works.
> >>>
> >>
> >>OK, but just give me moment since I'm making a change to use lhs/rhs.
> >>
> >>Garth
> >
> >On second thought, even if I like to think of it as modifying the test
> >function v -->  v + delta * beta . grad(v), this is problematic when
> >delta * beta . grad(v) hits div(grad(u)) since we run into
> >difficulties of either defining div(grad(u)) for an H^1 function or
> >integrating by parts onto v which is already differentiated.
> >
>
> I don't see the issue - we control H^1, not H^2, so the div(grad(u))
> term is fine.
>
> The div(grad(u)) is essential for P2 elements to maintain consistency.

My point is that div(grad(u)) is not well-defined globally for H^1
functions (in particular P1 and P2) so when you write

  r = u-u0 + dt*( dot(velocity, grad(u_mid)) -div(grad(u_mid)) - f)
  F += (h/2.0*vnorm)*dot(velocity, grad(v))*r*dx

that should be interpreted locally, as in computing the residual on
each cell, integrating over each cell and summing the results. This is
also what happens when DOLFIN assembles so everything is fine.

The div(grad(u)) term will be zero (element-wise, but globally a delta
function over edges) for P1 but nonzero for P2.

I'm just saying that the demo is ok and I think I understand why now. :-)

--
Anders



> Garth
>
> >So I agree the best formulation is this:
> >
> >   a(v, u_h) + sum_T<delta * beta . grad(v), R_T>_T = L(v)
> >
> >where R_T is the element-wise residual (where we may differentiate u_h
> >twice and get zero). R_T will be zero for u so the method is consistent.
> >
> >This is what is currently implemented in the advection-diffusion demo.
> >
>

Attachment: signature.asc
Description: Digital signature