dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Wednesday 04 March 2009 18:55:20 Hatef Monajemi wrote:
> Johan Hake wrote:
> > On Wednesday 04 March 2009 08:21:11 Anders Logg wrote:
> > > On Tue, Mar 03, 2009 at 08:08:56PM -0500, Hatef Monajemi wrote:
> > > > Hi
> > > >
> > > > I was just wondering if it is possible to implement the
>
> Petrov-Galerkin
>
> > > > type of weighting functions in dolfin. I am trying to solve a
>
> convection
>
> > > > dominated problem in which the use of standard Galerkin weighting
> > > > generates the oscillations which can be removed by upwinding. Is
>
> there
>
> > > > any demo which shows the implementation of Petrov-Galerkin type of
> > > > testfunctions.
> > > >
> > > > Thanks,
> > > >
> > > > Hatef
> > >
> > > Yes, look at
> > >
> > >   demo/pde/stokes/stabilized
> > >
> > > It's very straightforward, just write
> > >
> > >   v = TestFunction(V)
> > >   v = v + delta*A(v)
> > >
> > > where delta is a stabilization parameter and A is a suitable
>
> operator.
>
> That works fine when you have one Testfunction like "v", How about when
> you have the following set of Testfunction:
> -->
> P2 = VectorElement("Lagrange", "triangle", 1)
> P1 = FiniteElement("Lagrange", "triangle", 1)
> TH = P2 + P1
>
> (v,w) = TestFunctions(TH)
> -->
>
> The analogy with what anders said, one should be able to do something
> like :
>
> (v,w)=( v + delta*A(v), w)
>
> However, this creates problem in ffc. An alternative is to multiply the
> whole equation sets by the added stabilization term "delta*A(v)" and
> add the new terms to the unstabilized variational form. But it is much
> easier if one can just change the weighing functions v to v+delta*A(v).
> is there anyway to do that?
>
> > The stabilizing parameter showed in this demo is quite crude. If it
>
> wont work
>
> > for you, you can consider:
> >
> >   P = FunctionSpace(mesh,"CG",1)
> >   V = VectorFunctionSpace(mesh,"CG"2)
> >   v = TestFunction(P)
> >
> >   field = Function(V,"Some vector valued function that drives the
>
> advection")
>
> >   stab = Function(V,open("SUPG_stab.h").read()) # V is a Vector
> >   stab.field = field
> >
> >
> >   v = v + 0.2*dot(stab,grad(v))
> >
> > Here 0.2*stab relates to the delta Anders mention and A relates to grad
> >
> > We could consider adding the handed function as a SpecialFunction in
>
> DOLFIN.
>
> > It is a general SUPG stabilization term, which depends not only on the
>
> local
>
> > meshsize but also on the local field size.
> >
> >
> > Johan
>
> Dear Johan:
>
> Thanks for the better wight functions. I think these weight functions
> are those first introduced by Hughes and Brooks and also Kelly. 

It was at least from Brooks and Hughs paper from 1982 I got the term from.

> Just one question? where in the element do you calculate your field
> variable? 

Not sure I get what you aiming at, but I call eval of the attached field 
function. This is done in the eval function of the stabilizing term.

> this method somehow assumes that field variable (velocity in 
> my case) is constant over an element.

This method does not require a constant field, that's the whole point, but if 
you have a constant field you can just define a constant Function:

  field = Constant(mesh,(1,2,3))

and attach this function to the stabilizing function.

Johan

> Hatef
>
>
>
> ------------------------------------------------------------------------
>                    Message Attached :SUPG_stab.h
> ------------------------------------------------------------------------
>                 Message Attached :attached
>
> --------------------------------------
> Hatef Monajemi
> Graduate Student
> Dept of Civil and Environmental Eng
> 2032 MC Building, Carleton University
> 1125 Colonel By Drive, Ottawa
> Ontario, Canada, K1S 5B6