Hello fellow friends of PySPH,
i've been trying to implement Monaghan Style Heat Transfer in PySPH but ran into a few Problems. I have probably misunderstood something, and would be very thankful for any hints.
I've based my Code on this Equation:
It ends up looking something like this:
```Python
class HeatConduction(Equation):
def __init__(self, dest, sources):
super(HeatConduction, self).__init__(dest, sources)
def loop(self, d_idx, s_idx, d_T, s_T, DWIJ, d_dT, s_m, d_rho, s_rho, R2IJ, XIJ):
# i ^= d; j ^= s
dot_product = DWIJ[0] * XIJ[0] + DWIJ[1] * XIJ[1] + DWIJ[2] * XIJ[2]
if d_idx == s_idx:
d_dT[d_idx] += 0
else:
d_dT[d_idx] += ((s_m[s_idx] / (d_rho[d_idx] * s_rho[s_idx]))
* ((4 * 4500 * 4500) / (4500 * 4500))
* (d_T[d_idx] - s_T[s_idx])
* (dot_product / R2IJ))
def converged(self):
return 1 # assume convergence for now
```
With the acccording Integrator looking something like:
```Python
class HeatStep(IntegratorStep):
def initialize(self, d_idx, d_T):
pass
def stage1(self, d_idx, d_T, d_dT, dt):
pass
def stage2(self, d_idx, d_T, d_dT, dt):
d_T[d_idx] += dt * d_dT[d_idx]
```
I've included the full .py file as an attachment as well. It will run and visualize the simulation.
For the first few steps, it looks like what i would expect:
As you can see in the Image above, the rate of change of the Temperature "dT" after 10 timesteps is positive on the initially cold (upper) side and negative on the warm side, so the temperature should equalize (and it does for the first 8 steps or so).
However, after ~10 steps, the simulation becomes very unstable and the temperature and its rate of change start to oscillate back and forth between the two sides like this:
This Problem of waves of Heat travelling up and down can be seen more clearly in this 2d example:
I believe i must have made a mistake implementing the stepper or equation, but am neither smart nor experienced enough to figure out where.
Any hint in the right direction is highly appreciated!
Thank you in advance
Max