Time step is too strict in mpm88.py

[f1shel]

Describe the bug
Fine tuning the time step can cause program errors in mpm88.py

To Reproduce
Change code in mpm88.py:L9 from dt = 2e-4 to dt = 2.8e-4 will raise memory access errors.

Log/Screenshots

$python .\mpm88.py
[Taichi] version 1.6.0, llvm 15.0.1, commit f1c6fbbd, win, python 3.9.0
[Taichi] Starting on arch=cuda
...
RuntimeError: [taichi/rhi/cuda/cuda_driver.h:taichi::lang::CUDADriverFunction<void *,void *,unsigned __int64>::operator ()@92] CUDA Error CUDA_ERROR_ILLEGAL_ADDRESS: an illegal memory access was encountered while calling memcpy_host_to_device (cuMemcpyHtoD_v2)
[E 09/04/23 11:56:41.537 25652] [taichi/rhi/cuda/cuda_driver.h:taichi::lang::CUDADriverFunction<void *>::operator ()@92] CUDA Error CUDA_ERROR_ILLEGAL_ADDRESS: an illegal memory access was encountered while calling stream_synchronize (cuStreamSynchronize)

Additional comments
Perhaps the error is raised by illegal array indexing in both G2P and P2G part, since none boundary constraints is applied on particles and grids. To be specific:

# particle to grid
for i, j in ti.static(ti.ndrange(3, 3)):
    offset = ti.Vector([i, j])
    dpos = (offset - fx) * dx
    weight = w[i].x * w[j].y
    grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos)
    grid_m[base + offset] += weight * p_mass

for i, j in ti.static(ti.ndrange(3, 3)):
    offset = ti.Vector([i, j])
    dpos = (offset - fx) * dx
    weight = w[i].x * w[j].y
    g_v = grid_v[base + offset]

will throw out error when indexing base + offset is out of range. This is caused by particles that have wrongly moved inside the boundary predefined in L15 bound=3. Simply clamping particle position can get rid of illegal memory access, but still result in exploded simulation results due to negative J[p].

@ti.func
def clamp(val, l, r):
    return ti.min(ti.max(val, l), r)
...
		x[p] += dt * v[p]
		x[p] = clamp(x[p], (bound + 0.0) * dx, (n_grid - bound - 0.0)*dx)
        J[p] *= 1 + dt * new_C.trace()
        if J[p] <= 0:
            print(J[p])

[jim19930609]

Was able to reproduce locally with taichi-nightly. Seems that whenever dt gets larger than a threshold, it's gonna cause an error with CUDA.

Let me investigate that further.

[jim19930609]

The thing is that the value of dt will essentially influence the index of certain memory access:

When dt exceeds certain threshold, this indexing will get out-of-bound thus result in a CUDA illegal address error.

So it is expected behaviour though. The best you can improve, if any, would be to add a boundary for stress but that will degrade performance as well.

[f1shel]

@jim19930609 Discovering a new problem in mpm88, probably I should open a new issue, but I will post it here though.

According to slides of GAMES201 Lec08, constitutive model of water in mpm88 is defined as (page 17)

The internal force on the node is given by (page 11)

The cauchy stress sigma and PK1 stress P are related by

So the force can be written as

In mpm88, the force momentum is given in the following code

        stress = -dt * 4 * E * p_vol * (J[p] - 1) / dx**2
        affine = ti.Matrix([[stress, 0], [0, stress]]) + p_mass * C[p]

The proper form should consist of a extra J term due to conversion from P to sigma

        stress = -dt * 4 * E * J[p] * p_vol * (J[p] - 1) / dx**2
        affine = ti.Matrix([[stress, 0], [0, stress]]) + p_mass * C[p]

[Elijesan]

Hello, I am trying to do a simulation (SPH), but when I use a dt varibele (courant condition) I get that there is an error when accessing the memory, it only happens when I put the dt variable, thare is a bug with taichi ?