dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Mar 04, 2009 at 09:08:49AM +0100, 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.
> 
> 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.

Would that work? The added term must be linear in a TestFunction, so
it can't be a Function.

-- 
Anders


> 
> Johan

> // Copyright (C) 2008 Johan Hake.
> // Licensed under the GNU LGPL Version 2.1.
> 
> /* Streamline Upwind Petrow/Galerkin stabilizing coefficient function
> 
>   The perturbed testfunction, v', is given by
> 
>      v' = v + s*grad(v) 
> 
>   where the stabilizing term s is given by
>  
>     s = tau *h*tau_l/(2*|a|)*a
> 
>   where a is the field, tau a global tuning parameter, (not inlcuded here),
>    and tau_l the local stabilization parameter
> 
>  */
> 
> class Stab: public Function {
> public:
>   double D; 
>   Function *field;
>   Stab(const FunctionSpace& V):
>     Function(V),D(1.0){}
>   
>   void eval(double* v, const Data& data) const
>   {
>     if (!field)
>        error("Attach a field function.");
>     double field_norm = 0.0;
>     double tau = 0.0;
>     double h = data.cell().diameter();
>     UFCCell ufc_cell(data.cell());
>     field->eval(v, data.x, ufc_cell, data.cell().index());
>     for (uint i = 0;i < geometric_dimension(); ++i)
>       field_norm += v[i]*v[i];
>     field_norm = sqrt(field_norm);
>     double PE = 0.5*field_norm * h/D;
>     if (PE > DOLFIN_EPS)
>       tau = 1/tanh(PE)-1/PE;
>     for (uint i = 0;i < geometric_dimension(); ++i)
>       v[i] *= 0.5*h*tau/field_norm;
>   }
> };

Attachment: signature.asc
Description: Digital signature