### How to implement a PDE with curl operator

471

views

-1

Hello everybody,

I am trying to solve the simple PDE where the curl of the electric field equates the time derivative of the magnetic field:

$nablaxE=\frac{dB}{dt}$

I get stuck when I need to implement it in FEniCS.

In my script I have written:

F = curl(E)*v*dx - sigma*Bx*v*dx

a,L = lhs(F),rhs(F)

where Bx is the known expression of the time derivative of the magnetic field.

But I get this error:

Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().

Traceback (most recent call last):

File "/usr/lib/python2.7/dist-packages/IPython/core/interactiveshell.py", line 2883, in run_code

exec(code_obj, self.user_global_ns, self.user_ns)

File "<ipython-input-65-3832f10cd577>", line 1, in <module>

F = curl(Jed)*v*dx - sigma*Bx*v*dx

File "/usr/lib/python2.7/dist-packages/ufl/measure.py", line 429, in __rmul__

(integrand.ufl_shape, integrand.ufl_free_indices))

File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 172, in error

raise self._exception_type(self._format_raw(*message))

UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().

How can I solve it?

Thank you for your help.

Bye

Diego

I am trying to solve the simple PDE where the curl of the electric field equates the time derivative of the magnetic field:

$nablaxE=\frac{dB}{dt}$

`n``a``b``l``a``x``E`=`d``B``d``t`I get stuck when I need to implement it in FEniCS.

In my script I have written:

F = curl(E)*v*dx - sigma*Bx*v*dx

a,L = lhs(F),rhs(F)

where Bx is the known expression of the time derivative of the magnetic field.

But I get this error:

Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().

Traceback (most recent call last):

File "/usr/lib/python2.7/dist-packages/IPython/core/interactiveshell.py", line 2883, in run_code

exec(code_obj, self.user_global_ns, self.user_ns)

File "<ipython-input-65-3832f10cd577>", line 1, in <module>

F = curl(Jed)*v*dx - sigma*Bx*v*dx

File "/usr/lib/python2.7/dist-packages/ufl/measure.py", line 429, in __rmul__

(integrand.ufl_shape, integrand.ufl_free_indices))

File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 172, in error

raise self._exception_type(self._format_raw(*message))

UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().

How can I solve it?

Thank you for your help.

Bye

Diego

Community: FEniCS Project

### 3 Answers

4

curl(E) is a vector and presumably the test function v is vector-value as well, so you need

`inner(curl(E), v)*dx`

not `curl(E)*v*dx`

, and similarly for the second term.
sorry for the unrelated question, but how did you get that ttfont? I haven't been able to figure that out...

written
8 months ago by
pf4d

2

Select the text to mark as a code sample and then hit the '{:}' button: "Insert/Edit code sample".

written
8 months ago by
Douglas N Arnold

2

Your E and v are scalar-valued functions while curl(E) is a vector-valued function (in 2D, curl(E) = (d E / d y, -d E / d x) ). Your first attempt,

The PDE problem to find a scalar function E such that curl(E) equals a given scalar function B does not make sense, since curl(E) is a vector. The problem to find a vector-valued function E such that curl(E) equals a given scalar function B is not well-posed because there are many solutions E for any B, no matter what the boundary conditions. I would say that you need to formulate a mathematically meaningful problem before getting to its implementation.

`curl(E)*v*dx`

, fails because you can only integrate a scalar expression and curl(E)*v is vector-valued. Your second attempt, `inner(curl(E), v)*dx`

, fails, because you cannot take the inner product of a vector and a scalar.The PDE problem to find a scalar function E such that curl(E) equals a given scalar function B does not make sense, since curl(E) is a vector. The problem to find a vector-valued function E such that curl(E) equals a given scalar function B is not well-posed because there are many solutions E for any B, no matter what the boundary conditions. I would say that you need to formulate a mathematically meaningful problem before getting to its implementation.

Actually function E is a vector field E = (Ex, Ey, 0) and the curl(E) returns a component in the z-axis (out of plane)

curl(E) = (0, 0, dEx/dy-dEy/dx)

The function B is aligned along the z-axis, it is also a vector field but it has no component in x and y axes.

B = (0, 0, 10).

So the PDE equates the z-component of the curl(E) and z-component of B.

I see from your reply that probably I need to define the PDE problem to the solver in a different way when my unknown is a vector field and not a scalar function?

Thank you

curl(E) = (0, 0, dEx/dy-dEy/dx)

The function B is aligned along the z-axis, it is also a vector field but it has no component in x and y axes.

B = (0, 0, 10).

So the PDE equates the z-component of the curl(E) and z-component of B.

I see from your reply that probably I need to define the PDE problem to the solver in a different way when my unknown is a vector field and not a scalar function?

Thank you

written
8 months ago by
Diego

From the lines

```
V = FunctionSpace(mesh,'P',2)
...
E = TrialFunction(V)
```

your function E is a scalar valued function.
written
8 months ago by
klunkean

0

Thank you Douglas for your reply.

Actually I still get an error on shape matching when I use "inner" instead of "dot". This is the error:

F = inner(curl(E),v)*dx - inner(Bx,v)*dx

File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 145, in inner

return Inner(a, b)

File "/usr/lib/python2.7/dist-packages/ufl/tensoralgebra.py", line 165, in __new__

error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)

Below is my script. The domain is 2D.

# Problem data

W = 0.03 # m

H = 0.015 # m

w_cornice = 0.005 #m

h_cornice = 0.005 #m

w_plate = 0.02 #m

h_plate = 0.005 #m

freq = 700

T = 1/freq

num_steps = 12

dt = T/num_steps # s

k_max = 18

N_elem = 100

sigma_o = 1

sigma_FZ = 59523 # X52 steel liquid state

# ---------- MESH GENERATION -------------------

domain = Rectangle(Point(0,0),Point(W,H)) # subdomain 0

plate_dom = Rectangle(Point(w_cornice,h_cornice),Point(w_cornice+w_plate,h_cornice+h_plate))

domain.set_subdomain(1,plate_dom) # subdomain 1

mesh = generate_mesh(domain, N_elem)

coor = mesh.coordinates()

V = FunctionSpace(mesh,'P',2)

plot(mesh)

markers = MeshFunction('size_t',mesh,2,mesh.domains())

dx = Measure('dx',domain=mesh,subdomain_data=markers)

# ---------------- BOUNDARY CONDITIONS ----------------------

A_D = Expression('0',degree=1,t=0)

def boundary(x, on_boundary):

# ---------------- EXTERNAL LOADS ---------------------------

Bx = Constant(10)

Bx_int = interpolate(Bx,V)

plot(Bx_int)

interactive()

# -------------- ELECTRICAL CONDUCTIVITY --------------------

class conductivity(Expression):

def __init__(self,mesh,**kwargs):

self.markers = markers

def eval_cell(self,values,x,cell):

if markers[cell.index] == 0: # air

values[0] = sigma_o # vacuum

elif markers[cell.index] == 1: # ferromagnetic plate

values[0] = sigma_FZ

sigma = conductivity(mesh,degree=1)

#--------------- VARIATIONAL PROBLEM ------------------------

E = TrialFunction(V)

v = TestFunction(V)

F = inner(curl(E),v)*dx - sigma*inner(Bx,v)*dx

a,L = lhs(F),rhs(F)

E = Function(V)

t = 0

for n in range(num_steps):

t += dt

interactive
Your formulation still does not make sense. The inner product is defined for two arguments coming from the same hilbert space. I just repeat what Douglas already told you but if you write

Either you manually write down the desired component of your curl or you project onto the desired direction using a dot product with a unit base vector:

Actually I still get an error on shape matching when I use "inner" instead of "dot". This is the error:

F = inner(curl(E),v)*dx - inner(Bx,v)*dx

File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 145, in inner

return Inner(a, b)

File "/usr/lib/python2.7/dist-packages/ufl/tensoralgebra.py", line 165, in __new__

error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)

Below is my script. The domain is 2D.

# Problem data

W = 0.03 # m

H = 0.015 # m

w_cornice = 0.005 #m

h_cornice = 0.005 #m

w_plate = 0.02 #m

h_plate = 0.005 #m

freq = 700

T = 1/freq

num_steps = 12

dt = T/num_steps # s

k_max = 18

N_elem = 100

sigma_o = 1

sigma_FZ = 59523 # X52 steel liquid state

# ---------- MESH GENERATION -------------------

domain = Rectangle(Point(0,0),Point(W,H)) # subdomain 0

plate_dom = Rectangle(Point(w_cornice,h_cornice),Point(w_cornice+w_plate,h_cornice+h_plate))

domain.set_subdomain(1,plate_dom) # subdomain 1

mesh = generate_mesh(domain, N_elem)

coor = mesh.coordinates()

V = FunctionSpace(mesh,'P',2)

plot(mesh)

markers = MeshFunction('size_t',mesh,2,mesh.domains())

dx = Measure('dx',domain=mesh,subdomain_data=markers)

# ---------------- BOUNDARY CONDITIONS ----------------------

A_D = Expression('0',degree=1,t=0)

def boundary(x, on_boundary):

return on_boundary

bc = DirichletBC(V,Constant(0),boundary)# ---------------- EXTERNAL LOADS ---------------------------

Bx = Constant(10)

Bx_int = interpolate(Bx,V)

plot(Bx_int)

interactive()

# -------------- ELECTRICAL CONDUCTIVITY --------------------

class conductivity(Expression):

def __init__(self,mesh,**kwargs):

self.markers = markers

def eval_cell(self,values,x,cell):

if markers[cell.index] == 0: # air

values[0] = sigma_o # vacuum

elif markers[cell.index] == 1: # ferromagnetic plate

values[0] = sigma_FZ

sigma = conductivity(mesh,degree=1)

#--------------- VARIATIONAL PROBLEM ------------------------

E = TrialFunction(V)

v = TestFunction(V)

F = inner(curl(E),v)*dx - sigma*inner(Bx,v)*dx

a,L = lhs(F),rhs(F)

E = Function(V)

t = 0

for n in range(num_steps):

t += dt

A_D.t = t

Bx.t = t

solve(a == L, E, bc)

Jed_x, Jed_y = E.split()

plot(E, title='solution at time t = %.4f s' % (t))

interactive

`inner(curl(E),v) `

you basically want to calculate $\text{curl}\left(E\right)\cdot v$curl(`E`)·

`v`where the first argument is a vector and the second a scalar. That operation is not defined.

Either you manually write down the desired component of your curl or you project onto the desired direction using a dot product with a unit base vector:

```
mesh = UnitSquareMesh(3,3)
V = FunctionSpace(mesh,'Lagrange',1)
E = Function(V)
v = TestFunction(V)
ey = Constant((0., 1.))
F = inner(dot(curl(E),ey), v)*dx
```

The operation $\text{curl}\left(E\right)\cdot e_y$curl(`E`)·

`e`

_{y}returns the y component of curl(E).

written
8 months ago by
klunkean

Please login to add an answer/comment or follow this question.