Flow in porous media - 2D/single phase/Homogeneous - Neumann BC definition

[gaboperezayala]

Hello Dr. Lu, I'm trying to implement a 2D pressure equation to solve single-phase flow in homogeneous porous media in DeepXDE, but I have some doubts about the Neumann boundary conditions definition.

The equation to solve is

p_xx + p_yy = 0 in Ω = [0,1] x [0,1]

with the Dirichlet boundary conditions:

p = 10 on Ω = (x,Ly), x ε [0,1] p = 5 on Ω = (x,0), x ε [0,1]

and, with the Neumann boundary conditions:

u = 0 on ∂Ω = (0,y), y ε [0,1] u = 0 on ∂Ω = (Lx,y), y ε [0,1]

As Neumann boundary conditions are defined as: dy/dn(x) = func(x). In my study the first-order derivative dp/dx (p_x) is set as 0, as well as, dp/dy (p_y) is set as 0. However, dp/dy (p_y) cannot be 0, it has to be different than 0 (otherwise results has nonsense because it means that there is no pressure gradient in y) AND I HAVEN'T FOUND HOW TO IMPLEMENT THIS ON DeepXDE I'd really appreciate the help you can provide.

Implementation

This description goes through the implementation of a solver for the above described Pressure equation step-by-step.

First, import modules required.

import deepxde as dde
from deepxde.backend import tf
import numpy as np
import matplotlib.pyplot as plt
from numpy import ma
from matplotlib import ticker, cm

Define a computational geometry using a rectangle geometry.

geom = dde.geometry.Rectangle([0,0],[1,1])

Next, express the PDE residual of the equation.

def pde(x,y):
    p_xx = dde.grad.hessian(y, x, i=0, j=0)
    p_yy = dde.grad.hessian(y, x, i=1, j=1)
    return p_xx + p_yy

Next, boundary conditions are considered

def boundary_t(x, on_boundary):
    return on_boundary and np.isclose(x[0], 1)

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

def boundary_r(y, on_boundary):
    return on_boundary and np.isclose(y[0], 1)

def boundary_l(y, on_boundary):
    return on_boundary and np.isclose(y[0], 0)

bc_t = dde.DirichletBC(geom, lambda x: 10, boundary_t)
bc_b = dde.DirichletBC(geom, lambda x: 5, boundary_b)
bc_r = dde.NeumannBC(geom, lambda x: 0, boundary_r, component = 0)
bc_l = dde.NeumannBC(geom, lambda x: 0, boundary_l, component = 0)

Next, define the PDE problem as:

data = dde.data.PDE(geom,
                   pde,
                   [bc_t, bc_b, bc_r, bc_l],
                   num_domain = 1024,
                   num_boundary = 256,
                   num_test = 256)

Next, choose the network, here it is used a fully connected neural network of 1 hidden layer with 30 nodes and 30.000 iterations, using "tanh" as activation function and Xavier Initialization.

net = dde.nn.FNN([2] + [30]*1 + [1],
                "tanh",
                "Glorot normal")

Now, build a Model and choose the optimizer ("adam") and learning rate (0.001), then train the model using 30.000 iterations:

model = dde.Model(data,net)
model.compile("adam", lr = 0.001)
losshistory, train_state = model.train(epochs = 30000)

Plot results

dde.saveplot(losshistory, train_state, issave=True, isplot=True)

X_train, y_train, X_test, y_test, best_y, best_ystd = train_state.packed_data()
x = X_test[:, 0]
y = X_test[:, 1]
z = best_y[:, 0]
fig, ax = plt.subplots()
cs = ax.tricontourf(y, x, z, 50, cmap='rainbow')
cbar = fig.colorbar(cs)
ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
plt.show()

As you can see in the previous figure; if there is no gradient of pressure in y, then the results have nonsense. If Neuman bc are not implemented, then the result has sense as in the next figure. Nevertheless, the problem is that Neumann BC has to be implemented.

[lululxvi]

I am confused. In the figure, why "there is no gradient of pressure in y"?

[gaboperezayala]

Dr. Lulu thanks for your answer the problem was in the definition of the boundary conditions, now it is working as it's supposed, please refer to issue #584