DOLFIN

calculating partial derivatives

Asked by Souvik Roy

I am trying to solve a vorticity problem. I have got the exact expression for vorticity and I am trying to match the numerical vorticity. My exact velocity is represented as vel and my exact vorticity as vor. The expression for vor and vel are given below:

class exact(Expression):
     def eval(self, value, x):
     value[0] = x[1]**2
            value[1] = x[0]**2
     def value_shape(self):
        return (2,)

vel = exact()

class vortex(Expression):
     def eval(self, value, x):
     value[0] = 2*(x[0]-x[1])

vor = vortex()

The vorticity is obtained from the velocity by setting vor = v_x-u_y where vel=(u,v) and v_x is the partial derivative of v
wrt x and u_y is the partial derivative of u wrt y.

Just to check I put this in my code:

exvor = Function(Z)
exvor = interpolate(vor,Z)

uvor = Function(Z)
uvor = project(D(uob[1],0)-D(uob[0],1),Z)

curlerr = errornorm(exvor, uvor,'L2')

where Z = FunctionSpace(mesh,"CG", 1)

and mesh = UnitSquare(10,10)

In principle I am just using the partial derivative operator on the exact velocity to get the vorticity. So I expect the obtained vorticity to match with the Expression for vorticity given above.

But the curlerr = 0.0197852169412

I changed the number of points in the mesh from 10 to 100 but the curlerr is 0.000620929306036.
I expect it to reach machine precision (10^-15) as it has been achieved for the other quantities in the same code. I am
basically trying to do a patch test. What could be the problem?

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jan Blechta (blechta) said : 		#1

On Fri, 19 Apr 2013 19:06:13 -0000
Souvik Roy <email address hidden> wrote:
> New question #227070 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/227070
>
> I am trying to solve a vorticity problem. I have got the exact
> expression for vorticity and I am trying to match the numerical
> vorticity. My exact velocity is represented as vel and my exact
> vorticity as vor. The expression for vor and vel are given below:
>
> class exact(Expression):
> def eval(self, value, x):
> value[0] = x[1]**2
> value[1] = x[0]**2
> def value_shape(self):
> return (2,)
>
> vel = exact()
>
> class vortex(Expression):
> def eval(self, value, x):
> value[0] = 2*(x[0]-x[1])
>
> vor = vortex()
>
> The vorticity is obtained from the velocity by setting vor = v_x-u_y
> where vel=(u,v) and v_x is the partial derivative of v wrt x and u_y
> is the partial derivative of u wrt y.
>
> Just to check I put this in my code:
>
> exvor = Function(Z)
> exvor = interpolate(vor,Z)
>
> uvor = Function(Z)
> uvor = project(D(uob[1],0)-D(uob[0],1),Z)
>
> curlerr = errornorm(exvor, uvor,'L2')
>
> where Z = FunctionSpace(mesh,"CG", 1)
>
> and mesh = UnitSquare(10,10)
>
> In principle I am just using the partial derivative operator on the
> exact velocity to get the vorticity. So I expect the obtained
> vorticity to match with the Expression for vorticity given above.
>
> But the curlerr = 0.0197852169412
>
> I changed the number of points in the mesh from 10 to 100 but the
> curlerr is 0.000620929306036. I expect it to reach machine precision
> (10^-15) as it has been achieved for the other quantities in the same
> code. I am basically trying to do a patch test. What could be the
> problem?

Matching FE spaces. For example for CG2 velocity its derivatives are
from DG1.

Revision history for this message
Anders Logg (logg) said : 		#2

FEniCS no longer uses Launchapd for Questions & Answers. Please
consult the documentation on the FEniCS web page for where and
how to (re)post your question: http://fenicsproject.org/support/

Revision history for this message
Anders Logg (logg) said : 		#3

FEniCS no longer uses Launchapd for Questions & Answers. Please
consult the documentation on the FEniCS web page for where and
how to (re)post your question: http://fenicsproject.org/support/

Can you help with this problem?

Provide an answer of your own, or ask Souvik Roy for more information if necessary.

To post a message you must log in.