dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

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.


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;
  }
};