Inverse problem: Heat convection-diffusion equation

[saluisto]

[Hello] community,

i have written a PINN using Deepxde to solve the 1D heat convection-diffusion equation (https://en.wikipedia.org/wiki/Convection–diffusion_equation) where the convection velocity is inversely solved. To estimate the convection velocities over time I use observations for different depths.As boundary conditions i am using observed temperatures that do change in time. I post this code here to help other members of the community that want to do a similar problems with Deepxde and maybe to get some constructive feedback from the community in how i can improve the code. Especially suggestions to increase the efficiency are highly appreciated.

Loss (any ideas to increase the convergence?):

Observed vs. simulated for different depths

Inversely solved convective velocity q in m/day as a function of time (time is in days not minutes):

DATA:temp_profiles.csv

Code

import deepxde as dde
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
```

```
#%% This section reads the measured temperature in different depths to setup the boundary conditions and
# to solve for the advective flow velocity

def read_data(filename):
    # Read the CSV file
    df = pd.read_csv(filename)
    
    obs_data = df.dropna()
    obs_data=obs_data.drop(index=range(4))
    
    return obs_data

df=read_data(filename='temp_profiles.csv')
df['time'] = pd.to_datetime(df['time'], format="%d.%m.%Y %H:%M:%S")
#df.set_index('time', inplace=True)
df = df.drop('depth1', axis=1)
first_time = df['time'].min()
df['time'] = df['time'] - first_time

#Here time is converted from seconds into days. I noticed that with seconds the efficiency was not very good
df['time'] = df['time'].dt.total_seconds()/(60*60*24)  #time is converted in 

time = df['time'].values.reshape(-1,1)

df.set_index('time', inplace=True)

df.plot()
#%%
# Here the heat advection-conduction PDE is beeing defined 

def pde(x, y):
    ke=3.8*60*60*24 #ke is converted from J/(s m K) to J/(day m K)
    roh_sed=2650
    c_sed=1078
    roh_water=1024
    c_water=3986
    D=ke/(roh_sed*c_sed)
    a=roh_water*c_water/(roh_sed*c_sed)
    
    T, q = y[:, 0:1], y[:, 1:2] #Output of the NN is the temperature field T and q fluxes in time
    dy_t = dde.grad.jacobian(y, x, i=0, j=1)
    
    #Not sure why I need the component=0 statement here, but wihtout it the code does not run
    dy_xx = dde.grad.hessian(y, x, component=0, i=0, j=0) 
    dy_dx = dde.grad.jacobian(y, x, i=0, j=0)
    return (
        dy_t
        - D * dy_xx
        + a * q * dy_dx
    )



#%%
#Here I prepare the observations for the BC and the inverse problem
observe_x = np.array([0,0.21,0.42,0.62,0.82,1.01,1.22]).reshape(-1,1)
observe_T = df[['depth2','depth3','depth4','depth5','depth6','depth7','depth8']].values
xx, tt = np.meshgrid(observe_x,time)
X = np.vstack((np.ravel(xx), np.ravel(tt))).T #x,t
observe_T=np.ravel(observe_T).reshape(-1,1)
observe_u = dde.icbc.PointSetBC(X,observe_T,component=0)
#%% Define gometry
L=1.22
max_time=time[len(time)-1]
geom = dde.geometry.Interval(0, L)
#geom = geom.add_points(observe_x)
timedomain = dde.geometry.TimeDomain(0, max_time)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

#Initially I used more points distributed in the domain but i figured out that i can work with the observations only
#as anchors without the same accuracy. Altough i am a bit unsure here what the better strategy is. Maybe you have some
#suggestions.

##%
# data = dde.data.TimePDE(
#     geomtime,
#     pde,
#     [observe_u],
#     num_domain=1000,
#     num_boundary=1000,
#     num_initial=100,
#     anchors=X,
# )

data = dde.data.TimePDE(
    geomtime,
    pde,
    [observe_u],
    anchors=X,
)


net = dde.nn.PFNN([2, [40, 30], [40, 30], [40, 30], 2], "tanh", "Glorot uniform")

model = dde.Model(data, net)
#played around with the learning rate a bit. This value seems to work good but the NN seems to have
#a slow convergence. If you have some ideas to improve the learning procedure I would be grateful 
model.compile("adam", lr=0.001)


losshistory, train_state = model.train(iterations=50000)
dde.saveplot(losshistory, train_state, issave=True, isplot=True)


Y = model.predict(X)  # Y is 3 x 6 matrix, i.e., [[y1(x1), ..., y6(x1)], [y1(x2), ..., y6(x1)], [y1(x3), ..., y6(x3)]]
#%%
#Plotting the different predictions vs observations

for i in range(len(observe_x)): 
    dp = ['depth2','depth3','depth4','depth5','depth6','depth7','depth8']
    columns = ['X[m]', 'Time[min]', 'T[C]', 'q[m/day]']
    sim=np.concatenate((X,Y),axis=1)
    result = pd.DataFrame(sim, columns=columns)
    
    rslt_d1 = result[result['X[m]'] == float(observe_x[i])]
    rslt_d1= rslt_d1.drop('X[m]',axis=1)
    
    rslt_d1.plot(x='Time[min]', y='T[C]', kind='scatter', title='Temp_simulated')
    df[dp[i]].plot(kind='line', title=dp[i])


rslt_d1.plot(x='Time[min]', y='q[m/day]', kind='line', title='Temp_simulated')

```