- #1
Telemachus
- 835
- 30
Hi. I was trying to test a code for the diffusion equation, using the analytical solution for infinite media, with a Dirac delta source term:
##q(\mathbf{r},t)=\delta (\mathbf{r}) \delta (t)##.
The code is not giving the analytical solution, and there might be several reasons why this is so. One reason could be that my code has something wrong, a bug, something I did wrong. Another reason might be that my modelling of the analytic situation is incorrect, for example, the boundary conditions might be affecting my solution (my numerical model is of course not infinite).
I am solving in a cube, and looking at the solution at the center of the domain in order to avoid border effects.
The equation I am solving is: ##\displaystyle \frac{1}{c}\frac{\partial \phi(\mathbf{r},t)}{\partial t}-D\nabla^2 \phi(\mathbf{r},t)=q(\mathbf{r},t)##.
I have tried with a cube of side two, and of side one. And I plot the solution at a distance of 0.2 units from the source.
So, to model the Dirac delta source I was using:
There is also, of course, a loop in time.
However, I have noted that this is very sensitive to the sigma values I use. I am also using a few points in the spatial coordinates, but I don't think that is an issue (perhaps it could be affecting the distance I am considering to the source, but not more than that, the solution looks converged, and if I use many more points I don't get considerable changes).
I am a bit worried about this, because I thought the code was working fine, the solutions looks physically as I was specting, and I get curves similar to the one reported in the papers, but when I've tried to see in detail how it behaves with respect to an analytic solution, I've found with this issue.
So first, I would like to know if my modelling of the Dirac delta function is appropriate. Second, how big should be my domain in order to avoid any effects from the boundaries? I think that, due to the small distance to the source, and considering the speed of propagation of the wave, the borders shouldn't be bothering. Should I look my problem somewhere else in the code? or should I focus in the proper modelling of the analytic situation, with an infinite homogeneous media, and the Dirac delta source term?
A possible way to model a one dimensional Dirac delta functions is to use that:
##\delta (x)=\displaystyle \lim_{\sigma \to 0} \frac{e^{-x^2/\sigma^2}}{\sigma\sqrt{\pi}}##
For a 3D dirac delta, with an impulse in time, I thought of using:
##\delta (\mathbf{r})\delta (t)=\displaystyle \lim_{\sigma \to 0} \frac{e^{-[(x-x̣_0)^2+(y-y_0)^2+(z-z_0)^2+(t-t_0)^2]/\sigma^2}}{\sigma^4 \pi^2}##
However, dividing for a small sigma seemed to cause trouble. I am not totally sure about this. But I think that the modelling of the Dirac delta function causes some problem, and that might be the reason why I am not being able to get the analytical situation. When I do this numerically, I take a small ##\sigma##, and I've noted that small variations on it produces big changes in the solution.
##q(\mathbf{r},t)=\delta (\mathbf{r}) \delta (t)##.
The code is not giving the analytical solution, and there might be several reasons why this is so. One reason could be that my code has something wrong, a bug, something I did wrong. Another reason might be that my modelling of the analytic situation is incorrect, for example, the boundary conditions might be affecting my solution (my numerical model is of course not infinite).
I am solving in a cube, and looking at the solution at the center of the domain in order to avoid border effects.
The equation I am solving is: ##\displaystyle \frac{1}{c}\frac{\partial \phi(\mathbf{r},t)}{\partial t}-D\nabla^2 \phi(\mathbf{r},t)=q(\mathbf{r},t)##.
I have tried with a cube of side two, and of side one. And I plot the solution at a distance of 0.2 units from the source.
So, to model the Dirac delta source I was using:
Fortran:
sigg=0.1d0
siggan=0.01d0
Power=20.d0
rxs=0.5d0 ! x position of source
rys=0.d0 ! y position of source
rzs=0.5d0 ! z position of source
rts=0.3d0 ! at this time the pulse is emitted
do k=1,zrows
do j=1,yrows
do i=1,xrows
rq(i,j,k)= Power*dexp(-((x-rxs)**2 +
+ (y-rys)**2+(z-rzs)**2)/sigg**2)*dexp(-(t-rts)**2/siggan**2)
enddo
enddo
enddo
There is also, of course, a loop in time.
However, I have noted that this is very sensitive to the sigma values I use. I am also using a few points in the spatial coordinates, but I don't think that is an issue (perhaps it could be affecting the distance I am considering to the source, but not more than that, the solution looks converged, and if I use many more points I don't get considerable changes).
I am a bit worried about this, because I thought the code was working fine, the solutions looks physically as I was specting, and I get curves similar to the one reported in the papers, but when I've tried to see in detail how it behaves with respect to an analytic solution, I've found with this issue.
So first, I would like to know if my modelling of the Dirac delta function is appropriate. Second, how big should be my domain in order to avoid any effects from the boundaries? I think that, due to the small distance to the source, and considering the speed of propagation of the wave, the borders shouldn't be bothering. Should I look my problem somewhere else in the code? or should I focus in the proper modelling of the analytic situation, with an infinite homogeneous media, and the Dirac delta source term?
A possible way to model a one dimensional Dirac delta functions is to use that:
##\delta (x)=\displaystyle \lim_{\sigma \to 0} \frac{e^{-x^2/\sigma^2}}{\sigma\sqrt{\pi}}##
For a 3D dirac delta, with an impulse in time, I thought of using:
##\delta (\mathbf{r})\delta (t)=\displaystyle \lim_{\sigma \to 0} \frac{e^{-[(x-x̣_0)^2+(y-y_0)^2+(z-z_0)^2+(t-t_0)^2]/\sigma^2}}{\sigma^4 \pi^2}##
However, dividing for a small sigma seemed to cause trouble. I am not totally sure about this. But I think that the modelling of the Dirac delta function causes some problem, and that might be the reason why I am not being able to get the analytical situation. When I do this numerically, I take a small ##\sigma##, and I've noted that small variations on it produces big changes in the solution.
Last edited: