How to resolve indices assertion error received in Adios4dolfinx integration to advection diffusion reaction dolfin example?

General
,
1

I have added adios4dolfinx in ns_code2.ipynb by forking dolfinx-tutorial repo. In particular, I am storing mesh and function at each time step.

adios4dolfinx.write_mesh(cp_filename, mesh)
adios4dolfinx.write_function(cp_filename, u_, time=t, name="u_function")
adios4dolfinx.write_function(cp_filename, p_, time=t, name="p_function")

Commit shown on this link: (Adios coupled · hassaniqbal209/dolfinx-tutorial@d8ecf22 · GitHub)).

The goal is to run with dolfinx the following code written in dolfin: https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft09_reaction_system.py.

Made some changes corresponding to migration to dolfinx (as shown in commit above) but currently receiving error in adios4dolfinx.read_function() in advdiffreac.ipynb

Can you please guide how to go about resolving this?

AssertionError Traceback (most recent call last)
Cell In[38], line 6
4 w = fem.Function(W)
5 # adios4dolfinx.read_function(“./Flows/flow_function_0”, w)
----> 6 adios4dolfinx.read_function(cp_filename, w, time=0, name=“u_function”)
7 # mkdir -p Flows/flow_function_0
8 # touch Flows/flow_function_0/md.0 Flows/flow_function_0/md.idx

File ~/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/adios4dolfinx/checkpointing.py:358, in read_function(filename, u, engine, time, legacy, name)
354 # Compute owner of dofs in dofmap
355 num_dofs_global = (
356 u.function_space.dofmap.index_map.size_global * u.function_space.dofmap.index_map_bs
357 )
→ 358 dof_owner = index_owner(comm, input_dofmap.array, num_dofs_global)
360 # --------------------Step 4-----------------------------------
361 # Read array from file and communicate them to input dofmap process
362 if legacy:

File ~/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/adios4dolfinx/utils.py:127, in index_owner(comm, indices, N)
119 “”"
120 Find which rank (local to comm) which owns an index, given that
121 data of size N has been split equally among the ranks.
(…)
124 processes gets an extra value.
125 “”"
126 size = comm.size
→ 127 assert (indices < N).all()
128 n = N // size
129 r = N % size

AssertionError:

if the degree is changed from 2 to 3 in the following

W = fem.functionspace(
    mesh=msh,
    element=element(  # For Vector Function Space using Vector Element
        family="Lagrange", cell=msh.topology.cell_name(), degree=2
    ),
)

The assertion error is resolved. However, a new error is received on

F = (
    (u1 - u_n1) / Delta_t * v1
    + ufl.dot(w, ufl.grad(u1)) * v1
    + eps * ufl.dot(ufl.grad(u1), ufl.grad(v1))
) * ufl.dx

as follows

ValueError Traceback (most recent call last)
Cell In[42], line 20
15 f2.interpolate(lambda x: source(x, np.array([0.1, 0.3]), 0.05))
16 f3.x.array[:] = 0.0
18 F = (
19 (u1 - u_n1) / Delta_t * v1
—> 20 + ufl.dot(w, ufl.grad(u1)) * v1
21 + eps * ufl.dot(ufl.grad(u1), ufl.grad(v1))
22 ) * ufl.dx
23 F += (
24 (u2 - u_n2) / Delta_t * v2
25 + ufl.dot(w, ufl.grad(u2)) * v2
26 + eps * ufl.dot(ufl.grad(u2), ufl.grad(v2))
27 ) * ufl.dx
28 F += (
29 (u3 - u_n3) / Delta_t * v3
30 + ufl.dot(w, ufl.grad(u3)) * v3
31 + eps * ufl.dot(ufl.grad(u3), ufl.grad(v3))
32 ) * ufl.dx

File ~/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/operators.py:225, in dot(a, b)
223 if a.ufl_shape == () and b.ufl_shape == ():
224 return a * b
→ 225 return Dot(a, b)

File ~/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/tensoralgebra.py:212, in Dot.new(cls, a, b)
210 # Checks
211 if not ((ar >= 1 and br >= 1) or scalar):
→ 212 raise ValueError(
213 “Dot product requires non-scalar arguments, "
214 f"got arguments with ranks {ar} and {br}.”
215 )
216 if not (scalar or ash[-1] == bsh[0]):
217 raise ValueError(“Dimension mismatch in dot product.”)

ValueError: Dot product requires non-scalar arguments, got arguments with ranks 0 and 1.

2

This is not the correct space for a velocity field/vector element.
This could be written as

W = fem.functionspace(msh, ("Lagrange", 2, (msh.geometry.dim,))

Ie a second order space with number of components equal to the dimension of the mesh geometry (1, 2 or 3). In your case this number will be 2.