I want to see if Dolfinx can solve the following Poisson Equation for a 2D fluid flow
\nabla^2 \chi = \nabla\cdot\textbf{V}
\lim_{\textbf{|x|}\rightarrow \infty}\chi(x,y) = 0
and in order to simulate the infinite space, I opted to try what was used in the electromagnetics example by increasing the size of the domain and setting the values on the boundary on the larger domain to 0. To test this result, I compared it to the result from convolution with the free space Green’s function
\chi(x,y) = \int_\Omega G(\textbf{x},\textbf{x}_0) \nabla \cdot \textbf{V}(\textbf{x}_0) d\textbf{x}_0
G(\textbf{x},\textbf{x}_0) = - \frac{1}{2\pi}\ln(|\textbf{x}- \textbf{x}_0|)
Below is a comparison of the results using real divergence data on a [0\le x\le 1; 0\le y \le 1] Cartesian square (ignore the first 3 solution plots as they used the discrete sine transform-I, fast Fourier transform, and fenicsx’s LU solver methods to get results for homogeneous Dirichlet boundary conditions on the original domain; the last two plots are our focus as they compare the results over the assumption that the solution goes to 0 over infinite space)
The lower-right plot which used fenicsx’s LU solver with homogeneous Dirichlet boundary conditions set on the boundary of the domain 5 times larger than the original to simulate an infinite space should look similar to the result from convoluting with Green’s function. Now, let’s observe a manufactured case where the solution is a Gaussian function (a function which naturally decays to 0 as space goes to infinity) centered at (0.5,0.5)
\chi(x,y) = e^{-a\big((x-0.5)^2 + (y-0.5)^2\big)}
the Poisson Equation is then
\nabla^2 \chi = -4a\Big(a\big((x-0.5)^2 + (y-0.5)^2\big) - 1 \Big) e^{-a\big((x-0.5)^2 + (y-0.5)^2\big)}
For the solution to be accurate, we need the Gaussian source to be relatively compared to the size of the domain (observe that each local maxima and minima of divergence in the data driven example were also relatively small compared to the overall area of the grid), so I have set a = 50 to produce a sharp concentric peak. Below are plots comparing the results
where it is difficult to visually see the difference, so I have checked the L2 errors. For the L2 errors of the Green’s function solution, I have:
Error_L2: 1.61e-02
Error_max: 2.55e-02
and for the LU solver which attempted to simulate an infinite space, I have
Error_L2: 8.62e-02
Error_max: 1.61e-01
so the overall the errors produced by Fenicsx’s LU solver are relatively much higher with extended space, but not ruinously so for this case. From what we can see for the case involving data, these results probably would be an issue.
Unless I messed something up in my code somewhere, it looks like the method of extending the space might not be sufficient for all cases, and overall produces significantly more error. So, I want to know if there is an alternative method to simulate infinite space which fenicsx can handle to produce more accurate results? Or maybe a different solver than LU would handle this better? Let me know anything that might work!
Here is my code:
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import fftconvolve
from scipy.fft import dst, idst
from scipy.interpolate import RegularGridInterpolator
import xarray as xr
import dolfinx as df
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import gmsh
# Helper function
def error_L2(y_ex, y_num):
difference = y_ex - y_num
l2_error = np.linalg.norm(difference, ord=2)
abs_max_error = np.linalg.norm(difference, ord=np.inf)
print(f"\nError_L2: {l2_error:.2e}")
print(f"Error_max: {abs_max_error:.2e}\n")
# Matplotlib plotter functions ############################################################
def cells_at_xy(X_grid,Y_grid, mesh, return_points):
points = np.vstack((X_grid.ravel(), Y_grid.ravel(), np.zeros(X_grid.size))).T
bb_tree = df.geometry.bb_tree(mesh, mesh.topology.dim)
cell_indices = []
for point in points:
cells_candidates = df.geometry.compute_collisions_points(bb_tree, point)
colliding_cells = df.geometry.compute_colliding_cells(mesh, cells_candidates, point)
cell_index = colliding_cells.links(0)[0]
cell_indices.append(cell_index)
if return_points:
return points, cell_indices
else:
return cell_indices
def xy_grid(mesh):
#Nx = len(mesh.geometry.x[:,0])
#Ny = len(mesh.geometry.x[:,1])
Nx = np.arange(mesh.geometry.x[:,0].min(),mesh.geometry.x[:,0].max(),h).shape[0]
Ny = np.arange(mesh.geometry.x[:,1].min(),mesh.geometry.x[:,1].max(),h).shape[0]
x_ = np.linspace(mesh.geometry.x[:, 0].min(), mesh.geometry.x[:, 0].max(), Nx+1)
y_ = np.linspace(mesh.geometry.x[:, 1].min(), mesh.geometry.x[:, 1].max(), Ny+1)
X_, Y_ = np.meshgrid(x_, y_)
return X_, Y_
def fenics_2_numpy(var):
X_,Y_ = xy_grid(var.function_space.mesh)
points, cell_indices = cells_at_xy(X_,Y_, var.function_space.mesh,True)
numpy_var = var.eval(points, cell_indices)
numpy_var = numpy_var[:, 0].reshape(X_.shape)
return numpy_var
# Solver functions ############################################################################################
def solve_poisson_green(f, h, pad = True):
Ny, Nx = f.shape
if pad: # 2N x 2N
dist_y = np.linspace(-(Ny-1)*h, (Ny-1)*h, 2*Ny-1)
dist_x = np.linspace(-(Nx-1)*h, (Nx-1)*h, 2*Nx-1)
else: # N x N
dist_y = np.linspace(0, (Ny-1) * h, Ny-1)
dist_x = np.linspace(0, (Nx-1) * h, Nx-1)
DX, DY = np.meshgrid(dist_x, dist_y)
R = np.sqrt(DX**2 + DY**2)
with np.errstate(divide='ignore'):
G = -1/(2*np.pi) * np.log(R)
G[np.where(dist_y == 0)[0], np.where(dist_x == 0)[0]] = -1/(2*np.pi) * np.log(h) + 0.16889 # Singularity correction: cell average
chi = fftconvolve(f, G, mode='same') * h**2
return chi
def solve_poisson_fft(f, h):
f_int = f[1:-1,1:-1] # Truncate off the values at x,y = 0 & L
m = f_int.shape[0]
Ne = (m+1)*2
# Boundary Conditions
u0y = 0;
u1y = 0;
ux0 = 0;
ux1 = 0;
# Apply BCs
f[:,0] = f[:,0] - np.transpose(ux0) /h**2;
f[:,m-1] = f[:,m-1] - np.transpose(ux1) /h**2;
f[0,:] = f[0,:] - u0y/h**2;
f[m-1,:] = f[m-1,:] - u1y/h**2;
# Odd extension
g = np.hstack((np.zeros((m,1)), f_int, np.zeros((m,1)), -f_int[:,::-1])) # Odd extension of f in x
g = np.vstack((np.zeros((1,Ne)), g, np.zeros((1,Ne)), -g[::-1,:])) # Odd extension of g in y
ghat = np.fft.fft2(g);
# wavenumbers/eigenvalues
gridx, gridy = np.arange(0,Ne), np.arange(0,Ne)
I, J = np.meshgrid(gridx,gridy) # New finite indices after extension
omega_x, omega_y = 2*np.pi/Ne * I, 2*np.pi/Ne * J # Convert finite space into Fourier Space:
# omega_indice = 2pi/(number of elements along dimension) * indice
mu = -(4/h**2)*(np.sin(omega_x/2)**2 + np.sin(omega_y/2)**2)
mu[0,0] = 1
v = np.real(np.fft.ifft2(ghat/mu))
VP = v[:(m+2),:(m+2)] # Extract u from its odd extension
return VP
def solve_poisson_DSTI(f):
'''
This Process prevents the need to perform an odd extension
by instead using an Discrete Sine Transform technique.
'''
L = x.max()
f_int = f[1:-1,1:-1]# Truncate off the values on the boundary where x,y = 0 & L
m = f_int.shape[0]
fhat = dst(dst(f_int, axis = 0, type = 1, norm = 'ortho'), axis = 1, type = 1, norm = 'ortho')
# Wavenumbers/eigenvalues
omega_x, omega_y = np.arange(1, m+1, 1) * np.pi/L, np.arange(1, m+1, 1) * np.pi/L # kx = pi x / L, ky = pi y / L
Omega_X, Omega_Y = np.meshgrid(omega_x, omega_y, indexing= 'ij')
mu = Omega_X**2 + Omega_Y**2
uhat = -fhat/mu
VP_int = idst(idst(uhat, axis = 0, type = 1, norm = 'ortho'), axis = 1, type = 1, norm = 'ortho')
# Pad result with Dirichlet boundary values
VP = np.pad(VP_int, pad_width = 1, mode = 'constant', constant_values = 0)
return VP
def solve_poisson_fenics(f, inf_space = False, NAME = "uh"):
'''
Solves the Poisson Equation using the standard LU solver and Lagrangian functionspace.
For inf_space = True, from the Electromagnetics example,
"Since we cannot solve the problem on an infinite domain, we will truncate the domain"
So we embed the domain of interest inside of a larger square and set the solution to 0 on the boundary of the large mesh
to simulate a solution going to 0 as space goes to 0.
'''
family = "Lagrange"
order = 1
def scipy_to_fenics_fun(x):
'''
moves from our scipy interpolated function
to a Dolfinx function
'''
# To return a scipy function value, we need to pass in a tuple (x,y)
point = (x[0], x[1])
return f(point)
if inf_space:
gmsh.initialize()
#gmsh.option.setNumber("General.Terminal", 0)# prevent output
# ---- Create Graded mesh ----
L = 1 # Length of actual domain
L_outer = 5 # Length of outer domain (where space goes to "infinity")
coarse_h = h*10
fine_h = h
# lower-left corner coordinates of full and inner domain
x_i, y_i, z_i = x.min(), y.min(), 0 # (0,0,0)
x_f, y_f, z_f = x_i-(L_outer-L)/2, y_i-(L_outer-L)/2, 0 # (-2.5,-2.5,0)
full_square = gmsh.model.occ.add_rectangle(x_f, y_f, z_f, dx = L_outer, dy = L_outer) # full domain
inner_square = gmsh.model.occ.add_rectangle(x_i, y_i, z_i, dx = L, dy = L) # inner domain
# overlap the squares
gmsh.model.occ.fragment([(2, full_square)], [(2, inner_square)])
gmsh.model.occ.synchronize()
surfaces = gmsh.model.getEntities(dim=2)
surface_ids = [s[1] for s in surfaces]
surface_tag = 1
# Tag physical entities
gmsh.model.addPhysicalGroup(2, surface_ids, tag=surface_tag)
gmsh.model.setPhysicalName(2, surface_tag, "msh_domain")
gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "FacesList", [inner_square])
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "InField", 1)
gmsh.model.mesh.field.setNumber(2, "SizeMin", fine_h)
gmsh.model.mesh.field.setNumber(2, "SizeMax", coarse_h)
gmsh.model.mesh.field.setNumber(2, "DistMin", 0) # Grid starts coarsening immediately outisde of inner square
gmsh.model.mesh.field.setNumber(2, "DistMax", L) # Grid hits target coarseness resolution at L distance away from boundaries
gmsh.model.mesh.field.setAsBackgroundMesh(2)
gmsh.model.mesh.generate(2)
# Create dolfinx mesh
model_rank = 0
mesh = df.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, rank=0, gdim=2)[0]
gmsh.finalize()
prefix = "infinite_space"
else:
mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, N-1,N-1)
prefix = "standard"
tdim = mesh.topology.dim
fdim = tdim-1
mesh.topology.create_connectivity(fdim, tdim)
S = df.fem.functionspace(mesh, (family, order)) # Scalar space to store divergence
# boundary condition
boundary_facets = df.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = df.fem.locate_dofs_topological(S,fdim, boundary_facets)
bc = df.fem.dirichletbc(df.default_scalar_type(0),boundary_dofs,S)
# forcing term
forcing = df.fem.Function(S)
forcing.interpolate(scipy_to_fenics_fun)
# Solve
u, v = ufl.TrialFunction(S), ufl.TestFunction(S)
a = ufl.dot(ufl.grad(u),ufl.grad(v)) * ufl.dx
L = -forcing*v*ufl.dx
problem = LinearProblem(a=a,L=L,bcs = [bc],
petsc_options_prefix=f"{prefix}_Poisson",petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
return uh
USE_DATA = False
PLOT_ONLY_GREEN = True
if USE_DATA:
file = "Heavy_rainfall_north_China_20120712_850mb_1800Z.grib"
ds = xr.open_dataset(file,engine = "cfgrib")
res = 1
geopotential = ds.variables['z'][::res,::res].values
gpot_height = geopotential/9.8
lats = ds.variables['latitude'][::res].values
lons = ds.variables['longitude'][:len(lats)].values
div = ds.variables['d'][::res,:len(lats)].values
center_lat = np.deg2rad(lats[(len(lats)-1)//2])
center_lon = np.deg2rad(lons[(len(lons)-1)//2])
LONS, LATS = np.meshgrid(lons,lats)
# ---- map scale factor ----
# define map scale factor
radlats = np.deg2rad(LATS)
n = np.sin(center_lat)
m = np.cos(center_lat)/np.cos(radlats) * (np.tan(np.pi/4 - radlats/2)/np.tan(np.pi/4 - center_lat/2))**n
# scale forcing term
D = div/m**2
N = D.shape[0]
L = 1
h = L/(N-1) # dx = dy = 1/60
x = np.arange(0,h*N,h)
y = np.arange(0,h*N,h)
X,Y = np.meshgrid(x,y)
else:
N = 60
L = 1
h = L/(N-1)
x = np.arange(0,h*N,h)
y = np.arange(0,h*N,h)
X,Y = np.meshgrid(x,y)
a = 50
D = 4 * a * (a * ((X-1/2)**2 + (Y-1/2)**2) - 1) * np.exp(-a*((X-1/2)**2+(Y-1/2)**2))
u_ex = np.exp(-a*((X-1/2)**2 + (Y-1/2)**2))
D_interp = RegularGridInterpolator(
(y, x),
D,
bounds_error=False, fill_value=0.0)
# Collect solutions
u_green = solve_poisson_green(-D, h, pad = True)
u_fft = solve_poisson_fft(D, h)
u_dst = solve_poisson_DSTI(D)
u_LU = solve_poisson_fenics(D_interp, False)
u_LU_inf_space = solve_poisson_fenics(D_interp, True, NAME = "uh_inf")# Create functional expressions to put into Dolfinx
# extract infinite space solution over desired domain
x_ext = np.arange(u_LU_inf_space.function_space.mesh.geometry.x[:,0].min(),u_LU_inf_space.function_space.mesh.geometry.x[:,0].max()+h,h)
y_ext = np.arange(u_LU_inf_space.function_space.mesh.geometry.x[:,1].min(),u_LU_inf_space.function_space.mesh.geometry.x[:,1].max()+h,h)
left_bound = np.where(np.isclose(x_ext,x.min()))[0][0]
right_bound = np.where(np.isclose(x_ext,x.max()))[0][0]
lower_bound = np.where(np.isclose(y_ext,y.min()))[0][0]
upper_bound = np.where(np.isclose(y_ext,y.max()))[0][0]
u_LU_inf_space = fenics_2_numpy(u_LU_inf_space).T[lower_bound:upper_bound+1,left_bound:right_bound+1]
fig = plt.figure(figsize=(12,15))
if PLOT_ONLY_GREEN:
ax1,ax2,ax3 = fig.add_subplot(321),fig.add_subplot(322),fig.add_subplot(323)
if USE_DATA:
D = D[::-1,::-1]
u_green = u_green[::-1,::-1]
u_LU_inf_space = u_LU_inf_space[::-1,::-1]
contour_f = ax1.contourf(X,Y,D, levels = 11, cmap = 'bwr', vmin = -abs(D).max(), vmax = abs(D).max())
ax1.set_title(r"divergence s^{-1}")
ax2.contour(X,Y,D, levels = [-6e-4,-4.5e-4,-3e-4,-1.5e-4,-1e-4,1e-4,1.5e-4,2e-4,2.5e-4,3e-4,3.5e-4], cmap = 'bwr', vmin = -6e-4, vmax = 6e-4, linestyles = '-', linewidth = 1.5)
ax3.contour(X,Y,D, levels = [-6e-4,-4.5e-4,-3e-4,-1.5e-4,-1e-4,1e-4,1.5e-4,2e-4,2.5e-4,3e-4,3.5e-4], cmap = 'bwr', vmin = -6e-4, vmax = 6e-4, linestyles = '-', linewidth = 1.5)
plt.colorbar(contour_f)
else:
contour_u = ax1.contourf(X,Y,u_ex, levels = 11, cmap = 'bwr', vmin = -abs(u_ex).max(), vmax = abs(u_ex).max())
ax1.set_title("exact solution")
plt.colorbar(contour_u)
contour_green = ax2.contourf(X,Y,u_green, levels = 11, cmap = 'bwr', vmin = -abs(u_green).max(), vmax = abs(u_green).max())
contour_LU_inf = ax3.contourf(X,Y,u_LU_inf_space, levels = 11, cmap = 'bwr', vmin = -abs(u_LU_inf_space).max(), vmax = abs(u_LU_inf_space).max())
ax1.set_xlabel('x');ax1.set_ylabel('y')
ax2.set_xlabel('x');ax2.set_ylabel('y')
ax2.set_title("inf. space green sol.")
ax3.set_xlabel('x');ax3.set_ylabel('y')
ax3.set_title("LU inf. space sol.")
else:
# make u_LU plottable in matplotlib
u_LU = fenics_2_numpy(u_LU).T
fig = plt.figure(figsize=(12,15))
ax1,ax2,ax3,ax4,ax5,ax6 = fig.add_subplot(321),fig.add_subplot(322),fig.add_subplot(323),fig.add_subplot(324),fig.add_subplot(325),fig.add_subplot(326)
if USE_DATA:
D = D[::-1,::-1]
u_dst = u_dst[::-1,::-1]
u_fft = u_fft[::-1,::-1]
u_LU = u_LU[::-1,::-1]
u_green = u_green[::-1,::-1]
u_LU_inf_space = u_LU_inf_space[::-1,::-1]
contour_f = ax1.contourf(X,Y,D, levels = 11, cmap = 'bwr', vmin = -abs(D).max(), vmax = abs(D).max())
ax2.contour(X,Y,D, levels = [-6e-4,-4.5e-4,-3e-4,-1.5e-4,-1e-4,1e-4,1.5e-4,2e-4,2.5e-4,3e-4,3.5e-4], cmap = 'bwr', vmin = -6e-4, vmax = 6e-4, linestyles = '-', linewidth = 1.5)
ax3.contour(X,Y,D, levels = [-6e-4,-4.5e-4,-3e-4,-1.5e-4,-1e-4,1e-4,1.5e-4,2e-4,2.5e-4,3e-4,3.5e-4], cmap = 'bwr', vmin = -6e-4, vmax = 6e-4, linestyles = '-', linewidth = 1.5)
ax4.contour(X,Y,D, levels = [-6e-4,-4.5e-4,-3e-4,-1.5e-4,-1e-4,1e-4,1.5e-4,2e-4,2.5e-4,3e-4,3.5e-4], cmap = 'bwr', vmin = -6e-4, vmax = 6e-4, linestyles = '-', linewidth = 1.5)
ax5.contour(X,Y,D, levels = [-6e-4,-4.5e-4,-3e-4,-1.5e-4,-1e-4,1e-4,1.5e-4,2e-4,2.5e-4,3e-4,3.5e-4], cmap = 'bwr', vmin = -6e-4, vmax = 6e-4, linestyles = '-', linewidth = 1.5)
ax6.contour(X,Y,D, levels = [-6e-4,-4.5e-4,-3e-4,-1.5e-4,-1e-4,1e-4,1.5e-4,2e-4,2.5e-4,3e-4,3.5e-4], cmap = 'bwr', vmin = -6e-4, vmax = 6e-4, linestyles = '-', linewidth = 1.5)
ax1.set_title(r"divergence s^{-1}")
plt.colorbar(contour_f)
else:
contour_u = ax1.contourf(X,Y,u_ex, levels = 11, cmap = 'bwr', vmin = -abs(u_ex).max(), vmax = abs(u_ex).max())
plt.colorbar(contour_u)
contour_dst = ax2.contourf(X,Y,u_dst, levels = 11, cmap = 'bwr', vmin = -abs(u_dst).max(), vmax = abs(u_dst).max())
contour_fft = ax3.contourf(X,Y, u_fft, levels = 11, cmap = 'bwr', vmin = -abs(u_fft).max(), vmax = abs(u_fft).max())
contour_LU = ax4.contourf(X,Y,u_LU, levels = 11, cmap = 'bwr', vmin = -abs(u_LU).max(), vmax = abs(u_LU).max())
contour_green = ax5.contourf(X,Y,u_green, levels = 11, cmap = 'bwr', vmin = -abs(u_green).max(), vmax = abs(u_green).max())
contour_LU_inf = ax6.contourf(X,Y,u_LU_inf_space, levels = 11, cmap = 'bwr', vmin = -abs(u_LU_inf_space).max(), vmax = abs(u_LU_inf_space).max())
ax1.set_xlabel('x');ax1.set_ylabel('y')
ax2.set_xlabel('x');ax2.set_ylabel('y')
ax2.set_title("DST-I solution")
ax3.set_xlabel('x');ax3.set_ylabel('y')
ax3.set_title("FFT solution")
ax4.set_xlabel('x');ax4.set_ylabel('y')
ax4.set_title("LU Solution")
ax5.set_xlabel('x');ax5.set_ylabel('y')
ax5.set_title("inf. space green sol.")
ax6.set_xlabel('x');ax6.set_ylabel('y')
ax6.set_title("LU inf. space sol.")
plt.colorbar(contour_dst)
plt.colorbar(contour_fft)
plt.colorbar(contour_LU)
plt.colorbar(contour_green)
plt.colorbar(contour_LU_inf)
plt.show()
error_L2(u_ex,u_green)
error_L2(u_ex,u_LU_inf_space)