SUPG implementation

Asked by Marcos Samudio on 2012-02-25

Hello everyone,

I am trying to add the SUPG terms to a steady NS solver form (PSPG already working, from CBC.PDESys package).

The input objects are

u_= Previous step velocity
u = Unknown actual velocity
p = Unknown actual pressure
v = Velocity Test Function
esu = SUPG stabilization parameter (a constant)

So I define my SUPG stabilization like this

def supg_stabilization(self, u_, u, p, v, esu, n):
        """Add SUPG stabilization to Navier Stokes solver. Bilinear shapefunctions"""
        esu = esu*u_
        F = - esu*inner(grad(q), grad(p))*dx \
              - esu*inner(grad(v), dot(grad(u), u_))*dx
        for bc in self.bcs:
            if bc.type() in ('ConstantPressure', 'Outlet',
                             'VelocityInlet'):
                F = F + esu*inner(v, dot(grad(p), n))*ds(bc.bid)
        return F

I am having trouble with the inner products. I know tensor dimensions are not coherent, but I do not know how to fix it.

Any ideas will be appreciated!

Thanks,

Marcos

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Marcos Samudio
Solved:
Last query:
Last reply:

Hi Marcos,

You can ask questions specific to CBC.PDESys through their own launchpad
website.

On Saturday, 25 February 2012, Marcos Samudio wrote:

> New question #188834 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/188834
>
> Hello everyone,
>
> I am trying to add the SUPG terms to a steady NS solver form (PSPG
> already working, from CBC.PDESys package).
>
> The input objects are
>
> u_= Previous step velocity
> u = Unknown actual velocity
> p = Unknown actual pressure
> v = Velocity Test Function
> esu = SUPG stabilization parameter (a constant)
>
> So I define my SUPG stabilization like this
>
> def supg_stabilization(self, u_, u, p, v, esu, n):
> """Add SUPG stabilization to Navier Stokes solver. Bilinear
> shapefunctions"""
> esu = esu*u_
> F = - esu*inner(grad(q), grad(p))*dx \
> - esu*inner(grad(v), dot(grad(u), u_))*dx
> for bc in self.bcs:
> if bc.type() in ('ConstantPressure', 'Outlet',
> 'VelocityInlet'):
> F = F + esu*inner(v, dot(grad(p), n))*ds(bc.bid)
> return F
>
> I am having trouble with the inner products. I know tensor dimensions are
> not coherent, but I do not know how to fix it.

Did you have a look at the stabilization routine in NSCoupled? It's
similar. Looks like the ranks of things are wrong here. If v is vector and
p scalar then inner(v, dot(grad(p), n)) is trying to take the inner product
of a vector with a scalar. Likewise, inner(grad(v), dot(grad(u), u_)) tries
to multiply a second rank tensor (grad(v)) with a vector. Not knowing all
that you're trying to do I think perhaps it should be grad(q) here first.

Mikael

>
> Any ideas will be appreciated!
>
> Thanks,
>
> Marcos
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Marcos Samudio (marcossamudio) said : #2

Hi Mikael,

Yes, I am having trouble with the ranks. I just can not figure out how to make it work in Dolfin. I have to calculate the integral u.grad(v)(u.grad(u) + grad(p)) in order to stabilize the NS equations. I'm unable to express this correctly in Dolfin's syntax.

Any suggestions? I am searching the demos for something similar, but I can't find this type of formulation yet.

Thanks a lot,

Marcos

Marcos Samudio (marcossamudio) said : #3

I think this is the right way to define F:

stab = dot(grad(v), u_)
F = - esu*inner(stab, grad(p))*dx - esu*inner(stab, dot(grad(u), u_))*dx

And so now I have to fix it for the natural BCs. Is there a way to do it properly, the way I posted it at the beginning of the thread?

Marcos

Den Feb 25, 2012 kl. 9:15 PM skrev Marcos Samudio:

> Question #188834 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/188834
>
> Marcos Samudio posted a new comment:
> I think this is the right way to define F:
>
> stab = dot(grad(v), u_)
> F = - esu*inner(stab, grad(p))*dx - esu*inner(stab, dot(grad(u), u_))*dx
>

That should work.

> And so now I have to fix it for the natural BCs. Is there a way to do it
> properly, the way I posted it at the beginning of the thread?

I may be wrong, but I donĀ“t think there should be a natural BC for this form. You have simply multiplied the steady momentum equation (neglecting the laplacian) with the test function dot(grad(v), u_) right? So I guess there's no integration by parts and no ds term.

Mikael

>
> Marcos
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Marcos Samudio (marcossamudio) said : #5

Yes, I think you are right. Indeed it is only the steady momentum equation, and no laplacian since I'm only using linear elements.

Thanks Mikael