### Plasma-turbulence

I am attempting to get a simple 2-d plasma turbulence model running to see how fenics treats it. The pdf below gives the initial system and the weak form, for this example I am just using neumann boundary conditions along with an initial particle density (variable n) gradient. The code below that shows what I currently have, I used examples from the fenics website as a guide on how to solve coupled systems. The code currently runs it just doesn't give the zonal flow behavior that I expect to see. I am hoping that someone will be kind enough to take a look at what I've done and tell me where I've gone wrong (or if I'm completely off base in the approach that I've used). Thank you for your time and advice.

File attached: math.pdf (62.22 KB)

import random

import mshr

from dolfin import *

from fenics import *

from mshr import *

# initial condition class

class InitialCondition(UserExpression):

def eval_cell(self, value, x, ufc_cell):

value[0] = (x[0] + 100*(.5-random.random()))/10

value[1] = 10*(.5-random.random())

value[2] = 10*(.5-random.random())

def value_shape(self):

return(3,)

T = 5.0 # final time

num_steps = 50000 # number of time steps

dt = T / num_steps # time step size

eps = 0.01 # diffusion coefficient

K = 10.0 # reaction rate

## create mesh

channel = Rectangle(Point(0,0),Point(64,64))

domain = channel

mesh = generate_mesh(domain,64)

# Define function space

P1 = FiniteElement('P', triangle, 2)

element = MixedElement([P1, P1, P1])

V = FunctionSpace(mesh, element)

Q = FunctionSpace(mesh, "DG", 0) #for derivatives

# Define test functions

v_1, v_2, v_3 = TestFunctions(V)

# Define functions for velocity and concentrations

u = Function(V)

u_n = Function(V)

u.interpolate(InitialCondition())

u_n.interpolate(InitialCondition())

# Split system functions to access components

u_1, u_2, u_3 = split(u)

u_n1, u_n2, u_n3 = split(u_n)

# Define expressions used in variational forms

k = Constant(dt)

# Define variational problem

F = ((u_1 - u_n1) / k)*v_1*dx + project(u_3.dx(0), Q)*project(u_1.dx(1), Q)*v_1*dx - project(u_3.dx(1), Q)*project(u_1.dx(0), Q)*v_1*dx \

+ u_1*v_1*dx - dot(grad(u_1), grad(v_1))*dx - u_3*v_1*dx - .5*project(u_3.dx(1), Q)*v_1*dx \

+ ((u_2 - u_n2) / k)*v_2*dx + project(u_3.dx(0), Q)*project(u_2.dx(1), Q)*v_2*dx - project(u_3.dx(1), Q)*project(u_2.dx(0), Q)*v_2*dx \

+ dot(grad(u_2), grad(v_2))*dx - u_3*v_2*dx + u_1*v_2*dx \

- dot(grad(u_3), grad(v_3))*dx - u_2*v_3*dx

# Create VTK file for visualization output

vtkfile_u = File('u4.pvd')

# Time-stepping

t = 0

for n in range(num_steps):

# Update current time

t += dt

# Solve variational problem for time step

solve(F == 0, u)

# Save vector solution

vtkfile_u << (u, t)

# Update previous solution

u_n.assign(u)

# Hold plot

interactive()

As the code is currently written, the projections will only be executed once, when the residual is defined, and will generate new Function objects in the space Q. These new Functions are no longer connected to the unknown solution, and will still contain the projections of the initial condition's derivatives when the solution is updated.