Issue for assigning initial condition in Fenicsx

General
1

Hey. I have a really simple question. I am trying to convert my code from fenics to fenicsx but I am not very familiar with the new syntax even after reading the tutorial.

I am trying to assign an initial value to u_old, which has two elements, but I got error from this

import ufl
import random
from dolfinx.mesh import CellType, create_interval
from dolfinx.fem import Function, FunctionSpace
from ufl import dx
from mpi4py import MPI

mesh = create_interval(MPI.COMM_WORLD, 1000, [0, 1])
P1 = ufl.FiniteElement("CG", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1 * P1)
q, v = ufl.TestFunctions(ME)
u_old = Function(ME)
u_new = Function(ME)
theta_old, vy_old = ufl.split(u_old)
theta_new, vy_new = ufl.split(u_new)

# Interpolate initial condition
u_old.sub(0).interpolate(1.57+0.01*random.gauss(0, 1))
u_old.sub(1).interpolate(0)

And the error is

IndexError: invalid axis: 0 (ndim = 0)

How could I fix this? Any help and suggestions are appreciated.

2

You need to provide a function within interpolate. Try this:

u_old.sub(0).interpolate(lambda x: numpy.full((x.shape[1],), 1.57+0.01*random.gauss(0,1)))
u_old.sub(1).interpolate(lambda x: numpy.full((x.shape[1],), 0))
1 Like
3

Thank you very much! Can I ask you one more about the boundary condition? Is there a quick way to set the end point a constant instead of using the topology part from tutorial? Here I only need u_old.sub(0)= 1.57 and u_old.sub(1) = 0 at two end points.

4

You can use dolfinx.fem.locate_dofs_topological or dolfinx.fem.locate_dofs_geometrical on a collapsed subspace to get the dofs in question, and use the map from the collapsed sub-space to assign the value directly.

from dolfinx.mesh import create_interval
from dolfinx.fem import Function, FunctionSpace, locate_dofs_geometrical
from dolfinx.io import XDMFFile
import ufl
import numpy as np
from mpi4py import MPI

mesh = create_interval(MPI.COMM_WORLD, 1000, [0, 1])
P1 = ufl.FiniteElement("CG", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1 * P1)
u_old = Function(ME)

M0, map0 = ME.sub(0).collapse()
M1, map1 = ME.sub(1).collapse()

dof_left = locate_dofs_geometrical(M0, lambda x: np.isclose(x[0], 1))
val_left = 1
assert len(dof_left == 1)
u_old.x.array[map0[dof_left[0]]] = 3

with XDMFFile(MPI.COMM_WORLD, "u0.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(u_old.sub(0).collapse())
1 Like
5

Thank you very much for your reply. but I still have some stupid questions. Sorry about this but I really want to figure out what is going on here.

My idea is since u-old is p1*p1 function and it has two parts, one is theta_old and the other is vy-old. For theta_old I try to assign 1.57 at x=0 and x=1. For vy_old I try to assign 0 at x=0 and x=1.

In your code, I am not sure how to distinguish which value is assigned to which function. Or are you just treating u_old as a whole?

Thanks again for your help!

6

There was a typo in my initial code that Ive corrected now. It should have been

You find dofs of theta-old with M0 and map0, while dofs from vy-old with M1 and map1

7

Aha. That makes sense. I have one more question. In your code you define

val_left = 1

but later you don’t use it. Would you explain a little more on this one? Thank you very much.

8

Just ignore it. The value set is 3

9

I see. Thank you very much! Have a good day!