Difficulty in Setting Initial Conditions

[davidsohutskay]

Hello,

I have been trying to write a code to solve a reaction-diffusion type problem. I wanted to compare deepXDE and SciANN but was having difficulty getting my initial conditions set up in SciANN.

In deepXDE, I can use

rho0*tf.cast(tf.math.greater(((x-x_center_scaled)**2)/(x_axis_scaled**2),1),tf.float64)

To set a circle (line segment in 1D for now) in the center of the domain to zero while keeping everything else at a constant value. The result is shown below (sorry for the unlabeled axes - that is x from 0 to 75 and t from 0 to 10).

In SciANN, this doesn't seem to work since the Functionals cannot be acted on by TF operations. In accordance with the Burger eq example, I tried using

0.25*(1 - sign(tt - TOL)) * ((1 - sign(x - (left_boundary+TOL))) + (1 + sign(x - (right_boundary-TOL)))) * ( rho0)

as well as

0.5*(1 - sign(t - TOL)) * rho0 * sign(((x-x_center)**2)/(x_axis**2))

In numpy, both of these give the correct shape

I also tried using

0.5*(1 - sign(t - TOL)) * rho0 * tanh(((x-x_center)**2)/(x_axis**2))

Which should give a smoother boundary. In all cases the whole domain appears as solid rho0, with no area of 0 in the center, as shown below.

I also note that the initial condition loss seems to immediately get stuck at a particular value.

So I am wondering: is there a problem with the way I have set up my boundary conditions, a way to use TF operations or a "greater_than" function in SciANN, or is there another way to set up my boundary conditions?

Thanks,
David

[davidsohutskay]

I will note that the code does work with a "smoother" initial condition like a linear function

0.5*(1 - sign(t - TOL)) * (rho - x/7.5)

which (even with the weird boundary conditions) gives the expected

[sciann]

Hi - thanks for your notes. It is certainly not what I expect. Is it possible to share your code? You can email me privately and I will check that out.
It could be the combination of initializer, optimizer, and batch size that you get faster convergence in one package and not the other. I

[davidsohutskay]

Sure, here is the file. Thank you for your reply.

import numpy as np
import matplotlib.pyplot as plt 
import sciann as sn
from sciann.utils.math import diff, sign, step, tanh, sigmoid
import tensorflow as tf

# Constants
rho0 = 10
D_rho = 0.833
#K_rhorho = 10000
#p_rho = 0.034
#d_rho = p_rho*(1-rho0/K_rhorho)

x_axis = 7.5
y_axis = 7.5
x_center = 37.5
y_center = 37.5

# Variables
x = sn.Variable('x')
t = sn.Variable('t')
rho = sn.Functional('rho', [x,t], 6*[32], 'tanh')

# Loss function (PDE)
L1 = diff(rho, t) - D_rho*diff(rho, x, order=2)

# Initial conditions
TOL = 0.001

# First check t=0
# Second check if x > 30
# Third check if x < 45

# Quick check in Numpy
tt = 0     
xx = np.linspace(0, 75, 100)
#yy =  0.5*0.5*(1 - np.sign(tt - TOL)) * ((1 - np.sign(xx - (30+TOL))) + (1 + np.sign(xx - (45-TOL)))) * ( rho0)
yy =  0.5*(1 - np.sign(tt - TOL)) * rho0 * np.tanh(((xx-x_center)**2)/(x_axis**2))
#yy =  0.5*0.5*(1 - np.sign(tt - TOL)) * ((xx<30) + (xx>45)) * ( rho0)
plt.plot(xx,yy)

# This doesn't work
C1 = 0.5*(1 - sign(t - TOL)) * rho0 * tanh(((x-x_center)**2)/(x_axis**2))
#C1 = 0.5*(1 - sign(t - TOL)) * rho0 * sign(((x-x_center)**2)/(x_axis**2)) 
#C1 = 0.5*0.5*(1 - sign(t - TOL)) * ((1 - sign(x - (30+TOL))) + (1 + sign(x - (45-TOL)))) * ( rho0)

# This does work but is not what I want
#C1 = 0.5*(1 - sign(t - TOL)) * (rho - x/7.5)

# Boundary conditions
C2 = 0.5*(1 - sign(x - (0+TOL))) * (rho - rho0)
C3 = 0.5*(1 + sign(x - (75-TOL))) * (rho - rho0)

# Model
m = sn.SciModel(inputs=[x, t], 
                targets=[L1, C1, C2, C3],
                loss_func="mse", 
                optimizer="adam")

# Training data
x_data, t_data = np.meshgrid(
    np.linspace(0, 75, 100), 
    np.linspace(0, 1, 100)
)

# Set targets (RHS)
data_L1 = 'zeros'
data_C1 =  'zeros'
data_C2 = 'zeros'
data_C3 = 'zeros'
target_data = [data_L1, data_C1, data_C2, data_C3]

# Train the model
h = m.train([x_data, t_data], 
            target_data, 
            learning_rate = 0.001, 
            epochs=1000,
            batch_size = 2**6, 
            save_weights_to='weights.dat', 
            verbose=True)

# Test data
x_test, t_test = np.meshgrid(
    np.linspace(0, 75, 200), 
    np.linspace(0, 1, 100)
)
rho_pred = rho.eval(m, [x_test, t_test])

# Plot the results
fig = plt.figure(figsize=(6, 8))
plt.pcolor(x_test, t_test, rho_pred, cmap='seismic', vmin=0, vmax=10)
plt.xlabel('x')
plt.ylabel('t')
plt.colorbar()

# Plot a surface
x_plot = x_test.reshape(-1,1).flatten() 
t_plot = t_test.reshape(-1,1).flatten() 
rho_plot = rho_pred.reshape(-1,1).flatten() 

my_cmap = plt.get_cmap('hot') 
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_zlim(0,10)
trisurf = ax.plot_trisurf(x_plot, t_plot, rho_plot, 
                     cmap = my_cmap,
                     vmin = 0,
                     vmax = rho0,
                     linewidth = 0.2,  
                     antialiased = True, 
                     edgecolor = 'grey') 
fig.colorbar(trisurf, ax = ax, shrink = 0.5, aspect = 5) 

I tried modifying a few of the conditions - I am just surprised that the loss stays completely constant over the epochs and there is no evidence it picked up the initial condition at all (just flat everywhere). Especially when the linear one with the discontinuous boundary does work.

[sciann]

I spot an issue in the C1:

C1 = 0.5*(1 - sign(t - TOL)) * rho0 * tanh(((x-x_center)**2)/(x_axis**2))

This relation does not impose any condition on the network (rho). Therefore it does not change at all as there is no parameter gets affected by this condition.

Is it possible to have a typo in C1, i.e. rho0 -> rho?

[davidsohutskay]

Ah, you are correct! Good catch, I missed that! The condition should have been

C1 = 0.5*(1 - sign(t - TOL)) * (rho - rho0 * tanh(((x-x_center)**2)/(x_axis**2)))

It works now! Thanks a lot for your help!