Assigning initial conditions based on radial coordinates

martamestroni · October 4, 2022, 12:25pm

Hi all,

I am trying to assign initial conditions in a circular mesh which I am then using for solving the diffusion equation. What I want to obtain is a constant initial value in the whole mesh except for the most external nodes lying on the external-most circumference, to which I would like to prescribe a specific different value (B_value).

I am using a circular mesh created on Gmsh saved in .msh ASCII version 2 file format ($MeshFormat: 2.2 0 8).

This is a working MWE:

# Initializing parameters
t = 0 # Start time
T = 300 # End time
num_steps = 20 # Number of time steps
dt = (T-t)/num_steps # Time step size
alpha = 3
beta = 1.2
R=1
B_value =10000000
# Diffusion coefficient
D_=0.001

from re import X
import numpy as np
import meshio
import gmsh
import sys
import pygmsh
from dolfinx import mesh, fem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.fem import FunctionSpace, Function, Constant, locate_dofs_geometrical,locate_dofs_topological, dirichletbc, form, Expression
import dolfinx.mesh
import ufl
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
domain, cell_markers, facet_markers = gmshio.read_from_msh("circle.msh", MPI.COMM_WORLD, gdim=2)

V = fem.FunctionSpace(domain, ("CG", 1))
u_n = fem.Function(V)
u_n.interpolate(lambda x: (B_value*np.trunc(np.sqrt(x[0]**2+x[1]**2))))


# BOUNDARY condtion settng Neumann at the center and boarders
boundaries =[(1, lambda x: np.isclose(np.sqrt(x[0]**2+x[1]**2),R)),
            (2, lambda x: np.isclose(np.sqrt(x[0]**2+x[1]**2),0))]

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = dolfinx.mesh.locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with XDMFFile(domain.comm, "facet_tags.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tag)

ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

f_=0
f= Constant(domain, ScalarType(f_))

D= Constant(domain, ScalarType(D_))

# creating file to store solution
xdmf = XDMFFile(domain.comm, "2D_Diffusion_Circle.xdmf", "w")

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
F = u*v*ufl.dx + D*dt*ufl.dot(ufl.grad(u), ufl.grad(v))*ufl.dx - (u_n + dt*f)*v*ufl.dx

# Finishing defining boundaries and apppending to F


class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == "Dirichlet":
            u_no = Function(V)
            u_no.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            self._bc = dirichletbc(u_no, dofs)
        elif type == "Neumann":
                self._bc = ufl.inner(values, v) * ds(marker)
        elif type == "Robin":
            self._bc = values[0] * ufl.inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# Define the Dirichlet condition
boundary_conditions = [BoundaryCondition("Neumann", 1, 0),
                        BoundaryCondition("Neumann", 2, 0)
]

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc


a = fem.form(ufl.lhs(F))
L = fem.form(ufl.rhs(F))

A = fem.petsc.assemble_matrix(a)
A.assemble()
b = fem.petsc.create_vector(L)
uh = fem.Function(V)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

# write value of initial condition
xdmf.write_mesh(domain)
xdmf.write_function(u_n, t)

for n in range(num_steps):
    # Update Diriclet boundary condition 
    # u_exact.t+=dt
    # u_D.interpolate(u_exact)
    
    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    fem.petsc.assemble_vector(b, L)
    
    # Apply Dirichlet boundary condition to the vector
    # fem.petsc.apply_lifting(b, [a], [[bcs]])
    # b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    # fem.petsc.set_bc(b, [bcs])

    # Solve linear problem
    solver.solve(b, uh.vector)
    uh.x.scatter_forward()

    # Update solution at previous time step (u_n)
    u_n.x.array[:] = uh.x.array
    
    # write value of computed solution at each step
    t += dt
    xdmf.write_function(u_n, t)

xdmf.close()

The way I am defining the initial condition interpolating the lambda function is the only way I was able to run the code, but it is not optimal (it assigns the B_value only to half of the boundary) and I would like a more elegant way.

I have tried some other ways to assign initial conditions, all based on the tutorials and some other examples posted here on the forum, with no success. I will briefly list them:

(1)

class initial_condition():
    def __init__(self, B_value, R):
        self.B_value = B_value
        self.R = R
    def __call__(self, x):
                if np.isclose(x[0]**2 + x[1]**2, self.R*self.R):
                    return self.B_value 
                else:
                    return 0

u_init = initial_condition(B_value, R)
u_n = fem.Function(V)
u_n.interpolate(u_init)

error: File “/root/shared/2D_Diffusion_Circle.py”, line 162, in call
if np.isclose(x[0]**2 + x[1]**2, self.R*self.R):
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

(2)

x= ufl.SpatialCoordinate(domain)
def initial_condition(x):
    if ufl.exp(np.isclose(np.sqrt(x[0]**2 + x[1]**2), 1)):
        return B_value
    else:
        return 0
expr=fem.Expression(initial_condition(x),V.element.interpolation_points())
u_n = fem.Function(V)
u_n.interpolate(expr)

Error: AttributeError: ‘Sum’ object has no attribute ‘sqrt’
(3)

def inside(x):
     return (x[0]*x[0] + x[1]*x[1] < R*R)

def boundary(x):
     return np.isclose(x[0]*x[0] + x[1]*x[1],R*R)

cells_inside = dolfinx.mesh.locate_entities(domain, domain.topology.dim, inside)
cells_boundary = dolfinx.mesh.locate_entities(domain, domain.topology.dim, boundary)

initial=fem.Function(V)
initial.x.array[cells_inside] = np.full_like(cells_inside, 0, dtype=ScalarType)
initial.x.array[cells_boundary] = np.full_like(cells_boundary, B_value, dtype=ScalarType)
u_n = fem.Function(V)
u_n.interpolate(initial)

Error: File “/root/shared/2D_Diffusion_Circle.py”, line 166, in
initial.x.array[cells_inside] = np.full_like(cells_inside, 0, dtype=ScalarType)
IndexError: index 455237 is out of bounds for axis 0 with size 455237

I guess all the error messages somehow depend on the definition of points lying on the external circuference (x^2+y^2==R) and on the implicit definition of x.
Could anyone spot what I am doing wrong or please suggest in which direction to go?

Thank you!

Marta

hawkspar · October 12, 2022, 8:39am

For what it’s worth, I don’t think your error messages relate to circular geometry but more to a wrong understanding of the types you are handling.

Your post was not taggued as dolfinx, but it seems that’s what we’re using, wo we should speak somewhat the same language. You haven’t posted your version of our beloved library though.

For (1), x is an array that spans the entire domain, and interpolate expects an array too.

def __call__(self, x:np.array) -> np.array:
    res=np.zeros_like(x[0])
    res[np.isclose(x[0]**2 + x[1]**2, self.R*self.R)]=self.B_value
    return res

As for (2), you’re misxing ufl.Expression and numpy.array objects - these two are like water and oil, they don’t mix well. In dolfinx you can directly interpolate a function that takes an array and returns an array :

def initial_condition(x:np.array) -> np.array:
    res=np.zeros_like(x[0])
    res[np.isclose(x[0]**2 + x[1]**2, 1)]=self.B_value
    return res
u_n = fem.Function(V)
u_n.interpolate(initial_condition)

For (3), I think you’re confusing locate_entities with locate_entities_geometrical. One returns facets, the other quadrature points, so I wouldn’t expect their numbering to match. I think this should work better :

def boundary(x:np.array) -> np.array:
     return np.isclose(x[0]*x[0] + x[1]*x[1],R*R)
dofs = fem.locate_dofs_geometrical(V, boundary)
u_n = fem.Function(V) # By default everything is already 0
u_n.x.array[dofs]=B_value

Hope I’m not too late.

Topic		Replies	Views
Improve efficiency when using parallel computing and limited by formulation General	0	265	March 3, 2023
Issues with Defining Boundary Conditions for Stokes Problem Using GMSH Mesh Tags dolfinx boundary-conditions , gmsh	10	162	March 19, 2024
Initial condition on part of the domain variational formulation	1	635	December 9, 2021
Setting initial condition on irregular boundary from marked input mesh General	3	333	November 12, 2022
How to supply boundary condition on a surface of a sphere dolfinx	1	289	September 13, 2023

Assigning initial conditions based on radial coordinates

Related Topics