OperatorBC for varied PDE systems applied on different clouds of anchor points

[sashahexe]

Dear @lululxvi,

Following the issue #780 I also have a question about OperatorBC

I am trying to implement a problem of a flow in a porous medium where the medium is split in let's say three point clouds corresponding to pores, solid and edges in between. Each of the point clouds needs its own set of PDEs. So I want to train PINN on these three clouds using them as anchor points.

I was trying to set the physics conditions for each of the point clouds using following pieces:

def pde_edges(x, u):
    c, q = u[:, 0:1], u[:, 1:2] #have to be >= 0. 
   ...
    return [eq_a, eq_b]


def pde_pores(x, u):
    c, q = u[:, 0:1], u[:, 1:2] #have to be >= 0. 
   ...
    return [eq_a, eq_b]
...

Help function to define if a point x belongs to a certain points cloud. gas_anchors is a (N,3) numpy array of N points in (x1, x2, t) space:

def inside_pores(x, on_boundary):
    return np.isclose(x, gas_anchors)

Put condition to the points in each cloud:

condition_pores = dde.OperatorBC(geomtime, pde_pores, inside_pores)
...

Put all together and test if it can be digested:

data = dde.data.TimePDE(
    geomtime,
    None,
    [bc_a_left, bc_a_right, ..., ic1, ic2,...
        condition_pores],
    num_domain= 0,
    num_boundary=400,
    num_initial=0,
    #anchors = [gas_anchors, edge_anchors, solid_anchors]
)

Unfortunately, I keep getting errors at this last stage. It seems this condition_pores doesn't work well. Would you have a hint how to define this OperatorBC condition with a given pde on a set of points within the geomtime domain?

The error I get:

File ~/anaconda3/lib/python3.8/site-packages/deepxde/icbc/boundary_conditions.py:47, in BC.filter(self, X)
46 def filter(self, X):
---> 47 return X[self.on_boundary(X, self.geom.on_boundary(X))]

IndexError: too many indices for array: array is 2-dimensional, but 3 were indexed

[forxltk]

Use 2 OperatorBCs for eq_a and eq_b respectively？

[sashahexe]

Thank you for the suggestion, @forxltk.

interestingly it worked to some extend (still getting an error but now at a latter stage, see below) in the following modification.

• I've split the pde() into two so each has only one equation as an output:

def pde_pores_c(x, u, X=gas_anchors):
    c, q = u[:, 0:1], u[:, 1:2] 
      
    
    dc_xx = dde.grad.hessian(u, x, component=0, i=0, j=0)
    dc_yy = dde.grad.hessian(u, x, component=0, i=1, j=1)
    dc_t = dde.grad.jacobian(u, x, i=0, j=2)

    eq_a = ...

    return eq_a

def pde_pores_q(x, u, X=gas_anchors):
    c, q = u[:, 0:1], u[:, 1:2] 
        
    eq_b = q

    return eq_b

• I tried np.allclose() instead of np.isclose() as it :

def inside_pores(x, _):
    return np.allclose(x, gas_anchors)

split the conditions:

condition_pores_c = dde.OperatorBC(geomtime, pde_pores_c, inside_pores)
condition_pores_q = dde.OperatorBC(geomtime, pde_pores_q, inside_pores)

Now the data module executes without an error:

data = dde.data.TimePDE(
    geomtime,
    None,
    [bc_a_left, bc_a_right, bc_b_left, bc_b_right, bc_a_top, 
        bc_a_bottom, bc_b_top, bc_b_bottom, ic1, ic2, condition_pores_c, condition_pores_q],
    num_domain= 0,
    num_boundary=400,
    num_initial=0,
    anchors = gas_anchors
)

BUT after the first step the error occurs:

layer_size = [3] + [32] * 3 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
# net.apply_output_transform(modify_output)  
    

model = dde.Model(data, net)
model.compile('adam', lr=0.001)
model.train(epochs=10000)   

losshistory, train_state = model.train(epochs=1000) 
dde.saveplot(losshistory, train_state, issave=True, isplot=True)

Initializing variables...
Training model...

Step Train loss Test loss Test metric
0 [1.05e+00, 2.95e-01, 9.32e-01, 4.81e-01, 4.61e-01, 2.95e-01, 1.84e-01, 4.81e-01, 1.56e-02, 7.05e-02, nan, nan] [1.05e+00, 2.95e-01, 9.32e-01, 4.81e-01, 4.61e-01, 2.95e-01, 1.84e-01, 4.81e-01, 1.56e-02, 7.05e-02, nan, nan] []

Best model at step 0:
train loss: inf
test loss: inf
test metric:

'train' took 0.401152 s

Training model...

Step Train loss Test loss Test metric
1 [9.20e-01, 3.02e-01, 8.51e-01, 4.34e-01, 4.35e-01, 3.02e-01, 2.12e-01, 4.34e-01, 2.90e-02, 6.99e-02, nan, nan] [9.20e-01, 3.02e-01, 8.51e-01, 4.34e-01, 4.35e-01, 3.02e-01, 2.12e-01, 4.34e-01, 2.90e-02, 6.99e-02, nan, nan] []

Best model at step 0:
train loss: inf
test loss: inf
test metric:
...

Saving loss history to /home/.../loss.dat ...
Saving training data to /home/.../train.dat ...
Saving test data to /home/.../test.dat ...

File <array_function internals>:5, in concatenate(*args, **kwargs)

ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 1 dimension(s)

• it seems the output for the file test.dat doesn't fit the expectations after the step 1 in terms of dimensions.
  Would you have any ideas?

• second thing is: note that the NN doesn't calculate the residual for the OperatorBC conditions, why do you think it happens?

[sashahexe]

As an alternative I thought to come from the other end. Instead of using special conditions on the certain points via OperatorBC I could use a condition on x in the pde.

Now I wrote a function to check if a tensorflow tensor x is in a certain points cloud (numpy array for my anchor points):

tf.compat.v1.enable_eager_execution()

def is_in(anchors, x):
    return np.allclose(x.numpy(), anchors)

And implemented it into the pde()

def pde_pores(x, u):  
    c, q = u[:, 0:1], u[:, 1:2] 

    eq_c_pores = ...
    eq_q_pores = ...
    
    return [eq_c_pores, eq_q_pores]


def pde_edges(x, u):
    c, q = u[:, 0:1], u[:, 1:2] 

    eq_a = ...
    eq_b = ...
    
    return [eq_a, eq_b]

def pde_solid(x, u):
    c, q = u[:, 0:1], u[:, 1:2] 
          
    eq_b = ...
    eq_a = ...

    return [eq_a, eq_b]

def pde(x, u):
    if is_in(gas_anchors,x):
        return pde_pores(x,u)
    
    elif is_in(edge_anchors,x):
        return pde_edges(x,u)
    
    else:
        return pde_solid(x,u)
    

But tensorflow seems to be unhappy when trying to use pde():

---

RuntimeError Traceback (most recent call last)
/home/sasha/Documents/Nextcloud/RESEARCH/openPNM/Geometry-and-diffusion-in-paper/Scripts/diffusion_trials/pinns/2d_hetero_decomposed/drafts2.ipynb Cell 11' in <cell line: 4>()
2 activation = "tanh"
3 initializer = "Glorot uniform"
----> 4 net = dde.maps.FNN(layer_size, activation, initializer)
5 # net.apply_output_transform(modify_output)
8 model = dde.Model(data, net)

File ~/anaconda3/lib/python3.8/site-packages/deepxde/nn/tensorflow_compat_v1/fnn.py:27, in FNN.init(self, layer_sizes, activation, kernel_initializer, regularization, dropout_rate, batch_normalization, layer_normalization, kernel_constraint, use_bias)
15 def init(
16 self,
17 layer_sizes,
(...)
25 use_bias=True,
26 ):
---> 27 super().init()
28 self.layer_size = layer_sizes
29 self.activation = activations.get(activation)

File ~/anaconda3/lib/python3.8/site-packages/deepxde/nn/tensorflow_compat_v1/nn.py:12, in NN.init(self)
11 def init(self):
---> 12 self.training = tf.placeholder(tf.bool)
13 self.regularizer = None
15 self._auxiliary_vars = tf.placeholder(config.real(tf))
...
-> 3298 raise RuntimeError("tf.placeholder() is not compatible with "
3299 "eager execution.")
3301 return gen_array_ops.placeholder(dtype=dtype, shape=shape, name=name)

RuntimeError: tf.placeholder() is not compatible with eager execution.

Would you know how to fix it? Or if it's better to get back to separate OperatorBC conditions for each anchor points cloud?

[forxltk]

1. I think there are no sample points for OperatorBC. You should check inside_pores function.
2. As the error message, eager execution should be disable.

[sashahexe]

@forxltk

1. thanks I will check it.
2. would you know how to check if a certain tensorflow tensor is in a numpy array without using eager?

[forxltk]

No, but what about choose tf2 as backend so you can still use eager?

[sashahexe]

@forxltk thanks for your valuable suggestions, not everything is working but I have a feeling being closer to the end of the tunnel.

1. Cannot get it working with the modified pde(). I have activated tf2 and updated my Tensorflow, also made sure to activate eager by tf.config.run_functions_eagerly(True). Activating of eager seems to be a problem for deepxde. As soon as I activate it I get an error even in a simplest case without any conditions inside the pde() function. When eager is not explicitly activated by tf.config.run_functions_eagerly(True) I do not get this error but then I cannot check if an x is in a numpy array. Looks like a endless loop...

def pde(x, u):
    c, q = u[:, 0:1], u[:, 1:2] 

    eq_a = ...
    eq_b = ...
    
    return [eq_a, eq_b]

--> 561 self._test()
562 self.callbacks.on_train_begin()
563 if optimizers.is_external_optimizer(self.opt_name):

File ~/anaconda3/lib/python3.8/site-packages/deepxde/model.py:693, in Model._test(self)
...
--> 244 key = (y.ref(), xs.ref(), component)
245 elif backend_name in ["pytorch", "paddle"]:
246 key = (y, xs, component)

AttributeError: 'numpy.ndarray' object has no attribute 'ref'

Version of pde(x,u) with a condition and deactivated eager gives (of course) following error:

def is_in(anchors, x):
    return x.numpy() in anchors

def pde_cond(x, u):

    if is_in(gas_anchors,x):
        return pde_pores(x,u)
    
    elif is_in(edge_anchors,x):
        return pde_edges(x,u)
    
    else:
        return pde_solid(x,u)

AttributeError: 'Tensor' object has no attribute 'numpy'

I also tried auxiliary_var_function for the condition on x but comparing to a numpy array runs in the same problem with a tensor vs numpy. Not sure if there is a clever way to make a condition without switching to numpy arrays somehow?

2. Version with the OperatorBC. I have checked the inside_pores() function and modified it:

def inside_pores(x, _):
    return x in gas_anchors

Now the model compiles and trains. But in the end I see that only the anchor points at the time = 0 and time =1 were actually taken for the training (and testing). So here a refinement is still needed. I guess it comes from the nature of dde.OperatorBC(geomtime, pde_pores_q, on_boundary = inside_pores) where the last argument is actually an on_boundary function. I thought I could overwrite it somehow with my own inside_pores but it works weirdly. Would you have an idea?

[forxltk]

1. It seems that tf.config.run_functions_eagerly(True) convert the input x to ndarray, not Tensor or EagerTensor. Maybe you can modified the source code of ndarray or .ref?
2. No. It should work for all sample points if you only return something like np.isclose(x[0], 0) without on_boundary.

[sashahexe]

@forxltk

Here is my set up for sampling from anchor points.

gas_anchors is a (N,dim) np array. Analogously, the other two edge_anchors and edge_anchors have their own length:

gas_anchors.shape
>> (17140, 3)

1. help function to check from which cloud the given anchor point is:

def find_nearest_array(x):
    if x in gas_anchors:
        return 'gas'
    elif x in edge_anchors:
        return 'edge'
    elif x in solid_anchors:
        return 'solid'

2. function to return True if an anchor point belongs to pores:

def inside_pores(x, _):
    answer = find_nearest_array(x)

    if answer == 'gas':
        return True
    
    else:
        return False

3. Condition to use OperatorBC on the points which belong to pores:

condition_pores_c = dde.OperatorBC(geomtime, pde_pores_c, on_boundary = inside_pores)

4. geomtime domain:

L = 1.
time = 1.
geom = dde.geometry.Rectangle([0, 0], [L, L])
timedomain = dde.geometry.TimeDomain(0, time)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

6. put all together:

data = dde.data.TimePDE(
    geomtime,
    None,
    [bc_a_left, bc_a_right, ..., condition_pores_c, ...],
    num_domain= 0,
    num_boundary=400,
    num_initial=0,
    anchors = anchors_
)

After training I see that the model took only the anchors which have time = 0 and time = 1:

Although the gas_anchors look as following:

I don't use np.isclose() because I wanted to make it easier for the sampler, it only needs to take exactly the points which belong to the corresponding point clouds. Also np.isclose(x, gas_anchors) compares 1d values, and gives something like [True, False, True] for each [x1,x2,t] point so it needs further evaluation whether all are True or not. I could also use np.allclose(x, gas_anchors) to get True or False, but since I am only going to train for the point from a predefined cloud, I thought it would be better to return x in gas_pores.

Would you see a flaw in my setup?

[forxltk]

So you should still check inside_pores. For example, does the model take all the sample points if inside_pores always return true. Then you may find where the problem is.

[sashahexe]

Dear @forxltk thank you for the hints! Indeed I found the bottle neck with the inside_pores function and could finally get the code to the working stage with:

dim = anchors_.shape[1] #equals 3 for 2D pictures and 4 for 3D pictures

def inside_pores(x, _):
    check = np.sum(np.isclose(x, gas_anchors), axis = 1)
    return check[check == dim].size > 0

[lululxvi]

@sashahexe You cannot use not eager mode in DeepXDE. If you want to print tensor values, you may try PyTorch backend.

[sashahexe]

Thank you , will try it.

[sashahexe]

Dear @lululxvi, dear @forxltk, I have a question concerning pytorch implementation, maybe you will know it please:

I am trying to set different PDE systems for different collocation points. *I could already implement it using OperatorBC, but I hoped to do it over PDEs to spare on computational cost as my domains are pretty large. Also OperatorBC needs the PDE system to be split into individual equations, which makes it even more 'clumsy' when expanding to the gPINNs. So I hoped it would be possible to go from another end.

Here is my pde set up:

def classify_x(x):
    return knn.predict(x[:,0:2])

I am classifying my collocation points with a pretrained KNN classifier and check that the tensors prom PyTorch can be passed to my classify_x(x). This part works all right now (thanks to the pytorch background):

label = classify_x(x.detach().numpy())
print(label[label == 0.].shape)
print(label[label == 1.].shape)
print(label[label == 2.].shape)
>>>(926,)
>>>(1338,)
>>>(626,)

This is how I thought to set up the choice of a certain pde system depending on the label of x:

def pde(x, u):

    label = classify_x(x.detach().numpy())

    x_pores = x[label == 0.]
    u_pores = u[label == 0.]


    x_solid = x[label == 1.]
    u_solid = [label == 1.]

    x_edges = [label == 2.]
    u_edges = [label == 2.]

    
    return [pde_pores(x_pores, u_pores), pde_solid(x_solid, u_solid), pde_edges(x_edges, u_edges)]

The specific PDEs systems:

def pde_pores(x, u):  
    c, q = u[:, 0:1], u[:, 1:2] 

    eq_c_pores = ...
    eq_q_pores = ...
    
    return eq_c_pores, eq_q_pores


def pde_edges(x, u):
    c, q = u[:, 0:1], u[:, 1:2] 

    eq_a = ...
    eq_b = ...
    
    return eq_a, eq_b

def pde_solid(x, u):
    c, q = u[:, 0:1], u[:, 1:2] 
          
    eq_b = ...
    eq_a = ...

    return eq_a, eq_b

data:

data = dde.data.TimePDE(
    geomtime,
    pde,
    [],
    num_domain=1000,
    anchors = anchors_,
)

PINN:

layer_size = [3] + [32] * 3 + [2]
activation = 'tanh'
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
loss_weights = [1,1,1,1,1,1,1,1,1,1,1,1] 
model.compile('adam', lr=0.001,loss_weights=loss_weights)
model.train(epochs=2000)  

It compiles but seems to be unhappy about my way of picking the pdes depending on x label:

/content/drive/Othercomputers/My Laptop/DeepXDE/hetero_models/hetero_pytorch_colab/simulations/sim_6_PDEcond_0
Compiling model...
'compile' took 0.000180 s

Warning: epochs is deprecated and will be removed in a future version. Use iterations instead.
Training model...

[1. 1. 2. ... 0. 1. 1.]
(2890,)

---

RuntimeError Traceback (most recent call last)

in ()
12 loss_weights = [1,1,1,1,1,1,1,1,1,1,1,1]
13 model.compile('adam', lr=0.001,loss_weights=loss_weights)
---> 14 model.train(epochs=2000)
15
16 # loss_weights = [1, 100, 1, 100, 100, 100, 100, 1, 1, 1]

16 frames

/usr/local/lib/python3.7/dist-packages/torch/autograd/init.py in grad(outputs, inputs, grad_outputs, retain_graph, create_graph, only_inputs, allow_unused, is_grads_batched)
276 return Variable.execution_engine.run_backward( # Calls into the C++ engine to run the backward pass
277 t_outputs, grad_outputs, retain_graph, create_graph, t_inputs,
--> 278 allow_unused, accumulate_grad=False) # Calls into the C++ engine to run the backward pass
279
280

RuntimeError: One of the differentiated Tensors appears to not have been used in the graph. Set allow_unused=True if this is the desired behavior.

I tried to set up what is suggested namely allow_unused = True in the /usr/local/lib/python3.7/dist-packages/torch/autograd/__init__.py but the error persists:

--> 278 allow_unused=True, accumulate_grad=False) # Calls into the C++ engine to run the backward pass
279
280

RuntimeError: One of the differentiated Tensors appears to not have been used in the graph. Set allow_unused=True if this is the desired behavior.

Would you have an idea how to get around this problem?

[lululxvi]

I haven't tried this kind of thing. I guess the issue is x.detach(). If you detach x from the graph, then PyTorch cannot compute the gradients anymore.