Having trouble to understand if RBF-FD is able to deal with internal boundaries
Closed this issue · 3 comments
Hello,
I do not know if it is the best place to ask a question but I try. I have carefully read the docs but I still do not know if your library able to deal with internal constraints ? Such as in this example from Medusa, a C++ RBF-FD library,
By the way, is it possible to set a constraint to the inner boundary which differs from the outer one with the RBF library ?
Congratulations for your work,
Roman
Hi @RomanLF. This is where I prefer to get questions about this package, so thank you for posting here.
By internal constraints, do you mean constraints on an interior boundary (like the hole in the above domain)? The RBF package has all the tools needed to solve such a problem, although it probably requires more work to set up the problem than what Medusa has. The most relevant tools provided in this package are 1) a function to generate the sparse RBF-FD differentiation matrices and 2) a function to generate nodes in a domain defined in terms of vertices and simplices (i.e., how the vertices are connected). You should be able to generate nodes for any domain as long as the boundary does not intersect itself or have any gaps that would make it ambiguous to tell whether a point is inside the domain. Holes should not cause any issue.
I had a bit of fun tonight and put together a script to get the steady-state solution for the above PDE (i.e., assume du/dt=0). It is possible to solve the time dependent problem; it would just require more code (See the bottom example here for a time dependent solution: https://rbf.readthedocs.io/en/latest/fd.html). Hopefully the below script demonstrates how to solve the type of problem you are interested in.
import logging
import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse as sp
from rbf.pde.nodes import poisson_disc_nodes
from rbf.pde.fd import weight_matrix
from rbf.pde.geometry import contains
logging.basicConfig(level=logging.DEBUG)
spacing = 0.1
phi = 'phs3'
stencil_size = 30
order = 2
# Define exterior boundary vertices and how they are connected
ext_vertices = np.array([
[-2.0104849279161208, 0.9979079497907952],
[-1.003931847968545, 0.9979079497907952],
[-0.9986893840104851, 0.48012552301255274],
[0.6212319790301444, 0.9351464435146448],
[1.2293577981651373, 1.8713389121338917],
[2.2306684141546524, 0.1349372384937242],
[-0.35910878112712963, -1.3661087866108783],
[-1.1507208387942334, -0.0010460251046020552],
[-2, 0.004184100418410441],
[-2.0104849279161208, 0.9979079497907952]
])
ext_simplices = np.column_stack(
(np.arange(0, len(ext_vertices) - 1), np.arange(1, len(ext_vertices)))
)
# Define interior boundary vertices and how they are connected
int_vertices = np.array([
[0.9781840417671206, 0.6262905818998916],
[1.062457189241084, 0.6073737505811252],
[1.123555221159707, 0.5779475685297106],
[1.1741191096440846, 0.5422157760387072],
[1.2183625120679156, 0.48966902237546717],
[1.252071771057501, 0.4182054373934605],
[1.2647127431785954, 0.3362325016788057],
[1.2499649423706516, 0.24795395552456245],
[1.2267898268153115, 0.18069411083561482],
[1.180439595704632, 0.1218417467327857],
[1.123555221159707, 0.0734987333626047],
[1.0392820736857442, 0.04197068116466074],
[0.9634362409591772, 0.03566507072507186],
[0.8960177229800066, 0.03986881101813111],
[0.8222787189402894, 0.06719312292301582],
[0.7632875157085155, 0.1071286557070783],
[0.7190441132846845, 0.1554716690772593],
[0.6895485116687974, 0.2143240331800884],
[0.6705870534871563, 0.28158387786903605],
[0.668480224800307, 0.3467418524114536],
[0.6811211969214015, 0.41610356724693065],
[0.7232577706583831, 0.49597463281505605],
[0.7759284878296095, 0.554826996917885],
[0.8328128623745346, 0.5947625297019474],
[0.8981245516668559, 0.6178831013137731],
[0.9781840417671206, 0.6262905818998916]
])
int_simplices = np.column_stack(
(np.arange(0, len(int_vertices) - 1), np.arange(1, len(int_vertices)))
)
# Combine the exterior and interior boundary. Make sure to update the interior
# simplices to point to the correct vertex indices.
vertices = np.vstack((ext_vertices, int_vertices))
simplices = np.vstack((ext_simplices, int_simplices + len(ext_vertices)))
# Identify which simplices belong to the interior and which belong to the
# exterior
boundary_groups = {
'exterior': np.arange(len(ext_simplices)),
'interior': np.arange(len(ext_simplices), len(simplices)),
}
# Generate the nodes. `groups` identifies which group each node belongs to.
# `normals` contains the normal vectors corresponding to each node (nan if the
# node is not on a boundary). We are giving the interior boundary nodes ghost
# nodes to improve the accuracy of the free boundary constraint.
nodes, groups, normals = poisson_disc_nodes(
spacing,
(vertices, simplices),
boundary_groups=boundary_groups,
boundary_groups_with_ghosts=['interior']
)
# Enforce that the Lapacian is -5 at the interior nodes (not to be confused
# with the nodes on the interior boundary)
A_interior = weight_matrix(
nodes[groups['interior']],
nodes,
stencil_size,
[[2, 0], [0, 2]],
phi=phi,
order=order
)
b_interior = np.full(len(groups['interior']), -5)
# Enforce that u=x at the exterior boundary
A_boundary_exterior = weight_matrix(
nodes[groups['boundary:exterior']],
nodes,
stencil_size,
[0, 0],
phi=phi,
order=order
)
b_boundary_exterior = nodes[groups['boundary:exterior'], 0]
# Enforce a free boundary at the interior boundary nodes.
A_boundary_interior = weight_matrix(
nodes[groups['boundary:interior']],
nodes,
stencil_size,
[[1, 0], [0, 1]],
coeffs=[
normals[groups['boundary:interior'], 0],
normals[groups['boundary:interior'], 1]
],
phi=phi,
order=order
)
b_boundary_interior = np.zeros(len(groups['boundary:interior']))
# Use the ghost nodes to enforce that the Laplacian is -5 at the interior
# boundary nodes.
A_ghosts_interior = weight_matrix(
nodes[groups['boundary:interior']],
nodes,
stencil_size,
[[2, 0], [0, 2]],
phi=phi,
order=order
)
b_ghosts_interior = np.full(len(groups['ghosts:interior']), -5)
# Combine the LHS and RHS components and solve
A = sp.vstack(
(A_interior, A_boundary_exterior, A_boundary_interior, A_ghosts_interior)
)
b = np.hstack(
(b_interior, b_boundary_exterior, b_boundary_interior, b_ghosts_interior)
)
soln = sp.linalg.spsolve(A, b)
# The rest is just plotting...
fig1, ax1 = plt.subplots()
for smp in simplices:
ax1.plot(vertices[smp, 0], vertices[smp, 1], 'k-')
for name, idx in groups.items():
ax1.plot(nodes[idx, 0], nodes[idx, 1], '.', label=name)
ax1.grid(ls=':')
ax1.set_aspect('equal')
ax1.legend()
ax1.set_title('Node groups')
fig2, ax2 = plt.subplots()
for smp in simplices:
ax2.plot(vertices[smp, 0], vertices[smp, 1], 'k-')
p = ax2.scatter(nodes[:, 0], nodes[:, 1], s=10, c=soln, cmap='jet')
fig2.colorbar(p)
ax2.grid(ls=':')
ax2.set_aspect('equal')
ax2.set_title('Solution at nodes')
fig3, ax3 = plt.subplots()
for smp in simplices:
ax3.plot(vertices[smp, 0], vertices[smp, 1], 'k-')
points_grid = np.mgrid[
nodes[:, 0].min():nodes[:, 0].max():500j,
nodes[:, 1].min():nodes[:, 1].max():500j,
]
points_flat = points_grid.reshape(2, -1).T
from scipy.interpolate import LinearNDInterpolator
soln_interp_flat = LinearNDInterpolator(nodes, soln)(points_flat)
soln_interp_flat[~contains(points_flat, vertices, simplices)] = np.nan
soln_interp_grid = soln_interp_flat.reshape(points_grid.shape[1:])
p = ax3.contourf(*points_grid, soln_interp_grid, cmap='jet')
fig3.colorbar(p)
ax3.grid(ls=':')
ax3.set_aspect('equal')
ax3.set_title('Interpolated solution')
plt.show()
Hello,
First I am glad to read that I have posted this question in the right place.
This answers exactly to my question. Sorry for the title, I was asking on the dealing of internal boundary as you have precisely shown here. Thank you a lot for this new example it explains very well how to deal with this type of problem. It would be a great addition to the docs by the way.
Thank you for the link you provide about solving in a instationnary regime. I already saw you used the method of lines for the acoustic problem.
Thank you a lot !
I am glad I was able to help. I added a similar demo script to the bottom of https://rbf.readthedocs.io/en/latest/fd.html