initial variables values in the whole mesh in a reactional system

7 weeks ago by

Currently, I'm trying to integrate a system of PDEs which simulate the chemical species concentration at the end of the mesh. The system is composed of 3 chemical species and current is applied in a 1D domain. The data I receive (from another simulation) is a vector of the initial concentration of the chemical species every couple seconds (I have the value of the elapsed time) at the entry point of the mesh. I want to run the system in such way that I will be able to update the concentrations of the chemical species in the whole mesh prior to solving the system and returning the same structure for a posterior update. So far, I have the following code working properly:

from fenics import *

T = 5.0 # final time
num_steps = 5 # number of time steps
dt = T / num_steps # time step size
K = 10.0 # reaction rate

# Read mesh from file
mesh = UnitIntervalMesh(5)

# Define function space for system of concentrations
P1 = FiniteElement('P', interval, 1)
element = MixedElement([P1, P1, P1])
V = FunctionSpace(mesh, element)

# 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 source terms
f_1 = Constant(1)
f_2 = Constant(1)
f_3 = Constant(0)

# Define expressions used in variational forms
k = Constant(dt)
K = Constant(K)
w = Constant(1)

# Define variational problem
F = ((u_1 - u_n1) / k)*v_1*dx + w*u_1.dx(0)*v_1*dx + K*u_1*u_2*v_1*dx \
+ ((u_2 - u_n2) / k)*v_2*dx + w*u_2.dx(0)*v_2*dx + K*u_1*u_2*v_2*dx \
+ ((u_3 - u_n3) / k)*v_3*dx + w*u_3.dx(0)*v_3*dx - K*u_1*u_2*v_3*dx + K*u_3*v_3*dx \
- f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx

# Create VTK files for visualization output
vtk1 = File('reaction_system/u_1.pvd')
vtk2 = File('reaction_system/u_2.pvd')
vtk3 = File('reaction_system/u_3.pvd')

t = 0
for n in range(num_steps):

# Update current time
t += dt

# Solve variational problem for time step
solve(F == 0, u)

# Save solution to file (VTK)
_u_1, _u_2, _u_3 = u.split()

vtk1 << (_u_1, t)
vtk2 << (_u_2, t)
vtk3 << (_u_3, t)

# Update previous solution

Under the tutorial, it says that is possible to update the initial conditions at t=0 by replacing u_0 and projecting:

u_0 = Expression(('sin(x[0])', 'cos(x[0]*x[1])', 'exp(x[1])'), degree=1)
u_n = project(u_0, V)

However, my goal is to "fill" the whole mesh with the discrete values obtained in the previous simulation.
How could I manage that? How to assign u and u_n properly (with an array)?

Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »