Question about how to define a function

[guangyi_cai]

Im a rookie in Fenicsx, Im confused about how to define the function phi_0 in picture 3 so that I can put phi_0 into a ( u , v ) and do some advanced work.

my code is as followed:

# create a mesh
msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], CellType.triangle
)

#function to mark x=0
def in_boundary(x):
    return np.isclose(x[0],0.0)

#function to mark y=0,y=1
def wall_boundary(x):
    return np.logical_or(
        np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))
    )

#function to mark x=1
def in_boundary(x):
    return np.isclose(x[0],1.0)

# in_velocity
def in_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))

#define V and Q
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)
V, Q = functionspace(msh, P2), functionspace(msh, P1)

# No-slip condition on boundaries where y=0,y=1
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)
facets = locate_entities_boundary(msh, 1, wall_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)

# Driving velocity condition on boundary where x=0
in_velocity = Function(V)
in_velocity.interpolate(in_velocity_expression)
facets = locate_entities_boundary(msh, 1, in_boundary)
bc1 = dirichletbc(in_velocity, locate_dofs_topological(V, 1, facets))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1]

# alpha_0
alpha_0=Constant(mesh, PETSc.ScalarType(1))

# phi_0


# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))

a = form([[inner(grad(u), grad(v)) * dx + inner(alpha_0*(1-phi_0)*(1-phi_0)*u,v)*dx, -inner(p, div(v)) * dx], [inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx]) 

I would appreciate your help!

[dokken]

It seems like you want phi_0 to be a cellwise constant, i.e. Create a function in DG-0 that is 1 within the respective area and 0 elsewhere. See:
https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html

[guangyi_cai]

Appreciate your help！
So the code would be:

#define phi_0
W = functionspace(msh, ("DG", 0))

def condition_0(x):
    return np.logical_or(x[0] <= 0.25, np.logical_and(x[1] >= 1/3, x[1] <= 2/3))


def condition_1(x):
    return np.logical_or(np.logical_and(x[0] > 0.25,x[1] < 1/3),np.logical_and(x[0] > 0.25,x[1] > 2/3))

phi_0 = Function(W)
cells_0 = locate_entities(msh, msh.topology.dim, condition_0)
cells_1 = locate_entities(msh, msh.topology.dim, condition_1)

phi_0.x.array[cells_0] = np.full_like(cells_0, 1, dtype=default_scalar_type)
phi_0.x.array[cells_1] = np.full_like(cells_1, 0, dtype=default_scalar_type)

# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))

a = form([[inner(grad(u), grad(v)) * dx + inner(alpha_0*(Constant(msh, PETSc.ScalarType(1))-phi_0)*(Constant(msh, PETSc.ScalarType(1))-phi_0)*u,v)*dx, -inner(p, div(v)) * dx], [inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx]) 

I’m not sure if my code is correct.

[dokken]

You can check for correctness by visualizing phi_0 with either paraview or pyvista.

[guangyi_cai]

I have Defined variational problem:

# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))

a = form([[inner(grad(u), grad(v)) * dx + inner(alpha_0*(Constant(msh, PETSc.ScalarType(1.0))-phi_0)*(Constant(msh, PETSc.ScalarType(1.0))-phi_0)*u, v)*dx, -inner(p, div(v)) * dx], [inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx]) 

How can I solve the PDE in picture 1 with Taylor-Hood element,
is there some references for me?

[dokken]

It is unclear what your current issue is? What happens if you follow the demo
https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo_stokes.html

Do you get an error message.?

[guangyi_cai]

I have a few questions about the demo：

1.we can see “-inner(p, div(v)) * dx” in the weak formulation given in the document,
and the code from the demo:

a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx],
          [inner(div(u), q) * dx, None]])

why there is not a negative sign in front of " inner(p, div(v)) * dx" ?

2. Under this boundary conditions, I get the right output:

# Function to mark y = 0, y = 1, x = 1
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0)),np.isclose(x[0], 1.0))


# Function to mark the lid (x = 0)
def lid(x):
    return np.isclose(x[0], 0.0)


# Lid velocity
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))

# No-slip condition on boundaries where y = 0, y = 1, x = 1
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)

# Driving (lid) velocity condition on top boundary (x = 0)
lid_velocity = Function(V)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
bc1 = dirichletbc(lid_velocity, locate_dofs_topological(V, 1, facets))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1]

But once I want to remove the dirichlet boundary condition of x = 1 and I change the code into:

# Function to mark y = 0, y = 1
def noslip_boundary(x):
    return np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))


# Function to mark the lid (x = 0)
def lid(x):
    return np.isclose(x[0], 0.0)

(the rest parts of the code is as before)

I get an issue:

Traceback (most recent call last):
  File "/home/fenics/work/demo3.py", line 226, in <module>
    norm_u_0, norm_p_0 = nested_iterative_solver()
  File "/home/fenics/work/demo3.py", line 157, in nested_iterative_solver
    assert nsp.test(A)
AssertionError

I think I set the boundary conditions wrongly but I don’t know how to modify it.

[dokken]

@guangyi_cai

why there is not a negative sign in front of " inner(p, div(v)) * dx" ?

One usually solves for the negative pressure to ensure a symmetric system.

@guangyi_cai

But once I want to remove the dirichlet boundary condition of x = 1 and I change the code into:

# Function to mark y = 0, y = 1
def noslip_boundary(x):
    return np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))


# Function to mark the lid (x = 0)
def lid(x):
    return np.isclose(x[0], 0.0)

Just to make sure, here you are trying to enforce
u(0,y)=(1,0),
u(x,0)=u(x,1)=(0,0)
-du/dn+pn=0 at (1,y), i.e. A square channel with (1,0) at the left inlet, outlet at right boundary and top and bottom being zero?

(Im currently not at a computer, so I cannot test your code).

Is the error messages stemming from you modifying the referenced demo or your code above ?

[guangyi_cai]

The boundary conditions you mentioned is exactly what I want to set up!
And I just changed the code snippet from the referenced demo demo_stokes.py given in the document

# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)),
                         np.isclose(x[1], 0.0))


# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)

into:

# Function to mark y = 0, y = 1
def noslip_boundary(x):
    return np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))


# Function to mark the lid (x = 0)
def lid(x):
    return np.isclose(x[0], 0.0)

I tried to run the demo I modified and then I got the error messages I mentioned above.

[dokken]

The error message is for testing the null space.
Simply removing the null space (as you no longer have a nullspace, as you have a natural boundary condition. I.e. remove

    nsp = PETSc.NullSpace().create(vectors=[null_vec])
    assert nsp.test(A)
    A.setNullSpace(nsp)

this from the two places it is used.

[guangyi_cai]

Appreciate your help again！
There are still some problem:
I want to put the u that was solved out in the demo into my advanced work i.e.:

And there is also a constrain condition :
The solutions are mismatched between the velocity field and the phase field due to employing different polynomials as the basis functions on the same mesh. Projection and interpolationare simply used to transport datas:

I am confused about how to import u that was solved out in the demo into the weak formulation in my first picture.

[dokken]

@guangyi_cai

I am confused about how to import u that was solved out in the demo into the weak formulation in my first picture.

You can use adios4dolfinx, see
http://jsdokken.com/adios4dolfinx/README.html

[guangyi_cai]

I just installed adios4dolfinx 0.7.2 and tried to run the domo writing_functions_checkpoint.py downloaded from http://jsdokken.com/adios4dolfinx/docs/writing_functions_checkpoint.html
And then i got an error:

Traceback (most recent call last):
  File "/home/fenics/work/writing_functions_checkpoint.py", line 41, in <module>
    adios4dolfinx.write_mesh(filename, mesh)
  File "/usr/local/lib/python3.10/dist-packages/adios4dolfinx/checkpointing.py", line 86, in write_mesh
    num_xdofs_local = mesh.geometry.index_map().size_local
AttributeError: 'PosixPath' object has no attribute 'geometry'

[dokken]

Switch the order of the arguments. The webpage is based on v0.8.x of adios4dolfinx, and we changed the API to be consistent across all functions

Make API consistent (filename first) in read_ write_ operations. Vectorize dofmap reading from legacy

See api changes in Release v0.8.0 · jorgensd/adios4dolfinx · GitHub

[guangyi_cai]

I just referred to https://jsdokken.com/checkpointing-presentation/#/10
and just simply add

#save u with adios4dolfinx
    checkpoint_file = pathlib.Path("function_checkpoint.bp")
    adx.write_mesh(msh, checkpoint_file, engine="BP4")
    adx.write_function(u, checkpoint_file, engine="BP4")

behind

# Create finite element {py:class}`Function <dolfinx.fem.Function>`s
    # for the velocity (on the space `V`) and for the pressure (on the
    # space `Q`). The vectors for `u` and `p` are combined to form a
    # nested vector and the system is solved.
    u, p = Function(V), Function(Q)
    x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(u.x),
                                la.create_petsc_vector_wrap(p.x)])
    ksp.solve(b, x)

in demo_stokes.py
downloaded fromhttps://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo_stokes.html#:

# Create finite element {py:class}`Function <dolfinx.fem.Function>`s
    # for the velocity (on the space `V`) and for the pressure (on the
    # space `Q`). The vectors for `u` and `p` are combined to form a
    # nested vector and the system is solved.
    u, p = Function(V), Function(Q)
    x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(u.x),
                                la.create_petsc_vector_wrap(p.x)])
    ksp.solve(b, x)

    #save u with adios4dolfinx
    checkpoint_file = pathlib.Path("function_checkpoint.bp")
    adx.write_mesh(msh, checkpoint_file, engine="BP4")
    adx.write_function(u, checkpoint_file, engine="BP4")

I’m not sure if my operation is correct (since i have seen many very complex codes in other demos)

Furthermore,if i want to import u into my advanced work ,can i just simply use the code here？

from mpi4py import MPI
import adios4dolfinx
from pathlib import Path
import dolfinx
checkpoint_file = Path("function_checkpoint.bp")
msh = adios4dolfinx.read_mesh(MPI.COMM_WORLD, checkpoint_file, engine='BP4', ghost_mode=dolfinx.mesh.GhostMode.shared_facet)

from basix.ufl import element
from dolfinx.fem import  functionspace, Function
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
V = functionspace(msh, P2)
u = Function(V)
adios4dolfinx.read_function(u, checkpoint_file, engine="BP4")

[dokken]

Yes, that would be the general idea.

[guangyi_cai]

Appreciate for your patience！
Now I want to use projection to transport data (Maybe from P2 to P1) ,so that I can complete the follow-up calculation of the weak formulation below：

Im confused about how to do the projection work.Here is my code:

#read mesh
checkpoint_file = Path("function_checkpoint.bp")
msh = adios4dolfinx.read_mesh(MPI.COMM_WORLD, checkpoint_file, engine='BP4', ghost_mode=dolfinx.mesh.GhostMode.shared_facet)

#read u
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
U = functionspace(msh, P2)
u = Function(U)
adios4dolfinx.read_function(u, checkpoint_file, engine="BP4")

P1 = element("Lagrange", msh.basix_cell(), 1)
V = functionspace(msh, P1)

#projection
u_projected = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a_proj = inner(u_projected, v) * dx
L_proj = inner(u, v) * dx

A_proj = assemble_matrix(a_proj)
A_proj.assemble()
b_proj = assemble_vector(L_proj)
solver_proj = PETSc.KSP().create(msh.comm)
solver_proj.setOperators(A_proj)
solver_proj.solve(b_proj, u_projected.vector)
u_projected.x.scatter_forward()

And it failed with:

Traceback (most recent call last):
  File "/home/fenics/work/phase.py", line 33, in <module>
    L_proj = inner(u, v) * dx
  File "/usr/local/lib/python3.10/dist-packages/ufl/operators.py", line 167, in inner
    return Inner(a, b)
  File "/usr/local/lib/python3.10/dist-packages/ufl/tensoralgebra.py", line 162, in __new__
    raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
ValueError: Shapes do not match: <Coefficient id=140108284108736> and <Argument id=140108318196832>

It seems that there is a shape mismatch between the operands in the inner function …? How should I modify the code ?
And I alse referred to L2 projection of fine mesh solution to a coarse mesh ,maybe I should interpolate the function u first? If so ,how should I interpolate the function?

[dokken]

@guangyi_cai

P1 = element("Lagrange", msh.basix_cell(), 1)
V = functionspace(msh, P1)

You need to supply a shape to: i.e. make it a vector space like in

@guangyi_cai

P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
U = functionspace(msh, P2)

[guangyi_cai]

Let’s return to the discussion of solving stokes equation,now I want to Update the velocity field u^(n+1) by the following state equation using the known phase-field function in the previous step.

where α(Φ) = α_0 (1 - Φ) (1 - Φ) and the boundary conditions are:

@dokken

u(0,y)=(1,0)u(0,y)=(1,0),
u(x,0)=u(x,1)=(0,0)u(x,0)=u(x,1)=(0,0)u(x,0)=u(x,1)=(0,0)
-du/dn+pn=0−du/dn+pn=0-du/dn+pn=0 at (1,y)

I just modified demo_stokes.py downloaded from https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo_stokes.html

The modifications I made were:
1.changed boundary condituons in the demo
2.defined the Φ function
3.set

(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))  # type: ignore
a = form([[inner(grad(u), grad(v)) * dx + inner(alpha_0*(Constant(msh, PETSc.ScalarType(1.0))-phi_0)*(Constant(msh, PETSc.ScalarType(1.0))-phi_0)*u, v)*dx, inner(p, div(v)) * dx], [inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx , inner(Constant(msh, PETSc.ScalarType(0)), q) * dx]) 

4.simply removed

@dokken

    nsp = PETSc.NullSpace().create(vectors=[null_vec])
    assert nsp.test(A)
    A.setNullSpace(nsp)

in the demo

Actually I got a solution by running the code I modified but there are some differences between the result of the output and the expected result…
Maybe I set a wrongly or maybe there are more modifications needed to be done in the defined function nested_iterative_solver(): ?

[dokken]

Please provide a minimal reproducible example that uses only one of the four implementation given in

@guangyi_cai

I just modified demo_stokes.py downloaded from Stokes equations using Taylor-Hood elements — DOLFINx 0.7.3 documentation