indices of boundary dofs

[erexshah]

Dear developers,

I'm dealing with a diffusion-reaction problem in which the diffusion equation is a pde that is solved in a rectangle domain for concentrations. Reaction equations are non-pde and fsolve from scipy.optimize is used to solve the system of algebraic equations.

The reactions only occur at the bottom and top boundaries; therefore, I need the indices of boundary dofs to loop only over them but not the whole domain. I tried to use par_loop but could not use the external fsolve.c function in the loopy kernel. Now, I transfer the ufl function to ndarray NumPy and loop over all dofs. It works, but the run time is not efficient; e.g., the CPU run time for the solution of diffusion equation without reaction is about 20 seconds but for the diffusion-reaction is about 9 mins. I thought at least the loop size could be decreased and limited to the boundary dofs instead of looping over the whole dofs. Any other ideas to make the algorithm more efficient would be appreciated.

A simple example of the indices of boundary dofs (bottom and top) is as follows. I need a general method to catch these indices.

import firedrake as fd
import numpy as np

mesh = fd.UnitSquareMesh(5, 5, quadrilateral=True)
V = fd.FunctionSpace(mesh, "CG", 1)
x, y = fd.SpatialCoordinate(mesh)
f1 = fd.Function(V).interpolate(x*0+y)

indices_top = np.array([24,25,29,32,34,35])
top = f1.dat.data[indices_top]
print(top)

indices_bottom = np.array([0,3,5,9,14,20])
bottom = f1.dat.data[indices_bottom]
print(bottom)

And an example of a diffusion-reaction code is as follows:

import os
os.environ["OMP_NUM_THREADS"] = "1"
from firedrake import *
from scipy.optimize import fsolve
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation


# Diffusion and physical constants
D_m = 1e-9 # diffusion coefficient, m**2/s
pH = 7
c_H = 10**(-pH)*1e3
c_bulk = c_H/100

CFL = 0.5

Lx = 25e-9
Ly = 5e-9
nx = 50
ny = 10
h_grid = Lx/nx

# reaction constants
N_s_N = 1.67*1e-6 # surface site density [mol/m]
K_N = 1e-6 # equilibrium constant of >NH+ <--> >N + H+

# mesh generation
mesh = RectangleMesh(nx, ny, Lx, Ly, quadrilateral=True)

x, y = SpatialCoordinate(mesh)
n_v = FacetNormal(mesh)
Id = Identity(mesh.geometric_dimension()) # 2x2 Identity tensor

V = FunctionSpace(mesh, "CG", 2)
v = TestFunction(V)

# UFL functions
c_n1 = Function(V, name="u^{n+1}")
c_n = Function(V, name="u^{n}")
J_H = Function(V, name="J_{H+}")

# Initialization
alpha_NH = Function(V, name="[>NH+]").interpolate(x*0 + 1/11)
alpha_N = Function(V, name="[>N]").interpolate(x*0 + 10/11)
c_n.interpolate(x*0 + c_H)
c_init = Function(V).assign(c_n)
J_init = Function(V).assign(J_H)
alpha_N_init = Function(V).assign(alpha_N)
alpha_NH_init = Function(V).assign(alpha_NH)

# time step
dt = h_grid**2*CFL/D_m
dtc = Constant(dt)

# UFL form
F = ((c_n1-c_n)/dtc)*v*dx + D_m*dot(grad(c_n1),grad(v))*dx +\
    J_H*v*ds(3) + J_H*v*ds(4)
bc_c = DirichletBC(V, c_bulk, 2)

problem = NonlinearVariationalProblem(F, c_n1, bcs=bc_c)

solver = NonlinearVariationalSolver(problem)

# reaction function
def Rate_N(rr_, cchh_, nn_, nnhh_):
    """
    This function calculates reaction on >NH site to be used 
    with fsolve
    ---------------------
    parameters:
        cchh_: molar concentration of H+ [M]
        nn_: activity of >N [mole fraction]
        nnhh_: activity of >NH+ [mole fraction]
    ---------------------    
    Target:
        rr_: mole exchange based on the reaction [mol]
    """
    eqi_ = K_N*(nnhh_+rr_/(N_s_N*h_grid))-\
        ((cchh_*1e-3-rr_/h_grid**2)*(nn_-rr_/N_s_N*h_grid))
        
    return eqi_

def SR_N(conc_, nfrac_, nhfrac_):
    """
    This function takes the UFL functions of H+ concentration [mM],
    activity of >N [mole fraction], and activity of >NH+ 
    [mole fraction] and first calculates surface mole exchange
    in all dofs. Then, it calculates the H+ flux,
    [>N], and [>NH+] for the next time step and returns the
    UFL function of these values.
    """
    ch_ = conc_.dat.data
    n_ = nfrac_.dat.data
    nh_ = nhfrac_.dat.data
    nsize = np.size(ch_)
    r_assume = 0.0
    J_ret = np.zeros_like(ch_)
    n_ret = np.zeros_like(ch_)
    nh_ret = np.zeros_like(ch_)
    r_ret = np.zeros_like(ch_)
    for i in range(nsize):
        r_ret[i] = fsolve(Rate_N, r_assume,
                          (ch_[i], n_[i], nh_[i]))
        J_ret[i] = r_ret[i]/(dt*h_grid)
        n_ret[i] = n_[i] - r_ret[i]/(N_s_N*h_grid)
        nh_ret[i] = nh_[i] + r_ret[i]/(N_s_N*h_grid)
    JJ_ = Function(V)
    NN_ = Function(V)
    NNHH_ = Function(V)
    JJ_.dat.data[:] = J_ret
    NN_.dat.data[:] = n_ret
    NNHH_.dat.data[:] = nh_ret
    
    return [JJ_, NN_, NNHH_]   
    
cs = []
nhs = []
ns = []
jhs = []

cs.append(c_init.copy(deepcopy=True))
nhs.append(alpha_NH_init.copy(deepcopy=True))
ns.append(alpha_N_init.copy(deepcopy=True))
jhs.append(J_init.copy(deepcopy=True))

t_steps = 5e3
t_end=t_steps*dt

step=0
output_freq=5e1

for t in ProgressBar("Time step").iter(np.linspace(0.0, t_end,
                                                   int(t_end/dt))):
    solver.solve()
    SR_solution = SR_N(c_n1, alpha_N, alpha_NH)
    J_H.assign(SR_solution[0])
    alpha_N.assign(SR_solution[1])
    alpha_NH.assign(SR_solution[2])
    c_n.assign(c_n1)
    
    step += 1
    
    if step % output_freq == 0:
        cs.append(c_n.copy(deepcopy=True))
        nhs.append(alpha_NH.copy(deepcopy=True))
        ns.append(alpha_N.copy(deepcopy=True))
        jhs.append(J_H.copy(deepcopy=True))

nsp = 2**4
fn_plotter = FunctionPlotter(mesh, num_sample_points=nsp)

fig, axes = plt.subplots()
axes.set_aspect('equal')
colors = tripcolor(c_init, num_sample_points=nsp,
                   vmin=c_bulk, vmax=c_H, axes=axes)
fig.colorbar(colors)

def animate(c):
    colors.set_array(fn_plotter(c))
    
animation = FuncAnimation(fig, animate,
                          frames=cs)

animation.save("reaction(start_equil).mp4", writer="ffmpeg")

Thanks in advance.

[wence-]

You can build a bc = DirichletBC(V, 0, (mesh, labels, you, want)) and then bc.nodes gives you indices into the .data array that are on the boundary.