Additional parameters for 3d Poisson's equation for electrostatic potential computation

[troyrock]

I'm trying to model a the electrostatic potential in a 3D volume with 3 electrodes on the boundary. I would like to include the voltage of the electrodes as parameters to the model so that I can quickly produce the potential in the volume for a variety of different values of voltages. Is there a way to specify the range over which I would like to have the model valid for? In this example, I would like to try to develop a model for values between -1 an 1. The input values to the neural network are: x, y, z, v_left, v_middle, v_right.

Here are some code snippets:

       def pde(x, y):
           dy_xx = dde.grad.hessian(y, x, i=0, j=0)
           dy_yy = dde.grad.hessian(y, x, i=1, j=1)
           dz_zz = dde.grad.hessian(y, x, i=2, j=2)
           return epsilon*(dy_xx + dy_yy + dz_zz) - sigma

    def top_bottom_boundary_condition(x):  # Top and bottom boundary conditions are Dirchlet
           result = []
           for coords in x:
               if coords[2] == 200:
                   # left circle
                   if (coords[0]-ldotx)**2 + (coords[1] - ldoty)**2 < ldotr**2:
                       result.append(coords[3])  # voltage of left electrode       << this doesn't work because the voltage is not a space coordinate, how do I access the neural network input?
                   # right circle
                   elif (coords[0]-rdotx)**2 + (coords[1] - rdoty)**2 < rdotr**2:
                       result.append(coords[5])  # voltage of the right electrode
   
                   # middle rectangle
                   elif coords[0] < rect_width/2.0 and coords[1] < rect_height/2.0:
                       result.append(coords[4])  # voltage of the middle electrode
   
                   else:
                       result.append(0)  # zero everywhere except on the electrodes.
   
               else:
                   result.append(0)  # zero at the bottom
   
               return np.array(result, ndmin=2).transpose()
 
     data = dde.data.PDE(geom, pde, [bc_top_bottom, bc_side], num_domain=1200, num_boundary=120, num_tes  t=1500)
     net = dde.maps.FNN([6] + [50] * 4 + [1], "tanh", "Glorot uniform")

Thank you.

[lululxvi]

You can specify the range when you define the geometry https://deepxde.readthedocs.io/en/latest/modules/deepxde.geometry.html#deepxde.geometry.geometry_nd.Hypercube

[troyrock]

Thank you! I appreciate your taking time to help me. Now I understand that the geometry can refer to more than the physical geometry of the PDE. I have modified my program to try to implement this but I'm getting an error which indicates that I don't have it quite right:

Traceback (most recent call last):
File "poisson3d.py", line 164, in
main()
File "poisson3d.py", line 139, in main
data = dde.data.PDE(geom,
File "/home/troyrock/anaconda3/envs/deepxde/lib/python3.9/site-packages/deepxde/data/pde.py", line 88, in init
self.train_next_batch()
File "/home/troyrock/anaconda3/envs/deepxde/lib/python3.9/site-packages/deepxde/utils.py", line 24, in wrapper
return func(self, *args, **kwargs)
File "/home/troyrock/anaconda3/envs/deepxde/lib/python3.9/site-packages/deepxde/data/pde.py", line 144, in train_next_batch
self.train_x_all = self.train_points()
File "/home/troyrock/anaconda3/envs/deepxde/lib/python3.9/site-packages/deepxde/data/pde.py", line 202, in train_points
tmp = self.geom.random_boundary_points(
File "/home/troyrock/anaconda3/envs/deepxde/lib/python3.9/site-packages/deepxde/geometry/geometry_3d.py", line 26, in random_boundary_points
x_corner = np.vstack(
File "<array_function internals>", line 5, in vstack
File "/home/troyrock/anaconda3/envs/deepxde/lib/python3.9/site-packages/numpy/core/shape_base.py", line 283, in vstack
return _nx.concatenate(arrs, 0)
File "<array_function internals>", line 5, in concatenate
ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 6 and the array at index 1 has size 3

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
import deepxde as dde
from deepxde.backend import tf

# geometry variables
ldotx = -45
ldoty = 0
ldotr = 50
min_l = -0.1
max_l = 0.1

rect_width = 30
rect_height = 200
min_m = -0.1
max_m = 0.1

rdotx = 45
rdoty = 0
rdotr = 50
min_r = -0.1
max_r = 0.1

bottom_z = -200
top_z = 200
left_x = -200
right_x = 200
back_y = -200
front_y = 200
sheet_z = 200 - 40  # sheet of charge 40 nm from the top
quantum_top = 200 - 134.3  # z coordinate for the top of the quantum layer
quantum_bottom = 200 - 156.8  # z coordinate for the bottom of the quantum layer
lower_left_corner = [left_x, back_y, bottom_z, min_l, min_m, min_r]
upper_right_corner = [right_x, front_y, top_z, max_l, max_m, max_r]

# material variables
sigma0 = 0.0036
epsilon_above = 0.769
epsilon_quantum = 0.697
epsilon_below = 0.769

def main():
    def boundary(_, on_boundary):
        return on_boundary

    def boundary_t_b(x, on_boundary):
        return on_boundary and np.isclose(abs(x[2]), upper_right_corner[2])

    def top_bottom_boundary_condition(x):
        # Top and bottom boundary conditions are Dirchlet
        result = []
        for coords in x:
            if coords[2] == 200:
                # left circle
                if (coords[0]-ldotx)**2 + (coords[1] - ldoty)**2 < ldotr**2:
                    result.append(net.inputs[3])

                # middle rectangle
                elif coords[0] < rect_width/2.0 and coords[1] < rect_height/2.0:
                    result.append(net.inputs[4])

                # right circle
                elif (coords[0]-rdotx)**2 + (coords[1] - rdoty)**2 < rdotr**2:
                    result.append(net.inputs[5])

                else:
                    result.append(0)  # zero everywhere except on the electrodes.

            else:
                result.append(0)  # zero at the bottom

            return np.array(result, ndmin=2).transpose()

    def boundary_side(x, on_boundary):
        return on_boundary and not np.isclose(abs(x[2]), upper_right_corner[2])

    above_sheet = dde.geometry.Cuboid(xmin=[left_x, back_y, sheet_z, min_l, min_m, min_r],
                                      xmax=upper_right_corner)
    between = dde.geometry.Cuboid(xmin=[left_x, back_y, quantum_top, min_l, min_m, min_r],
                                  xmax=[right_x, front_y, sheet_z, max_l, max_m, max_r])
    quantum = dde.geometry.Cuboid(xmin=[left_x, back_y, quantum_bottom, min_l, min_m, min_r],
                                  xmax=[right_x, front_y, quantum_top, max_l, max_m, max_r])
    below = dde.geometry.Cuboid(xmin=lower_left_corner,
                                xmax=[right_x, front_y, quantum_bottom, max_l, max_m, max_r])

    def boundary_above_sheet(x, on_boundary):
        # above sheet of charge 40 nm down from the top
        return on_boundary and above_sheet.inside(x)

    def boundary_between_sheet_and_quantum_layer(x, on_boundary):
        return on_boundary and between.inside(x)

    def boundary_in_quantum_layer(x, on_boundary):
        return on_boundary and quantum.inside(x)

    def boundary_below_quantum_layer(x, on_boundary):
        return on_boundary and below.inside(x)

    def pde_above(x, y):
        dy_xx = dde.grad.hessian(y, x, i=0, j=0)
        dy_yy = dde.grad.hessian(y, x, i=1, j=1)
        dz_zz = dde.grad.hessian(y, x, i=2, j=2)
        return epsilon_above*(dy_xx + dy_yy + dz_zz) - sigma0*x[:, 2:3]

    def pde_between(x, y):
        dy_xx = dde.grad.hessian(y, x, i=0, j=0)
        dy_yy = dde.grad.hessian(y, x, i=1, j=1)
        dz_zz = dde.grad.hessian(y, x, i=2, j=2)
        return epsilon_above*(dy_xx + dy_yy + dz_zz) - sigma0*sheet_z

    def pde_quantum(x, y):
        dy_xx = dde.grad.hessian(y, x, i=0, j=0)
        dy_yy = dde.grad.hessian(y, x, i=1, j=1)
        dz_zz = dde.grad.hessian(y, x, i=2, j=2)
        return epsilon_quantum*(dy_xx + dy_yy + dz_zz) - sigma0*sheet_z

    def pde_below(x, y):
        dy_xx = dde.grad.hessian(y, x, i=0, j=0)
        dy_yy = dde.grad.hessian(y, x, i=1, j=1)
        dz_zz = dde.grad.hessian(y, x, i=2, j=2)
        return epsilon_below*(dy_xx + dy_yy + dz_zz) - sigma0*sheet_z

    geom = dde.geometry.Cuboid(xmin=lower_left_corner, xmax=upper_right_corner)
    bc_top_bottom = dde.DirichletBC(geom, top_bottom_boundary_condition, boundary_t_b)
    bc_side = dde.PeriodicBC(geom, 0, boundary_side)
    condition_above = dde.OperatorBC(geom, lambda x, y, _: pde_above, boundary_above_sheet)
    condition_between = dde.OperatorBC(geom, lambda x, y, _: pde_between, boundary_between_sheet_and_quantum_layer)
    condition_quantum = dde.OperatorBC(geom, lambda x, y, _: pde_quantum, boundary_in_quantum_layer)
    condition_below = dde.OperatorBC(geom, lambda x, y, _: pde_below, boundary_below_quantum_layer)
    data = dde.data.PDE(geom,
                        None,
                        [bc_top_bottom,
                         bc_side,
                         condition_above,
                         condition_between,
                         condition_quantum,
                         condition_below],
                        num_domain=1200,
                        num_boundary=120,
                        num_test=1500)

    net = dde.maps.FNN([6] + [50] * 4 + [1], "tanh", "Glorot uniform")
    model = dde.Model(data, net)

    model.compile("adam", lr=0.001)
    model.train(epochs=50000, model_save_path="project_model")
    # model.train(epochs=50000)
    model.compile("L-BFGS-B")
    losshistory, train_state = model.train()
    dde.saveplot(losshistory, train_state, issave=True, isplot=True)
    model._print_model()


if __name__ == "__main__":
    main()

[lululxvi]

You should use Hypercube

[troyrock]

NotImplementedError: Hypercube.random_boundary_points to be implemented

[lululxvi]

Ah, yeah, the sampling on the boundary of a hypercube has not been implemented yet. You can set 'num_boundary=0' and manually generate the points on the boundary and then pass them using 'anchors'.

[harshil-patel-code]

@troyrock Can you please post your working code? I want to learn how to use Hypercube and sample boundary points as suggested by @lululxvi

Thank You :)

[troyrock]

@harshil-patel-code I wasn't able to get this working and ended up using Julia's SciML suite of tools instead. The above code is basically as far as I got. Sorry.

[harshil-patel-code]

Thanks @troyrock for your reply.

Hi @lululxvi, I have developed the missing random_boundary_points function for Hypercube. I have tested it on my problem. It works fine. Feel free to include it in the repo. If you like, I can send merge request.

def random_boundary_points(self, n, random="pseudo"):
        x = sample(n, self.dim, random)
        rand_dim = np.floor(sample(n, 1, random)*self.dim).astype('int') # randomly picking dimension
        x[np.arange(n), rand_dim[:,0]] = np.round(x[np.arange(n), rand_dim[:,0]]) # replacing value of the randomly picked dimension with the nearest boundary value (0 or 1)
        return (self.xmax - self.xmin) * x + self.xmin

[lululxvi]

@harshil-patel-code Great. Feel free to send a pull request.