How can I specify a point-like source for the wave equation?

[ShaikhaTheGreen]

Again thanks for your contribution!

I want to add a source to my wave function simulation at a specific location in the domain. I'm currently writing it in this way:

#Compose the PDE
#eq = pde.WavePDE(speed = 0.01, bc=[{"value": 0}, {"derivative": 0}])
speed = 0.01
frequency = 10
source = f"sin(2*pi*{frequency}*t)"
eq = PDE(
    {
        "u": f"v",
        "v": f"{speed}**2 * laplace(u) + {source}"
    }, bc=[{"value": 0}, {"derivative":0}]
)

where my source f(x) = sin(2*pi*f*t)
I came to know that if I wanted the source f(x) to be at location x=0, then I should implement something like:

using the impulse function.
How can I implement this? Is there any specific syntax to do it?

[david-zwicker]

There is no special implementation of a delta-function, but you can in principle use any spatial field in the definition of the pde. I thus suggest to create a field that mimics the delta-impulse:

grid = pde.CartesianGrid([[0, 16], [0, 16]], 64)
point = pde.ScalarField(grid)  # create empty field
point.add_interpolated([8, 8], 1)  # add a point source at a given location
point.plot()

You can then use this field in your definition of the PDE:

speed = 0.01
frequency = 10
source = f"sin(2*pi*{frequency}*t)"
eq = PDE(
    {
        "u": f"v",
        "v": f"{speed}**2 * laplace(u) + {source} * point"
    }, bc=[{"value": 0}, {"derivative":0}],
    consts={'point': point}
)

Note that you need to specify explicitly what constants are defined for PDE.

[ShaikhaTheGreen]

Is it possible to also test placing the source at the boundary?
I tried:

speed = 0.01
frequency = 10
source = f"sin(pi*x) * exp(-2*pi*{frequency}*t)"
eq = PDE(
    {
        "u": f"v",
        "v": f"{speed}**2 * laplace(u) " 
    }, bc=[{"value": f"sin(pi*x) * exp(-2*pi*{frequency}*t)"}, {"derivative":0}]
)

but I get this error:

 f"Arguments {args} were not defined in expression signature {signature}"
RuntimeError: Arguments {'t'} were not defined in expression signature ['x']

[david-zwicker]

Unfortunately, time-dependent boundary conditions are not (yet) supported. You may open an issue asking for this as a feature-request, but I cannot promise that I will have time to look at this anytime soon. Alternatively, you could of course try to implement it yourself, although the boundary conditions are handled quite deeply in the code.

[ghost]

@engsbk @david-zwicker for time-dependent boundary conditions I was able to botch together an (approximate) solution. The problem I was trying to solve was heat diffusion through a convective and radiatively cooled sphere, whereby the boundary condition is given by the surface heat flux. This heat flux depends on the temperature difference between the surface and the ambient air, such that as the sphere cools the heat flux is decreased (hence the changing boundary conditions).

I was not able to "unpack" the eq.solve function to change the boundary condition at every iteration of the solver. However, the solution I found was to put the eq.solve within a for loop and split the simulation into chunks (=100*dt). The diffusion equation is solved using constant boundary conditions within the loop, and then at the end of each iteration the final radial temperature grid is saved and the new boundary condition is calculated. In the next iteration of the loop, the same process is repeated except the field I start solving from is updated from the final field of the last loop. While this sampling of the boundary conditions is an approximation, for my project it was easily accurate enough.

Py-pde seems like a fantastic package, however it would be good to see this feature implemented at some point down the line. Hope that helps!

Code:

import pde
import numpy as np

R_p = 203.4 * 10**(-6)
t_f = 0.5
T_i = 1943 #initial melt temp in kelvin
T_e = 280 #atmospheric temperature

dt = 0.00001
res = 30 # the resolution of the radial grid
#thermal properties
k_g = 0.024 #gas conductivity
k = 1.5 #conductivity
sigma = 5.67 * 10**(-8) #steffan-boltzmann constant
p_g = 1.225 #density of air

K_p = 5.71 * 10**(-7) #thermal diffusivity of particle

#spherically symmetric diffusion equation
grid = pde.SphericalGrid(radius=[0, R_p], shape=res) #defines a grid in spherical coordinates with res points
field = pde.ScalarField(grid, data = T_i) #initial temp everywhere is T_i

T_vals = []
dt_large = dt*100
T_s = T_i #surface temp
for i in np.arange(0,t_f, dt_large): #discretize so that instead of boundary conditions changing every iteration of the pde solver, they change every 1000 steps (approximation - will result in faster cooling)

print("time = {:.2f} s".format(i*10**3))
Re = 10.27 *10**3 #reynolds number for 200 micron particle 
Pr = 0.71 #constant - from table 
N = 2 + 0.838 * Re**(1/2) * Pr**(1/3)

h = N*k_g/(2*R_p)  #heat transfer coefficient   

#radiative and convective heat flux
q_R = sigma* (T_e**4 - T_s**4) 
q_c = h*(T_e - T_s)
flux = (q_R + q_c)/k   #surface heat flux boundary condition 
print("q_R =",q_R)
print("q_c =",q_c)
print("flux =",flux)  


#solving the 1D spherically symmetric heat equation use py-pde
bc_r = [{"derivative":0},{"derivative":flux}]  #time varying boundary condition
eq = pde.DiffusionPDE(K_p,bc=bc_r)
result = eq.solve(field, t_range = dt_large, dt=dt) #returns radial temperature array after iterating 1 step only of the euler method 
T_r = result.data  #the radial temperature profile for this step
T_vals.append(T_r)
print(T_r)
T_s = T_r[-1]  #redefine the surface temperature as the last element of the radial temperature profile
print(T_s)

result.plot(kind="image",cmap='magma')
        
grid = pde.SphericalGrid(radius=[0, R_p], shape=res) #redefine grid 
field = pde.ScalarField(grid, data=result)  #input temeprature array from the previous step

print(T_vals)