10 days ago by
Hi there,

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):

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)


# 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

# Hold plot

Community: FEniCS Project
I'm not really sure what the desired flow pattern is, but one issue that jumps out at me is the use of project() in the definition of the residual, which is probably not doing what you want it to.  You should be able to just directly use the desired partial derivatives in the residual without any projection.  

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.  
written 10 days ago by David Kamensky  
Thank you, you are right, that was not my intention. The syntax for a partial derivative that will be updated each iteration is just, for example, u_2.dx(1) then?
written 10 days ago by kirk tolfa  
Yes, you should be able to use that directly in the weak form.  
written 10 days ago by David Kamensky  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »