Hello, I am simulating a system on a box mesh where my vector valued function, which initially has a component only in the x direction should oscillate between x and z without decaying - based on the below strong form, with A=0 (the relaxation term), with initial value on the left face of the system set to (10, 0, 0), and omega=(0, 10, 0). However, the actual behaviour when solved has the z component generated correctly, but both z and x components rapidly decay.
Any help, or general tips on what I’ve done wrong would be greatly appreciated. Thanks.
D\frac{\partial^2 u}{\partial x^2}- A u + \boldsymbol{\omega} \times \mathbf{u} = 0
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
mesh = BoxMesh(Point(0, 0, 0), Point(10, 3, 3), 100, 6, 6)
# Define the function space and trial/test functions
V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V) #u represents the spin polarisation of electric current in the system. It is a vector which should precess in amagnetic field.
v = TestFunction(V)
#A is a damping coefficient, set to zero for the minimum working example
A=0
#OMEGA is an externally applied magnetic field, in this case in the y direction
OMEGA=[0, 10, 0]
# Define the weak forms of the equations
a = inner(grad(u), grad(v))*dx(mesh) - A*inner(u,v)*dx(mesh) + inner(cross(as_vector(OMEGA), u), v)*dx(mesh)
# Define a SubDomain for the left side of the system
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
# Create the boundary condition
u_bc = Expression(('10', '0.0', '0.0'), degree=1)
bcs = DirichletBC(V, u_bc, LeftBoundary())
# Solve the system with the boundary condition
solve(a == 0, u, bcs)
# Extract the components of u
ux, uy, uz = u.split()
# Plot the components of u
for component, label in zip([ux, uy, uz], ['x-component', 'y-component', 'z-component']):
plt.figure()
p = plot(component)
plt.colorbar(p)
plt.title(label)