1D Chemical Reaction System


105
views
0
7 weeks ago by
Hello,
I am trying to solve the currently available example: https://fenicsproject.org/pub/tutorial/html/._ftut1010.html#ftut1:reactionsystem in a 1D domain (UnitIntervalMesh). I'm new to FEniCs, and therefore there is much I don't understand yet. My final goal is to solve the problem stated in this article: https://www.ncbi.nlm.nih.gov/pubmed/21956887.
Several assumptions are to be made: flow velocity is constant, a homogenous system with no diffusion limitations.
So far I have the following code:

from fenics import *

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

# Read mesh from file
mesh = UnitIntervalMesh(50)

# 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 + dot(w, grad(u_1))*v_1*dx + K*u_1*u_2*v_1*dx \
+ ((u_2 - u_n2) / k)*v_2*dx + dot(w, grad(u_2))*v_2*dx + K*u_1*u_2*v_2*dx \
+ ((u_3 - u_n3) / k)*v_3*dx + dot(w, grad(u_3))*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
vtkfile_u_1 = File('reaction_system/u_1.pvd')
vtkfile_u_2 = File('reaction_system/u_2.pvd')
vtkfile_u_3 = 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()
vtkfile_u_1 << (_u_1, t)
vtkfile_u_2 << (_u_2, t)
vtkfile_u_3 << (_u_3, t)

# Update previous solution
u_n.assign(u)

However, I'm receiving the following error:

Dot product requires non-scalar arguments, got arguments with ranks 0 and 1.
Traceback (most recent call last):
File "adrs.py", line 42, in <module>
- f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx
File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 163, in dot
return Dot(a, b)
File "/usr/lib/python2.7/dist-packages/ufl/tensoralgebra.py", line 204, in __new__
"got arguments with ranks %d and %d." % (ar, br))
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Dot product requires non-scalar arguments, got arguments with ranks 0 and 1.

I think that is something related to w variable which should be a vector, however, I don't know how to state it neither its correct dimensions .
Thank you for your lost minutes,
Andre

Community: FEniCS Project

2 Answers


2
7 weeks ago by
w is a constant scalar but you try to multiply it with a gradient of something (a vector) using the dot() operator. Change w to a vector, i.e., w = Constant((1.0, 0.0)). In 1D, instead of grad(u), you can use u.dx(0) which is already a scalar.
When replacing both w = Constant((1,1,1)) and grad(u_1) by u.dx(0) for u_1,2,3, I got the following error:

Index out of bounds.
Traceback (most recent call last):
File "adrs.py", line 42, in <module>
- f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx
File "/usr/lib/python2.7/dist-packages/ufl/exproperators.py", line 512, in _dx
return d.__getitem__((Ellipsis,) + ii)
File "/usr/lib/python2.7/dist-packages/ufl/exproperators.py", line 449, in _getitem
all_indices, slice_indices, repeated_indices = create_slice_indices(component, shape, self.ufl_free_indices)
File "/usr/lib/python2.7/dist-packages/ufl/index_combination_utils.py", line 166, in create_slice_indices
error("Index out of bounds.")
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Index out of bounds.
written 7 weeks ago by André Guerra  
I forgot to replace the dot operator by *. Just dumb me. Ty this worked.
written 7 weeks ago by André Guerra  
1
7 weeks ago by
vlm  
Hi Andre,

Your problem is that you write w as a scalar by writing
w = Constant(1)
and your trying to make a dot product between it and the gradient of your variables u_1, u_2, u_3
The rank mismatches.
If you want to make w a vector, simply add more dimension, for example
w = Constant((1,1,1))
(Identity vector)

Hope it helps !
Hello, vim.
That was my first assumption but when I replaced w with Constant((1,1,1)) the following error appeared:

Dimension mismatch in dot product.
Traceback (most recent call last):
File "adrs.py", line 42, in <module>
- f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx
File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 163, in dot
return Dot(a, b)
File "/usr/lib/python2.7/dist-packages/ufl/tensoralgebra.py", line 206, in __new__
error("Dimension mismatch in dot product.")
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Dimension mismatch in dot product.

Is there a way to find the right dimensions? Regards, Andre
written 7 weeks ago by André Guerra  
1
My bad ! My answer was not very clear and I did not notice you were in 1-D
w = Constant((1,1,1)) creates the 3-D vector Identity
For your problem, I suggest you to follow the advice of Adam. For example, the foolowing form compiled for me :

# 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

(I got rid of the dot product since all variables are in 1 D)

Regards
written 7 weeks ago by vlm  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »