Hi,
I have a similar question to this particular post. I’m wondering how to interpolate special functions into dolfinx.
Unlike other posts I’ve seen however I’m wondering how to make them dependent on the values of other fem.Function objects.
Let me show a motivating example. For ease of understanding I’ve copied most of the diffusion example and implemented a trivial function as a test case.
import scipy.constants as con
import ufl
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import (
apply_lifting,
assemble_matrix,
assemble_vector,
create_vector,
set_bc,
)
from mpi4py import MPI
from petsc4py import PETSc
# Define temporal parameters
t = 0 # Start time
T = 1.0 # Final time
num_steps = 50
dt = T / num_steps # time step size
# Define mesh
nx, ny = 50, 50
domain = mesh.create_rectangle(
MPI.COMM_WORLD,
[np.array([-2, -2]), np.array([2, 2])],
[nx, ny],
mesh.CellType.triangle,
)
V = fem.functionspace(domain, ('Lagrange', 1))
def ElectronHeatCapacity(T, Ne):
"""Imagine a more complex funciton here. In reality there are fermi-dirac integrals ect.
For a specific example look at npplus fermi code.
"""
return 3 / 2 * con.k
def initial_condition(x, a=5):
return np.exp(-a * (x[0] ** 2 + x[1] ** 2))
u_n = fem.Function(V)
u_n.name = 'u_n'
u_n.interpolate(initial_condition)
cp = fem.Function(V)
cp.name = 'cp'
ne = fem.Function(V)
ne.name = 'ne'
ne.interpolate(lambda x: np.ones_like(x[1]) * 4e21)
# what I'm currently doing.
cp.x.array[:] = ElectronHeatCapacity(uh.x.array, ne.x.array)
# Create boundary condition
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(
domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool)
)
bc = fem.dirichletbc(
PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V
)
xdmf = io.XDMFFile(domain.comm, 'diffusion.xdmf', 'w')
xdmf.write_mesh(domain)
# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = fem.Function(V)
uh.name = 'uh'
uh.interpolate(initial_condition)
xdmf.write_function(uh, t)
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))
a = cp * u * v * ufl.dx + dt * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = (u_n + dt * f) * v * ufl.dx
bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = create_vector(linear_form)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
for i in range(num_steps):
t += dt
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, linear_form)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])
# Solve linear problem
solver.solve(b, uh.vector)
uh.x.scatter_forward()
## update the heat capacity
cp.x.array[:] = ElectronHeatCapacity(uh.x.array, ne.x.array)
# Update solution at previous time step (u_n)
u_n.x.array[:] = uh.x.array
# Write solution to file
xdmf.write_function(uh, t)
# Update plot
xdmf.close()
The only real change I’ve implemented here is the Cp specific heat fem.Function. The cp depends on the Ne the electron density and the current temperature. I have a function I call that gets that takes in those variables and return the specific heat. (Right now that function is trivial but imagine I’m calling some fermi-dirac integral solver function)
My current approach is to just directly operate of the DOF values. I’m curious if this I the best approach? Would it scale properly when using multiple MPI processes for example? Also are there any advantages to using quadrature elements for the internal values instead?
Thanks for any and all suggestions.