dolfin team mailing list archive

[Hatef Monajemi <hmonajem@xxxxxxxxxxxxxxxxxxx>]

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.

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