CBC.PDESys

SUPG implementation

Asked by Marcos Samudio

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_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:
CBC.PDESys Edit question
Assignee:
No assignee Edit question
Solved by:
Marcos Samudio
Solved:
Last query:
Last reply:
Revision history for this message
Mikael Mortensen (mikael-mortensen) said : 		#1

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

Revision history for this message
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?

Revision history for this message
Mikael Mortensen (mikael-mortensen) said : 		#3

On Monday, 5 March 2012, Marcos Samudio <
<email address hidden>> wrote:
> Question #188855 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/188855
>
> 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.
>

Revision history for this message
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_stabilization(self, u_, u, mesh, p, q, nu, nut_, eps, n, **kwargs):
        """Add stabilization to Navier Stokes solver"""
        h = CellSize(mesh)
        eps = eps*h
        if type(nu) is Constant:
            info_green('Adding PSPG stabilization for constant nu')
            F = - eps*inner(grad(q), grad(p))*dx - eps*inner(grad(q), dot(grad(u), u_))*dx
            if self.prm['time_integration'] == 'Transient':
                u_1 = kwargs.get('u_1', None)
                F = F - (1/self.prm['dt'])*eps*inner(u - u_1, grad(q))*dx
# for bc in self.bcs:
# if bc.type() in ('ConstantPressure', 'Outlet',
# 'VelocityInlet'):
# F = F + eps*inner(q, dot(grad(p), n))*ds(bc.bid)
            #MODIFICACION DE LAS INTEGRALES DE CONTORNO, segun calculos MS
            for bc in self.bcs:
                if bc.type() in ('Wall', 'Outlet', 'VelocityInlet'):
                    F = F + eps*nu*inner(grad(q), dot(grad(u), n))*ds(bc.bid)
        else:
            info_green('Adding PSPG stabilization for variable nu')
            F = - eps*inner(grad(q), grad(p))*dx \
                - eps*inner(grad(q), dot(grad(u), u_))*dx \
                - eps*inner(grad(grad(q)*nut_), grad(u))*dx \
                + eps*inner(outer(grad(q), grad(nut_)), grad(u) + grad(u).T)*dx
            for bc in self.bcs:
                if bc.type() in ('ConstantPressure', 'Outlet',
                                 'VelocityInlet'):
                    F = F + eps*inner(q, dot(grad(p), n))*ds(bc.bid)
        return F

    #esu is the supg stabilization parameter
    def supg_stabilization(self, u_, u, mesh, p, v, nu, nut_, esu, n, convection_form,
                      **kwargs):
        """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:
            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*inner(grad(stab), grad(u))*dx
            if self.prm['time_integration'] == 'Transient':
                u_1 = kwargs.get('u_1', None)
                F = F + 1/self.prm['dt']*esu*inner(u - u_1, stab)*dx
            #MODIFICACION DE LAS INTEGRALES DE CONTORNO, segun calculos MS
            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)
        ##FIXME: NO ESTA IMPLEMENTADO SUPG PARA NU VARIABLE
        else:
            info_green('Adding SUPG stabilization for variable nu')
            F = - eps*inner(grad(q), grad(p))*dx \
                - eps*inner(grad(q), dot(grad(u), u_))*dx \
                - eps*inner(grad(grad(q)*nut_), grad(u))*dx \
                + eps*inner(outer(grad(q), grad(nut_)), grad(u) + grad(u).T)*dx
            for bc in self.bcs:
                if bc.type() in ('ConstantPressure', 'Outlet',
                                 'VelocityInlet'):
                    F = F + eps*inner(q, dot(grad(p), n))*ds(bc.bid)
        return F

Revision history for this message
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_stabilization(self, u_, u, mesh, p, q, nu, nut_, eps, n, **kwargs):
        """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['time_integration'] == 'Transient':
                u_1 = kwargs.get('u_1', None)
                F = F - (1/self.prm['dt'])*eps*inner(u - u_1, grad(q))*dx
         for bc in self.bcs:
               if bc.type() in ('Wall', 'Outlet', 'VelocityInlet'):
                    F = F + eps*nu*inner(grad(q), dot(grad(u), n))*ds(bc.bid)
        return F

    def supg_stabilization(self, u_, u, mesh, p, v, nu, nut_, esu, n, convection_form,
                      **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*inner(grad(stab), grad(u))*dx
            if self.prm['time_integration'] == 'Transient':
                u_1 = kwargs.get('u_1', None)
                F = F + 1/self.prm['dt']*esu*inner(u - u_1, stab)*dx
            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

Revision history for this message
Mikael Mortensen (mikael-mortensen) said : 		#6

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://answers.launchpad.net/cbcpdesys/+question/188855
>
> 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_stabilization(self, u_, u, mesh, p, q, nu, nut_, eps, n, **kwargs):
> """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['time_integration'] == 'Transient':
> u_1 = kwargs.get('u_1', None)
> F = F - (1/self.prm['dt'])*eps*inner(u - u_1, grad(q))*dx
> for bc in self.bcs:
> if bc.type() in ('Wall', 'Outlet', 'VelocityInlet'):
> F = F + eps*nu*inner(grad(q), dot(grad(u), n))*ds(bc.bid)
> return F
>
> def supg_stabilization(self, u_, u, mesh, p, v, nu, nut_, esu, n, convection_form,
> **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*inner(grad(stab), grad(u))*dx
> if self.prm['time_integration'] == 'Transient':
> u_1 = kwargs.get('u_1', None)
> F = F + 1/self.prm['dt']*esu*inner(u - u_1, stab)*dx
> 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.