lululxvi / deepxde

A library for scientific machine learning and physics-informed learning
https://deepxde.readthedocs.io
GNU Lesser General Public License v2.1
2.64k stars 739 forks source link

Bloch periodic boundary condition #624

Closed tuzu123 closed 2 years ago

tuzu123 commented 2 years ago

Hi Prof.Lu, I have a little problem with the BC implementation. The Bloch periodic boundary condition is a more general periodic BC, and the 1D mathematical formula shows as follows: $$u(0,t)cos(k)=u(1,t),\ \ \ u_x(0,t)cos(k)=ux(1,t),\ \ x{\in}[0,1]$$ k is wave vector($k{\in}[0,2\pi]$) and k is fixed value in the process of solving the equation. This means $u(0,t)=u(1,t)$ is a special case for Bloch periodic BC and it has been implemented in class PeriodicBC(BC). Could you please point the way to implementing the Bloch boundary condition in deepXDE? Shed tears of gratitude! For the sake of understanding, the complete PDE show follows: Wave Equ: $$u{tt}=c^2u_{xx},\ \ x{\in}[0,1],t{\in}[0,1]$$ IC: $$u(x,0)=cos(k),\ \ u_t(x,0)=sin(k)$$ BC: $$u(0,t)cos(k)=u(1,t),\ \ u_x(0,t)cos(k)=u_x(1,t)$$

athtareq commented 2 years ago

Hello, maybe try something in the spirit of


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

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

geomtime = dde.geometry.Interval(0, 1)
bc_u_0 = dde.PeriodicBC(geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=0, component=0)
bc_u_1 = dde.PeriodicBC(geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=1, component=0)

bc_u_0=bc_u_0/tf.cos(k) 
bc_u_1 =bc_u_1/tf.cos(k) 

data = dde.data.PDE(geomtime, pde, [bc_u_0,bc_u_1], 100, 2, num_test=200)```
tuzu123 commented 2 years ago

After understanding the source code of class PeriodicBC(BC), I try to write the code of class BlochPeriodicBC. But I am not sure whether there are some stupid mistakes.

class BlochPeriodicBC(BC):

    def __init__(self, geom, component_x, on_boundary, derivative_order=0, component=0, phase=0):
        super().__init__(geom, on_boundary, component)
        self.component_x = component_x
        self.derivative_order = derivative_order
        self.phase = phase
        if derivative_order > 1:
            raise NotImplementedError(
                "BlochPeriodicBC only supports derivative_order 0 or 1."
            )

    def collocation_points(self, X):
        X1 = self.filter(X)
        X2 = self.geom.periodic_point(X1, self.component_x)
        return np.vstack((X1, X2))

    def error(self, X, inputs, outputs, beg, end):
        mid = beg + (end - beg) // 2
        if self.derivative_order == 0:
            yleft = outputs[beg:mid, self.component : self.component + 1]
            yright = outputs[mid:end, self.component : self.component + 1]
        else:
            dydx = grad.jacobian(outputs, inputs, i=self.component, j=self.component_x)
            yleft = dydx[beg:mid]
            yright = dydx[mid:end]
        return yleft * tf.cos(self.phase) - yright

Then, BlochBCs define as follows:

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

bc1 = dde.icbc.BlochPeriodicBC(geom, 0, boundary_l,derivative_order=0, component=0, phase=k)

bc2 = dde.icbc.BlochPeriodicBC(geom, 0, boundary_l,derivative_order=1, component=0, phase=k)

Or there are some easier ways to implement this such as using class OperatorBC(BC)?

tuzu123 commented 2 years ago

Hello, maybe try something in the spirit of

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

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

geomtime = dde.geometry.Interval(0, 1)
bc_u_0 = dde.PeriodicBC(geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=0, component=0)
bc_u_1 = dde.PeriodicBC(geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=1, component=0)

bc_u_0=bc_u_0/tf.cos(k) 
bc_u_1 =bc_u_1/tf.cos(k) 

data = dde.data.PDE(geomtime, pde, [bc_u_0,bc_u_1], 100, 2, num_test=200)```

Hi @athtareq Thank you very much for your reply!

I am a little confused about PeriodicBC. In PeriodicBC, there is no need to distinguish the left boundary and right boundary because they are equal to each other. In this case, we can input all boundary points or just left/right boundary points and the self.geom.periodic_point() will convert the left/right points to the right/left points.

In my opinion, the code is

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

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

geomtime = dde.geometry.Interval(0, 1)
bc_u_0 = dde.PeriodicBC(geomtime, 0, boundary_l, derivative_order=0, component=0)
bc_u_1 = dde.PeriodicBC(geomtime, 0, boundary_l, derivative_order=1, component=0)

bc_u_0=bc_u_0/tf.cos(k) 
bc_u_1 =bc_u_1/tf.cos(k) 

data = dde.data.PDE(geomtime, pde, [bc_u_0,bc_u_1], 100, 2, num_test=200)

But the PeriodicBC returns the error which means bc_u_0 = yleft - yright, so bc_u_0 = bc_u_0/tf.cos(k) means bc_u_0 = (yleft - yright)/tf.cos(k). So it seems didn't match the bc $u(0,t)cos(k) - u(1,t) = 0$ or $u_x(0,t)cos(k) - u_x(1,t) = 0$. Hoping for your reply!

Kind regards, Peng

athtareq commented 2 years ago

Hi @athtareq Thank you very much for your reply!

I am a little confused about PeriodicBC. In PeriodicBC, there is no need to distinguish the left boundary and right boundary because they are equal to each other. In this case, we can input all boundary points or just left/right boundary points and the self.geom.periodic_point() will convert the left/right points to the right/left points.

Hello @tuzu123, I'm not really an expert so I'm sorry if I can't be of great assistance, but the code :

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

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

geomtime = dde.geometry.Interval(0, 1)
bc_u_0 = dde.PeriodicBC(geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=0, component=0)
bc_u_1 = dde.PeriodicBC(geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=1, component=0)

data = dde.data.PDE(geomtime, pde, [bc_u_0,bc_u_1], 100, 2, num_test=200)

worked well for defining the boundary conditions u(0,t)=u(1,t) and u'(0)=u'(1) that's why I improvised my reply. Good luck !

See this for reference: https://github.com/lululxvi/deepxde/issues/457

lululxvi commented 2 years ago

@tuzu123 Your code looks good.

tuzu123 commented 2 years ago

Hi Lu Sorry to bother you again. I tried a lot but I still can't find the right way to solve my problem. At First, the complete PDE of 1D phononic crystal shows as follows: 2022-05-05

  1. When the IC was treated as soft constraints, the loss of IC seemed didn't converge. soft constranits code:

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

a = 0.4 k = 0

muA = 9.4*1e4 muB = 0.468 delta_t = 0.01

position = np.linspace(0, a, 200) ex_input = np.zeros([200]) ex_input[0:100] = 1979.1 ex_input[100:200] = 6

def ex_func(x): spline = sp.interpolate.Rbf( position, ex_input, function="linear", smooth=0, episilon=0 ) return spline(x[0:, 0])

def get_initial_loss(model): model.compile("adam", lr=0.001) losshistory, train_state = model.train(0) return losshistory.loss_train[0]

def pde(x, y, ex): dy_tt = dde.grad.hessian(y, x, i=1, j=1) dy_xx = dde.grad.hessian(y, x, i=0, j=0) return dy_tt - ex * 2 dy_xx

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

def ini_func1(x): x1, t1 = np.split(x, 2, axis=1) ini_x = x1 0 for i in range(5): ini_x = ini_x + 0.001 np.abs(k + i2np.pi/a) np.cos((k + i2np.pi/a)x1) return ini_x

def ini_func2(x): x1, t1 = np.split(x, 2, axis=1) ini_x = x1 0 c = ex_func(x)[:, np.newaxis] for i in range(5): ini_x = ini_x + 0.001 np.abs(k + i2np.pi/a) np.cos((k + i2np.pi/a)x1 - c*t1) return ini_x

def feature_transform(x): return tf.concat( [(x[:, 0:1] - 0.2) 2 np.sqrt(3), (x[:, 1:2] - 0.5) 2 np.sqrt(3)], axis=1 )

geom = dde.geometry.Interval(0, a) timedomain = dde.geometry.TimeDomain(0, 1) geomtime = dde.geometry.GeometryXTime(geom, timedomain)

bc1 = dde.icbc.BlochPeriodicBC(geomtime, 0, boundary_l, derivative_order=0, component=0, phase=k)

bc2 = dde.icbc.BlochPeriodicBC(geomtime, 0, boundary_l, derivative_order=1, component=0, phase=k, muleft=muA, muright=muB)

ic1 = dde.icbc.IC(geomtime, inifunc1, lambda , on_initial: on_initial)

observe_x = np.vstack((np.linspace(0, a, num=360), np.full((360), delta_t))).T ic2 = dde.icbc.PointSetBC(observe_x, ini_func2(observe_x), component=0)

data = dde.data.TimePDE( geomtime, pde, [bc1, bc2, ic1, ic2], num_domain=360, num_boundary=360, num_initial=360, anchors=observe_x, num_test=10000, auxiliary_var_function=ex_func, )

layer_size = [2] + [100] 4 + [1] activation = "tanh" initializer = "Glorot uniform" net = dde.nn.STMsFFN( layer_size, activation, initializer, sigmas_x=[1], sigmas_t=[1, 10] ) net.apply_feature_transform(feature_transform) net.apply_output_transform(lambda x, y: y10)

model = dde.Model(data, net) initial_losses = get_initial_loss(model) loss_weights = 5 / initial_losses model.compile( "adam", lr=0.001, loss_weights=loss_weights, decay=("inverse time", 2000, 0.9), ) pde_residual_resampler = dde.callbacks.PDEResidualResampler(period=1) losshistory, train_state = model.train( epochs=20000, callbacks=[pde_residual_resampler], display_every=500 )

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

And the loss of train steps:

Initializing variables... Training model... Step Train loss Test loss Test metric 0 [2.43e+16, 2.46e+00, 9.59e+11, 1.26e+00, 1.57e+00] [3.31e+16, 2.46e+00, 9.59e+11, 1.26e+00, 1.57e+00] []
Best model at step 0: train loss: 2.43e+16 test loss: 3.31e+16 test metric: [] 'train' took 2.194823 s Compiling model... 'compile' took 1.472080 s Initializing variables... Training model... Step Train loss Test loss Test metric 0 [5.19e+00, 9.82e+00, 3.50e+00, 4.14e+00, 1.24e+00] [6.50e+00, 9.82e+00, 3.50e+00, 4.14e+00, 1.24e+00] []
500 [9.11e-03, 2.56e-03, 2.28e-03, 1.34e-02, 1.11e-02] [2.00e-02, 2.56e-03, 2.28e-03, 1.34e-02, 1.11e-02] []
1000 [4.37e-03, 1.07e-03, 8.96e-04, 1.34e-02, 1.09e-02] [8.72e-03, 1.07e-03, 8.96e-04, 1.34e-02, 1.09e-02] []
1500 [2.89e-03, 6.23e-04, 5.09e-04, 1.35e-02, 1.08e-02] [5.20e-03, 6.23e-04, 5.09e-04, 1.35e-02, 1.08e-02] []
2000 [2.28e-03, 4.10e-04, 3.59e-04, 1.35e-02, 1.08e-02] [3.52e-03, 4.10e-04, 3.59e-04, 1.35e-02, 1.08e-02] []
2500 [1.89e-03, 2.96e-04, 2.81e-04, 1.35e-02, 1.08e-02] [2.56e-03, 2.96e-04, 2.81e-04, 1.35e-02, 1.08e-02] []
3000 [1.63e-03, 2.20e-04, 2.28e-04, 1.35e-02, 1.08e-02] [1.94e-03, 2.20e-04, 2.28e-04, 1.35e-02, 1.08e-02] []
3500 [1.43e-03, 1.77e-04, 1.96e-04, 1.35e-02, 1.08e-02] [1.52e-03, 1.77e-04, 1.96e-04, 1.35e-02, 1.08e-02] []
4000 [1.29e-03, 1.44e-04, 1.72e-04, 1.35e-02, 1.08e-02] [1.23e-03, 1.44e-04, 1.72e-04, 1.35e-02, 1.08e-02] []
4500 [1.21e-03, 1.24e-04, 1.55e-04, 1.35e-02, 1.08e-02] [1.01e-03, 1.24e-04, 1.55e-04, 1.35e-02, 1.08e-02] []
5000 [1.20e-03, 1.08e-04, 1.43e-04, 1.35e-02, 1.07e-02] [8.69e-04, 1.08e-04, 1.43e-04, 1.35e-02, 1.07e-02] []
5500 [1.14e-03, 9.70e-05, 1.34e-04, 1.35e-02, 1.07e-02] [7.40e-04, 9.70e-05, 1.34e-04, 1.35e-02, 1.07e-02] []
6000 [1.07e-03, 9.00e-05, 1.27e-04, 1.35e-02, 1.07e-02] [6.41e-04, 9.00e-05, 1.27e-04, 1.35e-02, 1.07e-02] []
6500 [1.05e-03, 8.52e-05, 1.20e-04, 1.35e-02, 1.07e-02] [5.64e-04, 8.52e-05, 1.20e-04, 1.35e-02, 1.07e-02] []
7000 [1.04e-03, 8.08e-05, 1.17e-04, 1.35e-02, 1.07e-02] [4.94e-04, 8.08e-05, 1.17e-04, 1.35e-02, 1.07e-02] []
7500 [9.67e-04, 9.10e-05, 1.12e-04, 1.35e-02, 1.07e-02] [4.30e-04, 9.10e-05, 1.12e-04, 1.35e-02, 1.07e-02] []
8000 [9.97e-04, 8.78e-05, 1.11e-04, 1.35e-02, 1.07e-02] [3.78e-04, 8.78e-05, 1.11e-04, 1.35e-02, 1.07e-02] []
8500 [9.68e-04, 7.55e-05, 1.09e-04, 1.35e-02, 1.07e-02] [3.31e-04, 7.55e-05, 1.09e-04, 1.35e-02, 1.07e-02] []
9000 [9.62e-04, 7.63e-05, 1.07e-04, 1.35e-02, 1.07e-02] [3.02e-04, 7.63e-05, 1.07e-04, 1.35e-02, 1.07e-02] []
9500 [9.42e-04, 7.04e-05, 1.03e-04, 1.35e-02, 1.07e-02] [2.65e-04, 7.04e-05, 1.03e-04, 1.35e-02, 1.07e-02] []
10000 [8.91e-04, 7.04e-05, 1.01e-04, 1.36e-02, 1.08e-02] [2.40e-04, 7.04e-05, 1.01e-04, 1.36e-02, 1.08e-02] []
10500 [9.00e-04, 3.30e-04, 1.02e-04, 1.36e-02, 1.08e-02] [2.38e-04, 3.30e-04, 1.02e-04, 1.36e-02, 1.08e-02] []
11000 [9.18e-04, 6.75e-05, 1.00e-04, 1.35e-02, 1.07e-02] [1.97e-04, 6.75e-05, 1.00e-04, 1.35e-02, 1.07e-02] []
11500 [8.99e-04, 8.05e-05, 1.00e-04, 1.35e-02, 1.07e-02] [1.83e-04, 8.05e-05, 1.00e-04, 1.35e-02, 1.07e-02] []
12000 [9.55e-04, 9.47e-05, 9.89e-05, 1.35e-02, 1.07e-02] [1.72e-04, 9.47e-05, 9.89e-05, 1.35e-02, 1.07e-02] []
12500 [9.30e-04, 6.17e-05, 9.91e-05, 1.35e-02, 1.07e-02] [1.53e-04, 6.17e-05, 9.91e-05, 1.35e-02, 1.07e-02] []
13000 [8.65e-04, 1.08e-04, 9.60e-05, 1.37e-02, 1.09e-02] [1.39e-04, 1.08e-04, 9.60e-05, 1.37e-02, 1.09e-02] []
13500 [9.00e-04, 6.99e-05, 9.68e-05, 1.35e-02, 1.07e-02] [1.32e-04, 6.99e-05, 9.68e-05, 1.35e-02, 1.07e-02] []
14000 [9.12e-04, 6.33e-05, 9.61e-05, 1.35e-02, 1.07e-02] [1.22e-04, 6.33e-05, 9.61e-05, 1.35e-02, 1.07e-02] []
14500 [1.00e-03, 6.89e-05, 1.05e-04, 1.40e-02, 1.11e-02] [1.21e-04, 6.89e-05, 1.05e-04, 1.40e-02, 1.11e-02] []
15000 [8.88e-04, 7.27e-05, 9.65e-05, 1.35e-02, 1.08e-02] [1.10e-04, 7.27e-05, 9.65e-05, 1.35e-02, 1.08e-02] []
15500 [9.61e-04, 8.41e-05, 9.88e-05, 1.37e-02, 1.09e-02] [1.11e-04, 8.41e-05, 9.88e-05, 1.37e-02, 1.09e-02] []
16000 [1.06e-03, 4.24e-04, 1.15e-04, 1.51e-02, 1.19e-02] [1.59e-04, 4.24e-04, 1.15e-04, 1.51e-02, 1.19e-02] []
16500 [8.65e-04, 7.48e-05, 9.51e-05, 1.35e-02, 1.08e-02] [9.71e-05, 7.48e-05, 9.51e-05, 1.35e-02, 1.08e-02] []
17000 [8.88e-04, 8.65e-05, 9.41e-05, 1.35e-02, 1.07e-02] [9.15e-05, 8.65e-05, 9.41e-05, 1.35e-02, 1.07e-02] []
17500 [9.19e-04, 2.26e-04, 1.00e-04, 1.36e-02, 1.08e-02] [1.07e-04, 2.26e-04, 1.00e-04, 1.36e-02, 1.08e-02] []
18000 [9.00e-04, 6.26e-05, 9.28e-05, 1.35e-02, 1.07e-02] [8.50e-05, 6.26e-05, 9.28e-05, 1.35e-02, 1.07e-02] []
18500 [8.97e-04, 1.23e-04, 9.54e-05, 1.35e-02, 1.07e-02] [8.57e-05, 1.23e-04, 9.54e-05, 1.35e-02, 1.07e-02] []
19000 [8.96e-04, 6.58e-05, 9.28e-05, 1.35e-02, 1.07e-02] [7.80e-05, 6.58e-05, 9.28e-05, 1.35e-02, 1.07e-02] []
19500 [9.12e-04, 6.20e-05, 9.45e-05, 1.35e-02, 1.07e-02] [7.68e-05, 6.20e-05, 9.45e-05, 1.35e-02, 1.07e-02] []
20000 [9.24e-04, 6.04e-05, 9.50e-05, 1.35e-02, 1.07e-02] [7.50e-05, 6.04e-05, 9.50e-05, 1.35e-02, 1.07e-02] []
Best model at step 19500: train loss: 2.53e-02 test loss: 2.44e-02 test metric: []

2. After, I try to construct hard constraints for the IC, but `ex_func(x)` seemed don't work in `modify_output(X, y)`.
Hard constrains code:

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

a = 0.4 k = 0

muA = 9.4*1e4 muB = 0.468 delta_t = 0.01

position = np.linspace(0, a, 200) ex_input = np.zeros([200]) ex_input[0:100] = 1979.1 ex_input[100:200] = 6

def ex_func(x): spline = sp.interpolate.Rbf( position, ex_input, function="linear", smooth=0, episilon=0 ) return spline(x[0:, 0])

def get_initial_loss(model): model.compile("adam", lr=0.001) losshistory, train_state = model.train(0) return losshistory.loss_train[0]

def pde(x, y, ex): dy_tt = dde.grad.hessian(y, x, i=1, j=1) dy_xx = dde.grad.hessian(y, x, i=0, j=0) return dy_tt - ex * 2 dy_xx

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

def modify_output(X, y): x, t = X[:, 0:1], X[:, 1:2] c = ex_func(X)[:, np.newaxis] ini_x = x 0 for i in range(5): ini_x = ini_x + 0.01 tf.abs(k + i2np.pi/a) tf.cos((k + i2np.pi/a)x) ini_x_t = x 0 for i in range(5): ini_x_t = ini_x_t + 0.01 tf.abs(k + i2np.pi/a) tf.cos((k + i2np.pi/a)x - c*delta_t)

return ini_x * (t-delta_t) + ini_x_t * t + y * (t-delta_t) * t * 10

def feature_transform(x): return tf.concat( [(x[:, 0:1] - 0.2) 2 np.sqrt(3), (x[:, 1:2] - 0.5) 2 np.sqrt(3)], axis=1 )

geom = dde.geometry.Interval(0, a) timedomain = dde.geometry.TimeDomain(0, 1) geomtime = dde.geometry.GeometryXTime(geom, timedomain)

bc1 = dde.icbc.BlochPeriodicBC(geomtime, 0, boundary_l, derivative_order=0, component=0, phase=k)

bc2 = dde.icbc.BlochPeriodicBC(geomtime, 0, boundary_l, derivative_order=1, component=0, phase=k, muleft=muA, muright=muB)

data = dde.data.TimePDE( geomtime, pde, [bc1, bc2], num_domain=360, num_boundary=360, num_initial=360, num_test=10000, auxiliary_var_function=ex_func, )

layer_size = [2] + [100] * 4 + [1] activation = "tanh" initializer = "Glorot uniform" net = dde.nn.STMsFFN( layer_size, activation, initializer, sigmas_x=[1], sigmas_t=[1, 10] ) net.apply_feature_transform(feature_transform) net.apply_output_transform(modify_output)

model = dde.Model(data, net) initial_losses = get_initial_loss(model) loss_weights = 5 / initial_losses model.compile( "adam", lr=0.001, loss_weights=loss_weights, decay=("inverse time", 2000, 0.9), ) pde_residual_resampler = dde.callbacks.PDEResidualResampler(period=1) losshistory, train_state = model.train( epochs=20000, callbacks=[pde_residual_resampler], display_every=500 )

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

Report errors:

Using backend: tensorflow.compat.v1 WARNING:tensorflow:From D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\tensorflow\python\compat\v2_compat.py:101: disable_resource_variables (from tensorflow.python.ops.variable_scope) is deprecated and will be removed in a future version. Instructions for updating: non-resource variables are not supported in the long term WARNING:tensorflow:From D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\deepxde\nn\initializers.py:118: The name tf.keras.initializers.he_normal is deprecated. Please use tf.compat.v1.keras.initializers.he_normal instead. D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\skopt\sampler\sobol.py:250: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+362=362. total_n_samples)) Warning: 10000 points required, but 10048 points sampled. Compiling model... Building Spatio-temporal Multi-scale Fourier Feature Network... D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\keras\legacy_tf_layers\core.py:236: UserWarning: tf.layers.dense is deprecated and will be removed in a future version. Please use tf.keras.layers.Dense instead. warnings.warn('tf.layers.dense is deprecated and ' D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\keras\engine\base_layer_v1.py:1676: UserWarning: layer.apply is deprecated and will be removed in a future version. Please use layer.__call__ method instead. warnings.warn('layer.apply is deprecated and ' Traceback (most recent call last): File "D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\IPython\core\interactiveshell.py", line 3524, in run_code exec(code_obj, self.user_global_ns, self.user_ns) File "", line 1, in runfile('D:/BaiduNetdiskWorkspace/DeepXDE/Wave_prop_1D_2materials/Wave_1d_2m.py', wdir='D:/BaiduNetdiskWorkspace/DeepXDE/Wave_prop_1D_2materials') File "D:\software\Pycharm\PyCharm 2018.3.2\helpers\pydev_pydev_bundle\pydev_umd.py", line 198, in runfile pydev_imports.execfile(filename, global_vars, local_vars) # execute the script File "D:\software\Pycharm\PyCharm 2018.3.2\helpers\pydev_pydev_imps_pydev_execfile.py", line 18, in execfile exec(compile(contents+"\n", file, 'exec'), glob, loc) File "D:/BaiduNetdiskWorkspace/DeepXDE/Wave_prop_1D_2materials/Wave_1d_2m.py", line 156, in initial_losses = get_initial_loss(model) File "D:/BaiduNetdiskWorkspace/DeepXDE/Wave_prop_1D_2materials/Wave_1d_2m.py", line 33, in get_initial_loss model.compile("adam", lr=0.001) File "D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\deepxde\utils\internal.py", line 22, in wrapper result = f(*args, **kwargs) File "D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\deepxde\model.py", line 108, in compile self._compile_tensorflow_compat_v1(lr, loss_fn, decay, loss_weights) File "D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\deepxde\model.py", line 124, in _compile_tensorflow_compat_v1 self.net.build() File "D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\deepxde\nn\tensorflow_compat_v1\msffn.py", line 194, in build self.y = self._output_transform(self.x, self.y) File "D:/BaiduNetdiskWorkspace/DeepXDE/Wave_prop_1D_2materials/Wave_1d_2m.py", line 91, in modify_output c = ex_func(X)[:, np.newaxis] File "D:/BaiduNetdiskWorkspace/DeepXDE/Wave_prop_1D_2materials/Wave_1d_2m.py", line 30, in ex_func return spline(x[0:, 0]) File "D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\scipy\interpolate\rbf.py", line 279, in call args = [np.asarray(x) for x in args] File "D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\scipy\interpolate\rbf.py", line 279, in args = [np.asarray(x) for x in args] File "D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\numpy\core_asarray.py", line 83, in asarray return array(a, dtype, copy=False, order=order) File "D:\software\Anaconda\envs\DeepXDE_tf2\lib\site-packages\tensorflow\python\framework\ops.py", line 870, in array " a NumPy call, which is not supported".format(self.name)) NotImplementedError: Cannot convert a symbolic Tensor (strided_slice_7:0) to a numpy array. This error may indicate that you're trying to pass a Tensor to a NumPy call, which is not supported


After changing the version of NumPy, the problem still exists. I think the problem comes from `ex_func(x)` in `modify_output(X, y)`. But I don't know how to solve it.

In addition, BlochPeriodicBC defines:

class BlochPeriodicBC(BC):

def __init__(self, geom, component_x, on_boundary, derivative_order=0, component=0, phase=0, muleft=1, muright=1):
    super().__init__(geom, on_boundary, component)
    self.component_x = component_x
    self.derivative_order = derivative_order
    self.phase = phase
    self.muleft = muleft
    self.muright = muright
    if derivative_order > 1:
        raise NotImplementedError(
            "BlochPeriodicBC only supports derivative_order 0 or 1."
        )

def collocation_points(self, X):
    X1 = self.filter(X)
    X2 = self.geom.periodic_point(X1, self.component_x)
    return np.vstack((X1, X2))

def error(self, X, inputs, outputs, beg, end):
    mid = beg + (end - beg) // 2
    if self.derivative_order == 0:
        yleft = outputs[beg:mid, self.component: self.component + 1]
        yright = outputs[mid:end, self.component: self.component + 1]
    else:
        dydx = grad.jacobian(outputs, inputs, i=self.component, j=self.component_x)
        yleft = dydx[beg:mid]
        yright = dydx[mid:end]
    return yleft * np.cos(self.phase) * self.muleft - yright * self.muright


Thank you in advance.
lululxvi commented 2 years ago

In modify_output, you cannot use numpy. You should use TensorFlow/PyTorch etc.