Need help explaining a confusing error

Hi everyone!

Having some trouble debugging some unusual errors that are triggering on pretty innocuous looking lines that I can’t see the problems with! Would be amazing if someone could help me to interpret what the error is saying and see if maybe I can figure out why they’re throwing! I’ll paste the error here:

Traceback (most recent call last):
  File "wabisabi.py", line 87, in <module>
    solver.simulate()
  File "/home/josh/Desktop/FastIPC/projects/brittle/DFGMPMSolver.py", line 1719, in simulate
    self.substep()
  File "/home/josh/Desktop/FastIPC/projects/brittle/DFGMPMSolver.py", line 1465, in substep
    self.G2P()
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/kernel.py", line 605, in __call__
    return self._primal(self._kernel_owner, *args, **kwargs)
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/kernel.py", line 501, in __call__
    self.materialize(key=key, args=args, arg_features=arg_features)
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/kernel.py", line 370, in materialize
    taichi_kernel = taichi_kernel.define(taichi_ast_generator)
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/kernel.py", line 367, in taichi_ast_generator
    compiled()
  File "/home/josh/Desktop/FastIPC/projects/brittle/DFGMPMSolver.py", line 1118, in G2P
    self.x[p] += self.dt * new_v_PIC # advection, use PIC velocity for advection regardless of PIC, FLIP, or APIC
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/common_ops.py", line 240, in augassign
    self += x
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/common_ops.py", line 184, in __iadd__
    self.atomic_add(other)
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/common_ops.py", line 160, in atomic_add
    return ti.atomic_add(self, other)
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/ops.py", line 136, in wrapped
    return a.element_wise_writeback_binary(imp_foo, b)
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/matrix.py", line 176, in element_wise_writeback_binary
    ret.entries[i] = foo(self.entries[i], other)
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/ops.py", line 130, in imp_foo
    return foo(x, wrap_if_not_expr(y))
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/ops.py", line 31, in wrap_if_not_expr
    return Expr(a) if not is_taichi_expr(a) else a
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/expr.py", line 33, in __init__
    self.ptr = impl.make_constant_expr(arg).ptr
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/util.py", line 199, in wrapped
    return func(*args, **kwargs)
  File "/home/josh/.local/lib/python3.6/site-packages/taichi/lang/impl.py", line 280, in make_constant_expr
    return Expr(taichi_lang_core.make_const_expr_f64(val))
TypeError: make_const_expr_f64(): incompatible function arguments. The following argument types are supported:
    1. (arg0: float) -> taichi_core.Expr

Invoked with: array([<ti.Expr>, <ti.Expr>, <ti.Expr>], dtype=object)

To me, it seems the problematic line here looks to be:
self.x[p] += self.dt * new_v_PIC

And previously, I had this exact same looking error throwing on this line:
new_v_FLIP = self.v[p] + (self.dt * new_v_FLIP)

BUT, the error went away and instead became the one above when I changed this line to this:
new_v_FLIP *= self.dt
new_v_FLIP += self.v[p]

Could someone help me to interpret what might be going on here based on the error? Not sure I understand what ti.Expr is :confused:

Generally speaking, ti.Expr is a dummy class for recording operations so that we can generate C code for the kernel.

Let me show you a rough idea on how these magic works.

Fun fact: the Python code inside @ti.kernel is actually executed for only once. Later invocations into that kernel will be redirected to launch an executable binary compiled by GNU C compiler for high-performance execution. Therefore Python objects like list or np.array won’t be available, there are only basic types (int, float).

Q: So how do we “generate C code” from your Python code?
A: We use an ti.Expr object with operator overrides to replace all the int and floats in your code, so when the code is executed once, ti.Expr can record the operaions so that we get the C code generated.

When you are writing:

x = 1
y = x + 2

It is actually being transformed secretly by Taichi compiler into:

x = ti.Expr(1)
y = x.__add__(ti.Expr(2))

Then ti.Expr.__add__ will append some “C code” by something implemented like this:

def __init__(self, const_val):
  printf("int %s = %d;", self.var, const_val)

def __add__(self, other):
  out = ti.Expr()
  printf("int %s = %s + %s;", out.var, self.var, other.var)
  return out

...

Finally, the output of all these "printf"s will be dumped into a C source file, and we will invoke the gcc compiler to compile them into an executable file. When the kernel is called, it will call the coresponding C function instead of your Python function, that’s why Taichi is faster.
The same applies to CUDA, OpenGL and Metal backend, they generate GPU shader codes instead though.


For the error message, it seems that you are using np.array inside a Taichi kernel, this is not supported/allowed.
If you want vectors, consider use ti.Vector instead.

I was able to get rid of these errors by simply reorganizing some lines, could you explain to me why this would fix these errors? I changed the line
self.x[p] += self.dt * new_v_PIC

simply to:
new_v_PIC *= self.dt
self.x[p] += new_v_PIC

What might be causing the error in the first place such that this fixes it??

Could you give me the full code (or a link to your github repo) so that I can reproduce the error? Thank in advance!

What parameter did you specified to initializw the Solver? e.g.:

dt = np.float(0.1)
solver = MPMSolver(dt=dt)

And, who wrote the solver?

Hi again, and thanks for your responses! I wrote the solver myself, and it’s part of a research codebase so I don’t think I can share it yet, but here’s how I’m setting dt!

def suggestedDt(E, nu, rho, dx, cfl):
    elasticity = math.sqrt(E * (1 - nu) / ((1 + nu) * (1 - 2 * nu) * rho))
    return cfl * dx / elasticity

E = 1e4
nu = 0.15
rho = 1
ppd = 3
ppc = ppd**3
vertices, vol = sampleFromTetWild(filename, rho) //particle sampling routine
vertexCount = len(vertices)
pVol = vol / vertexCount
dx = (ppc * pVol)**0.5
cfl = 0.4
maxDt = suggestedDt(E, nu, rho, dx, cfl)
dt = 0.9 * maxDt

Then, dt is just passed into the constructor and set using:
self.dt = dt

Thank you for sharingthe code! Would you mind help me confirm the object type of dt by:

print(type(dt))
print(type(dx))
print(type(cfl))

?
If one if them is np.ndarray, it’s likely the problem.

woah this was the problem!! I used numpy.linalg.det to compute tetrahedral volume, so it turns out that volume was converting every quantity it touched to numpy.float64, including dx! Here’s the output from that:

Type of vol: <class 'numpy.float64'>
Type of Dt: <class 'float'>
Type of dx: <class 'numpy.float64'>
Type of cfl: <class 'float'>

Thanks so much for all the help!!!