### 1D Chemical Reaction System

105

views

0

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:

However, I'm receiving the following error:

Thank you for your lost minutes,

Andre

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

`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.
written
7 weeks ago by
André Guerra

1

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 !

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

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

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.

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.