How to solve a 2D problem with an obstacle with some interior nodes inside of it

[RomanLF]

Hello again, I meet difficulties to procceed diffusion on a domain occupied by a sort of court-yard, with the model :

\begin{align}
- v\Delta u &= s\delta(x-x_S), \quad &x \in \Omega \\
\nabla u.n &= 0 \quad &x \in \Gamma_N \\
u(x) &= 0, \quad &x \in \Gamma_D
\end{align}

on the domain :

The source is located at the outside of the courtyard. A neuman zero-flux condition is set on the edges of the obstacle. However, the solution inside the courtyard is not equal to zero which is not physical as visible below.

Do you know how I could avoid this issue ?

Here is my code :

# -*- coding: utf-8 -*-
"""
Created on 12 2021

@author: Roman Lopez-Ferber
"""
import numpy as np

from rbf.pde.geometry import simplex_outward_normals
from rbf.pde.fd import weight_matrix, weights
from rbf.sputils import expand_rows
from scipy.spatial import cKDTree
import time

import matplotlib
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.nodes import min_energy_nodes
from rbf.pde.geometry import contains

######################
###   Geometry    ###
######################

# vertices and simplices for a "courtyard" geometry.
# Five first vertices : Exterior domain

vertices = np.array([[  0.,   0.],
                     [125.,   0.],
                     [125., 100.],
                     [  0., 100.],
                     [  0.,   0.],
                     [ 40.,  30.],
                     [ 80.,  30.],
                     [ 80.,  70.],
                     [ 40.,  70.],
                     [ 40.,  30.],
                     [ 45.,  35.],
                     [ 75.,  35.],
                     [ 75.,  65.],
                     [ 45.,  65.],
                     [ 45.,  35.]])

simplices = np.array([[ 0,  1],
                     [ 1,  2],
                     [ 2,  3],
                     [ 3,  4],
                     [ 5,  6],
                     [ 6,  7],
                     [ 7,  8],
                     [ 8,  9],
                     [10, 11],
                     [11, 12],
                     [12, 13],
                     [13, 14]])

ext_simplices = np.array([[ 0,  1],
                          [ 1,  2],
                          [ 2,  3],
                          [ 3,  4]])

# 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)),
    }

##### Domain building ###########

N_nodes = 6000
nodes, groups, normals = min_energy_nodes(N_nodes, 
                                          (vertices, simplices),
                                          boundary_groups=boundary_groups,
                                          boundary_groups_with_ghosts=['interior'])

### diffusion parameters ###
# Diffusion
Diff_x = 100 # m^2/s
Diff_y = 100 # m^2/s

# Source term parameters
source_idx = 887
x_s_n, y_s_n = 130.23213463982668, 150.413853115689
As = 1e3 # amplitude source

pos_source = source_idx

########## RBF simulation ###########

phi = 'phs5'
stencil_size = 45
order = 7

print("Solving ...")
tic = time.perf_counter()

# Enforce that the Lapacian is 0 at the interior nodes (not to be confused
# with the nodes on the interior boundary)
D_laplacian = weight_matrix(
    nodes[groups['interior']],
    nodes,
    stencil_size,
    [[2, 0], [0, 2]],
        coeffs=[
        Diff_x,
        Diff_y
        ],
    phi=phi,
    order=order
    )

# source term marked by a cross in graph
b_interior = np.full(len(groups['interior']), 0)
b_interior[pos_source] = As

# Enforce that u = 0 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 =  np.full(len(groups['boundary:exterior']), 0) 

# Enforce a free boundary at the interior boundary nodes. NEUMANN : $\nabla u \cdot n = 0$
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 0 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']), 0)

# Combine the LHS and RHS components and solve
A = sp.vstack(
    (- D_laplacian, 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)

toc = time.perf_counter()

print(f"ellapsed time = {toc - tic} s")

#%% plot solution
# The rest is just plotting...

# %matplotlib qt

import matplotlib
import matplotlib.colors as colors

fig, ax = plt.subplots(figsize=(8, 8))

fontsize = 12
font = {'size'   : fontsize}
matplotlib.rc('font', **font)

for smp in simplices:
    ax.plot(vertices[smp, 0], vertices[smp, 1], 'k-')

p = ax.scatter(nodes[groups['interior']][:, 0], nodes[groups['interior']][:, 1],
               s=10,
               c=soln[groups['interior']], cmap='jet',
               norm=colors.LogNorm(vmin= 1e-1, vmax = int(np.max(1.1*soln))))

fig.colorbar(p)
ax.grid(ls=':')
ax.set_aspect('equal')

x_s, y_s = nodes[groups['interior']][pos_source]
ax.plot(x_s, y_s, marker = 'x', markersize=12, linewidth = 20, label = "source")
ax.legend(loc = 'upper right')

nb_nodes = 'N =' + str(N_nodes)
ax.set_title('           Solution at nodes' + ' $' + nb_nodes + '$ \n'
                          + '$\phi $ = ' f'{phi}, order = {order},\n'
                           f'stencil_size = {stencil_size}')

fig.show()

#%% plot domain

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(loc = 'upper right')
ax1.set_title('Node groups')

fig1.show()

Best regards, Roman

[treverhines]

Unless I am misunderstanding the problem, it seems like the problem is ill-posed since any constant value at the center would satisfy those conditions.

You should also keep in mind that derivatives are computed using the nearest neighboring nodes, regardless of any boundaries in the domain. So derivatives at the center of the domain may be computed using some nodes from outside of that center region, which I am guessing you would not want. I think this is why you are getting any solution at all rather than a singular matrix error.

[RomanLF]

Hello,

Thank you for your answer. As you said, these nodes seem to disturb the solution, as the results are different when the 2D-courtyard is not considered. I imagine that in 3D with a limited wall height the problem would disappear too as court-yard's nodes would be connected to the exterior by vertical axis.

Thank you a lot for your answer !