Second Order PDEs

[dangalea]

Hi,

I have been looking into DL with physics and came across this project. However, I am having some trouble using this for a dummy case in mind, the wave equation, which is a second-order PDE. In both examples given, first-order PDEs are used. Is this package meant to be used for second-order PDEs, and if so, would you be able to share some pointers?

Regards,
Daniel

[holl-]

Hi,

Any second-order PDE can be transformed into a first-order PDE. In your case, you just need to add the momentum field to your state.

The wave equation would then look something like this (using Euler integration which is not stable!):

def simulate(x, p, dt=.1):
    return x + p * dt, p + field.laplace(x) * dt

where we could initialize the values like so:

height = CenteredGrid(SoftGeometryMask(Sphere(x=50, y=50, radius=5)), extrapolation.PERIODIC, x=100, y=100)
momentum = height * 0

Cheers,
Philipp

[dangalea]

Hi,

thanks for your help. I am trying to learn how to use Phiflow and haven't done calculus for a while, so sorry if some of questions might seem trivial.

I am trying to set up a simulation for a 1D wave equation using Neumann boundary conditions, with a Gaussian bell as initial starting conditions. My code for a normal finite element simulation is as follows:

import numpy as np
import matplotlib.pyplot as plt
import copy

#Initial Condition
sigma = 0.5
c = 1

#Domain Specification
x_length = 5.
dx = 0.1
nx = int(x_length / dx)

#Time Specification
t_length = 250
dt = 0.1
nt = int(t_length / dt)

#Arrays for calculating solution
u_n_plus_1 = np.zeros(nx)
u_n = np.zeros(nx)
u_n_minus_1 = np.zeros(nx)

#Initial conditions
for x in range(len(u_n)):
    u_n[x] = np.exp(-0.5*((x*dx - x_length/2)/sigma)**2)
u_n[0] = 0
u_n[-1] = 0

u_n_minus_1 = copy.deepcopy(u_n)

u_n_hist = []
u_n_hist.append(u_n)

## Time Stepping
for t in range(1, nt):
    for x in range(1, nx-1):
        u_n_plus_1[x] = 2*u_n[x] - u_n_minus_1[x] + (c * dt / dx)**2 * (u_n[x+1] - 2*u_n[x] + u_n[x-1])

    u_n_minus_1 = copy.deepcopy(u_n)
    u_n = copy.deepcopy(u_n_plus_1)

    u_n_hist.append(u_n)

fig, axs = plt.subplots(2, 2)
axs[0, 0].plot(u_n_hist[0])
axs[0, 1].plot(u_n_hist[int(nt/4)])
axs[1, 0].plot(u_n_hist[int(3*nt/4)])
axs[1, 1].plot(u_n_hist[-1])

axs[0, 0].set_title("Initial Condition")
axs[0, 1].set_title("T = " + str(int(nt/4)*dx))
axs[1, 0].set_title("T = " + str(3*int(nt/4)*dx))
axs[1, 1].set_title("T = " + str(t_length))

plt.tight_layout()
plt.show()

Now, I'm trying to replicate this in the PBDL book, i.e. build a NN that takes in x and t as inputs and returns u. As a first step, I thought of getting a Phiflow only simulation running, then adding a NN later. To that end, I'm trying to get a Phiflow only simulation going with the following code:

from phi.flow import *
import numpy as np
import matplotlib.pyplot as plt
import copy

def simulate(x, p, dt=.1):
    return x + p * dt, p + field.laplace(x) * dt

#Initial Condition
sigma = 0.5
c = 1

#Domain Specification
x_length = 5.
dx = 0.1
nx = int(x_length / dx)

#Time Specification
t_length = 200
dt = 0.1
nt = int(t_length / dt)

#Initial conditions
initial_cond = np.zeros(nx)
for x in range(nx):
    initial_cond[x] = np.exp(-0.5*((x*dx - x_length/2)/sigma)**2) 
initial_cond[0] = 0
initial_cond[-1] = 0

initial_cond_phi = math.tensor(initial_cond, spatial('x') ) # convert to phiflow tensor

height = CenteredGrid(SoftGeometryMask(Sphere(x=int(nx/2), y=int(nx/2), radius=5)), extrapolation.PERIODIC, x=nx, y=nx)
momentum = height * 0

x_hist = [initial_cond_phi]
p_hist = [momentum]
for i in range(5):
    x, p = simulate(x_hist[-1], p_hist[-1], dt=.1)
    x_hist.append(x)
    p_hist.append(p)
        
wave_hist = [wave.values.numpy('x') for wave in wave_hist]

fig, axs = plt.subplots(2, 2)
axs[0, 0].plot(wave_hist[0])
axs[0, 1].plot(wave_hist[1])
axs[1, 0].plot(wave_hist[2])
axs[1, 1].plot(wave_hist[-1])

axs[0, 0].set_title("Initial Condition")
axs[0, 1].set_title("T = " + str(int(nt/4)*dx))
axs[1, 0].set_title("T = " + str(3*int(nt/4)*dx))
axs[1, 1].set_title("T = " + str(t_length))

plt.tight_layout()
plt.show()

However, I am getting the following error: TypeError: 'NoneType' object cannot be interpreted as an integer. I'm guessing I'm mixing something up when it comes to the initial conditions. Would you be able to help put me on the right path? Also, would you be able to explain how you transformed the second-order PDE to a first-order PDE?

Thanks.

[holl-]

Sorry for the late reponse. You are storing a Tensor in initial_cond_phi but field.laplace expects a grid.
Try this:

initial_cond_phi = CenteredGrid(math.tensor(initial_cond, spatial('x')), extrapolation.PERIODIC)

Also your initial_cond_phi has only one spatial dimension, x but momentum currently has x and y.