How to implement a PDE with curl operator


471
views
-1
8 months ago by
Diego  
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}$nablaxE=dBdt

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
8 months ago by
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
8 months ago by
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, 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
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
8 months ago by
Diego  
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):
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
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 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)·ey   returns the y component of curl(E).
written 8 months ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »