dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Wednesday 04 March 2009 20:28:36 Hatef Monajemi wrote:
> Johan Hake wrote:
> > 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,
>
> Well, may be I did not mention my point in a lucid way. What I am
> trying to ask is :
>
> for the following stabilization term :
>
> 0.5*alpha*h {u(i)/|u|} ( partial N_k / partial x_i)
>
> and Peclet number,
>
> Pe =0.5*|u| *h /k
>
> in 2D problem we will have:
> |u|=u_x^2+u_y2
>
> 0.5*alpha*h (     {u_x/|u|} ( partial N_k / partial x) + u_y/|u|} (
> partial N_k / partial y)
>
> So given the above formulas, and refer to Zienkewiecz page.455 vol2 ,
> "The method presuppose that the velocity components (i.e. u_x and u_y)
> in a PARTICULAR ELEMENT are substantially constant".
>
> So, I do not mean that the method assumes that velocity field is
> constant (of course it is not) , but when you are going  to calculate
> the stabilization term you should come up with  u_x and u_y for the
> element (cell) in two directions. The problem will be answered if I
> find out what does "data.x" do in eval function "eval(v, data.x,
> ufc_cell, data.cell().index())" because I think that is where you are
> calculating your solution in eval.

There I call the velocity function, with the coordinate given in data.x. The 
ufc_cell, data.cell().index() are just suger, making the call faster. 

Johan

> Thanks again.
> Hatef
>
> > 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
>
> --------------------------------------
> Hatef Monajemi
> Graduate Student
> Dept of Civil and Environmental Eng
> 2032 MC Building, Carleton University
> 1125 Colonel By Drive, Ottawa
> Ontario, Canada, K1S 5B6