# SUPG implementation

Hello everyone,

I have already posted this question in the Dolfin forum, but I guess it is rather pertinent to this forum.

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_stabilizat

"""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 ('ConstantPress

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:
- CBC.PDESys Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Marcos Samudio

- Solved:

- Last query:

- Last reply:

Hi

I just thought you should know the Laplacian is not zero just because you have linear elements. This is a common mistake and one of the main reasons why the G2 method obtains really bad results in the Naver-Stokes benchmark in the Fenics book. Just integrate by parts and you get a inner(grad, grad) term plus a ds term, just like any other Laplacian.

Best regatds

Mikael

Marcos Samudio (marcossamudio) said : | #2 |

Thanks a lot Mikael, I will take a look at that chapter on the Fenics book.

I am having another issue when trying to solve N-S with cbc.pdesys. Everything goes normal, convergence is relatively fast and smooth, with not large oscillations, and all of sudden I get a 'nan' (not a number) error. There are no signs of divergence or anyting similar. Has anyone else had this same type of problem? (should I start a new question?)

Once again, thanks for your help!

Marcos?

On Monday, 5 March 2012, Marcos Samudio <

<email address hidden>> wrote:

> Question #188855 on CBC.PDESys changed:

> https:/

>

> Marcos Samudio posted a new comment:

> Thanks a lot Mikael, I will take a look at that chapter on the Fenics

> book.

>

> I am having another issue when trying to solve N-S with cbc.pdesys.

> Everything goes normal, convergence is relatively fast and smooth, with

> not large oscillations, and all of sudden I get a 'nan' (not a number)

> error. There are no signs of divergence or anyting similar. Has anyone

> else had this same type of problem? (should I start a new question?)

Actually, i have had a similar problem with one of the RANS models and I

was not able to figure out what caused it. Just nans all of a sudden when

everything seemed to be ok. Are you running a NS solver or RANS? Newton? Is

it reproducible? I would be very interested in figuring this one out...

Mikael

>

> Once again, thanks for your help!

>

> Marcos?

>

> --

> You received this question notification because you are an answer

> contact for CBC.PDESys.

>

Marcos Samudio (marcossamudio) said : | #4 |

I am running a NS solver, with Picard iterations. The 'nan' error shows up in the same iteration every time I run the code. I have tried modifying parameters, refining the mesh and gradual start-up (slowly developing BCs over time) for the simple lid-driven cavity test, and this keeps on happening.

I have modified some of your code, and my stabilization functions look like this:

def pspg_stabilizat

"""Add stabilization to Navier Stokes solver"""

h = CellSize(mesh)

eps = eps*h

if type(nu) is Constant:

F = - eps*inner(grad(q), grad(p))*dx - eps*inner(grad(q), dot(grad(u), u_))*dx

if self.prm[

u_1 = kwargs.get('u_1', None)

F = F - (1/self.

# for bc in self.bcs:

# if bc.type() in ('ConstantPress

# 'VelocityInlet'):

# F = F + eps*inner(q, dot(grad(p), n))*ds(bc.bid)

for bc in self.bcs:

if bc.type() in ('Wall', 'Outlet', 'VelocityInlet'):

else:

F = - eps*inner(grad(q), grad(p))*dx \

- eps*inner(grad(q), dot(grad(u), u_))*dx \

- eps*inner(

+ eps*inner(

for bc in self.bcs:

if bc.type() in ('ConstantPress

return F

#esu is the supg stabilization parameter

def supg_stabilizat

"""Add SUPG stabilization to Navier Stokes solver. v is the Test Function

for velocity"""

h = CellSize(mesh)

esu = esu*h

stab = dot(grad(v), u_)

if type(nu) is Constant:

F = esu*inner(stab, grad(p))*dx \

+ esu*inner(stab, dot(grad(u), u_))*dx \

+ esu*nu*

if self.prm[

u_1 = kwargs.get('u_1', None)

F = F + 1/self.

for bc in self.bcs:

if bc.type() in ('Wall', 'Outlet', 'VelocityInlet'):

##FIXME: NO ESTA IMPLEMENTADO SUPG PARA NU VARIABLE

else:

F = - eps*inner(grad(q), grad(p))*dx \

- eps*inner(grad(q), dot(grad(u), u_))*dx \

- eps*inner(

+ eps*inner(

for bc in self.bcs:

if bc.type() in ('ConstantPress

return F

Marcos Samudio (marcossamudio) said : | #5 |

Sorry for that previous message, I accidentally hit the tab key and it just got published.(now I can not erase it...)

What I meant is:

I am running a NS solver, with Picard iterations. The 'nan' error shows up in the same iteration every time I run the code. I have tried modifying parameters, refining the mesh and adding gradual start-up (slowly developing BC over time) for the simple lid-driven cavity test, and this keeps on happening; these changes only modify the moment at which the error occurs.

I have modified some of your code, and my stabilization functions look like this:

def pspg_stabilizat

"""Add stabilization to Navier Stokes solver"""

h = CellSize(mesh)

eps = eps*h

if type(nu) is Constant:

F = - eps*inner(grad(q), grad(p))*dx - eps*inner(grad(q), dot(grad(u), u_))*dx

if self.prm[

u_1 = kwargs.get('u_1', None)

F = F - (1/self.

for bc in self.bcs:

if bc.type() in ('Wall', 'Outlet', 'VelocityInlet'):

return F

def supg_stabilizat

"""Add SUPG stabilization to Navier Stokes solver."""

h = CellSize(mesh)

esu = esu*h

stab = dot(grad(v), u_)

if type(nu) is Constant:

F = esu*inner(stab, grad(p))*dx + esu*inner(stab, dot(grad(u), u_))*dx \

+ esu*nu*

if self.prm[

u_1 = kwargs.get('u_1', None)

F = F + 1/self.

for bc in self.bcs:

if bc.type() in ('Wall', 'Outlet', 'VelocityInlet'):

return F

Those are the main changes I have made, and I am running it with stabilization parameters esu and eps of around 0.5.

If you see anything wrong or unexpected, please let me know. I will try to generate the simplest example that reproduces the error.

Thanks,

Marcos

I canĀ“t see anything obviously wrong. Perhaps you should register your branch on the launchpad site so I can test your case that goes nan?

Mikael

Den Mar 5, 2012 kl. 10:40 PM skrev Marcos Samudio:

> Question #188855 on CBC.PDESys changed:

> https:/

>

> Marcos Samudio posted a new comment:

> Sorry for that previous message, I accidentally hit the tab key and it just got published.(now I can not erase it...)

> What I meant is:

>

> I am running a NS solver, with Picard iterations. The 'nan' error shows

> up in the same iteration every time I run the code. I have tried

> modifying parameters, refining the mesh and adding gradual start-up

> (slowly developing BC over time) for the simple lid-driven cavity test,

> and this keeps on happening; these changes only modify the moment at

> which the error occurs.

>

> I have modified some of your code, and my stabilization functions look

> like this:

>

> def pspg_stabilizat

> """Add stabilization to Navier Stokes solver"""

> h = CellSize(mesh)

> eps = eps*h

> if type(nu) is Constant:

> F = - eps*inner(grad(q), grad(p))*dx - eps*inner(grad(q), dot(grad(u), u_))*dx

> if self.prm[

> u_1 = kwargs.get('u_1', None)

> F = F - (1/self.

> for bc in self.bcs:

> if bc.type() in ('Wall', 'Outlet', 'VelocityInlet'):

> F = F + eps*nu*

> return F

>

> def supg_stabilizat

> **kwargs):

> """Add SUPG stabilization to Navier Stokes solver."""

> h = CellSize(mesh)

> esu = esu*h

> stab = dot(grad(v), u_)

> if type(nu) is Constant:

> info_green('Adding SUPG stabilization for constant nu')

> F = esu*inner(stab, grad(p))*dx + esu*inner(stab, dot(grad(u), u_))*dx \

> + esu*nu*

> if self.prm[

> u_1 = kwargs.get('u_1', None)

> F = F + 1/self.

> for bc in self.bcs:

> if bc.type() in ('Wall', 'Outlet', 'VelocityInlet'):

> F = F - esu*nu*inner(stab, dot(grad(u), n))*ds(bc.bid)

> return F

>

> Those are the main changes I have made, and I am running it with

> stabilization parameters esu and eps of around 0.5.

>

> If you see anything wrong or unexpected, please let me know. I will try

> to generate the simplest example that reproduces the error.

>

> Thanks,

>

> Marcos

>

> --

> You received this question notification because you are an answer

> contact for CBC.PDESys.