Question about boundary condition with discontinuous derivatives

[forxltk]

Hi, Dr. Lu Lu! I am trying to solve a steady-state seepage through rectangular dam. Here is the definition of the problem with boundary conditions.

Everything is fine except the boundary condition in the right side.
For y<=h, phi is always h, while for y>h, phi=y.
The result is as following:

Although the values of phi are continuous, the derivatives of phi are discontinuous at y=h. And the result is not good.
Do you have any idea for it?
Many thanks!

[lululxvi]

Maybe try a larger network and more points?

[forxltk]

OK, I will follow your suggestion.

[MINE0126]

Hello!Did you solve the problem?

[forxltk]

Hello!Did you solve the problem?

The result improves, but still not perfect.

[MINE0126]

I'm working on a PDE solution recently, how are your boundary conditions set? Can I have a look at the relevant code? Maybe I can help you.

[forxltk]

@MINE0126 Of course. Here is the code.

import matplotlib.pyplot as plt
import numpy as np
import deepxde as dde

import tensorflow as tf
from tensorflow.keras import layers, models


def pde(x, y):
    dpxx = dde.grad.hessian(y, x, i=0, j=0, component=0)
    dpyy = dde.grad.hessian(y, x, i=1, j=1, component=0)
    main_pde = dpxx+dpyy

    s = y[:, 1:2]
    dpn = y[:, 2:3]
    phis = y[:, 3:4]

    return [main_pde, dpn, tf.abs(s-phis)]


def boundary_left(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0.0)


def boundary_right(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0.5)


def boundary_bot(x, on_boundary):
    return on_boundary and np.isclose(x[1], 0.0)


def func_right(x):
    X = x[:, 0:1]
    Y = x[:, 1:2]
    h = tf.math.less_equal(Y, 0.5)
    f1 = tf.ones_like(X) * 0.5
    f2 = tf.ones_like(X) * Y
    f3 = tf.where(h, f1, f2)
    return f3


def func_bot(x, y, _):
    return dde.grad.jacobian(y, x, i=0, j=1)


geom = dde.geometry.Rectangle([0.0, 0.0], [0.5, 1.0])

bc_left = dde.DirichletBC(geom, lambda x: 1, boundary_left, component=0)
bc_right = dde.DirichletBC(geom, func_right, boundary_right, component=0)
s_left = dde.DirichletBC(geom, lambda x: 1, boundary_left, component=1)
bc_bot = dde.OperatorBC(geom, func_bot, boundary_bot)

BC = [bc_right, bc_bot, bc_left, s_left]


data = dde.data.PDE(
    geom, pde, BC,
    num_domain=40000,
    num_boundary=5000,
    )

layer_size = [2, 2]
activation = "tanh"
initializer = "Glorot uniform"

net = dde.maps.FNN(layer_size, activation, initializer)


def modify_output(X, Y):
    x, y = X[:, 0:1], X[:, 1:2]
    x_y = tf.concat([x, y], axis=1)

    layer_size = 100

    model_phi = models.Sequential()
    model_phi.add(layers.Dense(layer_size, activation='tanh'))
    model_phi.add(layers.Dense(layer_size, activation='tanh'))
    model_phi.add(layers.Dense(layer_size, activation='tanh'))
    model_phi.add(layers.Dense(layer_size, activation='tanh'))
    model_phi.add(layers.Dense(layer_size, activation='tanh'))
    model_phi.add(layers.Dense(layer_size, activation='tanh'))
    model_phi.add(layers.Dense(layer_size, activation='tanh'))
    model_phi.add(layers.Dense(layer_size, activation='tanh'))
    model_phi.add(layers.Dense(1, activation=None))
    phi_output = model_phi(x_y)

    model_s = models.Sequential()
    model_s.add(layers.Dense(layer_size, activation='tanh'))
    model_s.add(layers.Dense(layer_size, activation='tanh'))
    model_s.add(layers.Dense(layer_size, activation='tanh'))
    model_s.add(layers.Dense(layer_size, activation='tanh'))
    model_s.add(layers.Dense(layer_size, activation='tanh'))
    model_s.add(layers.Dense(layer_size, activation='tanh'))
    model_s.add(layers.Dense(layer_size, activation='tanh'))
    model_s.add(layers.Dense(layer_size, activation='tanh'))
    model_s.add(layers.Dense(1, activation=None))
    s_output = model_s(x)

    s_concat = tf.concat([x, s_output], axis=1)

    fs_output = model_phi(s_concat)
    dpx = dde.grad.jacobian(fs_output, s_concat, j=0)
    dpy = dde.grad.jacobian(fs_output, s_concat, j=1)

    alpha = tf.math.atan(dpx) + np.pi / 2
    nx = tf.math.cos(alpha)
    ny = tf.math.sin(alpha)

    dpn = dpx*nx+dpy*ny

    final_output = tf.concat([phi_output, s_output, dpn, fs_output], axis=1)
    return final_output


net.apply_output_transform(modify_output)

model = dde.Model(data, net)

loss_weights = [1, 1000, 100, 1, 1, 1, 100]
model.compile('adam', lr=1e-3, loss_weights=loss_weights)
model.train(epochs=100000, display_every=1000)
model.compile('adam', lr=1e-4, loss_weights=loss_weights)
checkpointer = dde.callbacks.ModelCheckpoint(
     "model/model.ckpt", verbose=1, save_better_only=True, period=1000)
#model.restore("model.ckpt-122000.ckpt", verbose=1)
model.train(epochs=100000, display_every=1000, callbacks=[checkpointer])


x = np.linspace(0, 0.5, 100)[:, None]
s_star = np.linspace(0, 0.5, 100)[:, None]
x_star = np.concatenate((x, s_star), axis=1)

output = model.predict(x_star)
s_pred = output[:, 1:2]
phis_pred = output[:, 3:4]
exact_x = np.linspace(0, 0.5, 11)[:, None]
exact_y = np.array([1.0, 0.986242, 0.967625, 0.945590, 0.920382, 0.891939, 0.859969, 0.823876, 0.782493, 0.733142, 0.662382])[:, None]
fs = np.concatenate((exact_x, exact_y), axis=1)
phi_pred = model.predict(fs)[:, 0:1]

fig = plt.figure(2, figsize=(5, 5))
plt.plot(fs[:, 0:1], phi_pred, label='Phi', zorder=1)
plt.plot(x_star[:, 0:1], phis_pred, label='Phis', zorder=2)
plt.plot(x_star[:, 0:1], s_pred, '--', label='S', zorder=3)
plt.scatter(exact_x, exact_y, s=30, c='black', marker='o', zorder=4)
plt.xlabel(r'$t$')
plt.ylabel(r'$s(t)$')
plt.title('Moving Boundary')
plt.legend()
plt.show()

x = np.linspace(0.5, 0.5, 100)[:, None]
s_star = np.linspace(0.0, 1.0, 100)[:, None]
x_star = np.concatenate((x, s_star), axis=1)

s1_star = np.linspace(0.5, 0.5, 50)[:, None]
s2_star = np.linspace(0.5, 1.0, 50)[:, None]
s = np.concatenate((s1_star, s2_star), axis=0)
phi = model.predict(x_star)[:, 0:1]
error_s = np.abs(s - phi)

plt.subplot(1, 2, 1)
plt.plot(s_star, phi, label='Predict')
plt.plot(s_star, s, '--', label='Exact')
plt.xlabel(r'$y$')
plt.ylabel(r'$\phi$')
plt.title('Right Boundary')
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(s_star, error_s)
plt.xlabel(r'$t$')
plt.ylabel(r'Point-wise Error')
plt.title('Absolute Error')
plt.yscale('log')
plt.tight_layout()
plt.show()

[chang-change]

@forxltk I don't understand why layer_size = [2, 2], not layer_size = [4]. And does layer_size not need to set the input layer?

[forxltk]

layer_size only sets the size of inputs. The hidden layers and the size of outputs are defined in modify_output.

[chang-change]

Does layer_size = [2, 2] represent the input of two neural networks? Or a neural network representing 4 inputs?

[forxltk]

The first 2 represent the inputs, the other represent the original output. But the outputs have been modified by modify_output. So whether it is [2, 2] or [2, 4] is OK.

[chang-change]

If layer_size = [2] + [32]*6 + [4], can it be considered that only the input layer of the modified neural network is related to [2], and the hidden layer and output layer no longer use [32]*6 + [4] ? Instead, take the newly defined structure in modify_output()?

[MINE0126]

Hello, @forxltk. I apologize for the late reply.
Have you made any progress recently?
I tried to understand the code you posted above, but I can't understand the following code

def pde(x, y):
    dpxx = dde.grad.hessian(y, x, i=0, j=0, component=0)
    dpyy = dde.grad.hessian(y, x, i=1, j=1, component=0)
    main_pde = dpxx+dpyy
    s = y[:, 1:2]
    dpn = y[:, 2:3]
    phis = y[:, 3:4]

    return [main_pde, dpn, tf.abs(s-phis)]

What is the meaning of S 、 dpn、phis respectively?

[forxltk]

• Unfortunately, I no longer pay attention to this problem.
• Sorry for the ambiguous code. s is the free surface/boundary. dpn is the Neumann BC of the free surface. And phis is the value of phi in the free surface.

[MINE0126]

Thank you for your reply.
I have some questions about your code and would like to ask for your answers.
①Y[:,0:1] is PHI, which represents the value of water head in the domain, and the water head at each place in the domain can be known by PHI. Therefore, why do you still need S, PHIS as part of the control equation?
②I think the left and right boundary conditions are not a constant, they should be determined by the function PHI = pgh, where p is the density of water, g is the acceleration of gravity, and h is the height from the water surface.
Do you agree with me? Please correct me if there is an error. Thank you very much!

[forxltk]

1. Because phi should be controlled by the Neumann BC of the free surface. Without s and phis, the fs BC can't be well defined. Then you will have the wrong phi.
2. It depends on your problem. My problem follows this paper https://www.sciencedirect.com/science/article/pii/S0022169422006576

[MINE0126]

Thank you for your answer.