dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

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.

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.
>
> --
> Anders