Adding forcing term in fluid simulation and verifying the simulation results via residual of the pde

[hangzhou188]

Hi, I am interested in fluid simulation, and i am following the codes at https://colab.research.google.com/github/tum-pbs/PhiFlow/blob/develop/docs/Fluid_Simulation.ipynb#scrollTo=6V52HLkhqWff.
I have two questions:

1. Can we add a forcing term in the pde equation so as to simulate more complex situations?
2. After we get the simulated velocity and pressure field, can we verify the results by computing residual of the ns equation? If so, how can we do that? I have some difficulty computing the (𝑢⋅∇)𝑢 and the 𝜈∇^2 𝑢 term here.

Thank you very much! Your project truly helps me a lot.

[holl-]

1. Sure, for an example of fluid flow with forcing see here.
2. Generally yes but you need to define what exactly the residuals should measure, taking the discretization into account. In case of operator splitting, the term 𝜈∇^2 if computed by diffuse.explicit and (𝑢⋅∇)𝑢 by the chosen advect function. With higher-order spatial integration, the terms are computed within the momentum_equation function.

[hangzhou188]

Thank you a lot! I am trying to run the code you provided in the colab, and I am running into some errors here:

Traceback (most recent call last):
File "", line 1, in
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/_functional.py", line 1096, in iterate
x = f(*x, **f_kwargs)
File "", line 1, in
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/_functional.py", line 1086, in iterate
x = f(*x, **f_kwargs)
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/_functional.py", line 198, in call
native_result = self.traceskey
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/backend/torch/_torch_backend.py", line 1023, in call
self.f(*args) # records all autograd.Function / nn.Module calls with their args -> self.autograd_function_calls, self.called_modules
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/_functional.py", line 168, in jit_f_native
result = self.f(**kwargs, **in_key.auxiliary_kwargs) # Tensor or tuple/list of Tensors
File "", line 1, in rk4_step
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phi/physics/fluid.py", line 269, in incompressible_rk4
rhs_1 = pde(v_1, **pde_aux_kwargs) - field.spatial_gradient(p_1, type=StaggeredGrid, order=pressure_order)
File "", line 1, in momentum_equation
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phi/physics/diffuse.py", line 91, in finite_difference
return diffusivity * laplace(grid, order=order, implicit=implicit).with_extrapolation(grid.extrapolation)
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phi/field/_field_math.py", line 92, in laplace
result_components = solve_linear(_lhs_for_implicit_scheme, result_components, solve=implicit, values_rhs=values_rhs, needed_shifts_rhs=needed_shifts_rhs, stack_dim=channel('laplacian'))
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/_optimize.py", line 588, in solve_linear
return _matrix_solve(y - bias, solve, matrix)
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/_functional.py", line 807, in call
native_result = self.traceskey # With PyTorch + jit, this does not call forward_native every time
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/backend/torch/_torch_backend.py", line 227, in select_jit
output = torch_function.apply(*args)
File "/anaconda3/envs/torch/lib/python3.9/site-packages/torch/autograd/function.py", line 506, in apply
return super().apply(*args, **kwargs) # type: ignore[misc]
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/backend/torch/_torch_backend.py", line 1083, in forward
y = (jit_f or f)(*args, **kwargs)
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/_functional.py", line 763, in forward_native
result = self.f(**kwargs, **in_key.auxiliary_kwargs) # Tensor or tuple/list of Tensors
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/_optimize.py", line 584, in _matrix_solve_forward
result = _linear_solve_forward(y, solve, backend_matrix, pattern_dims_in, pattern_dims_out, preconditioner, backend, is_backprop)
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/_optimize.py", line 665, in _linear_solve_forward
result = SolveInfo(solve, x, residual, iterations, function_evaluations, converged, diverged, ret.method, msg, t)
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/optimize.py", line 168, in init
msg = map(_default_solve_info_msg, msg, converged.trajectory[-1], diverged.trajectory[-1], iterations.trajectory[-1], solve=solve, method=method, residual=residual, dims=converged.shape.without('trajectory'))
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/functional.py", line 1141, in map
idx_kwargs = {k: v[idx] for k, v in sliceable_kwargs.items()}
File "/anaconda3/envs/torch/lib/python3.9/site-packages/phiml/math/_functional.py", line 1141, in
idx_kwargs = {k: v[idx] for k, v in sliceable_kwargs.items()}
File "/anaconda3/envs/torch/lib/pythonz3.9/site-packages/phi/field/_grid.py", line 423, in getitem
raise AssertionError("Cannot slice StaggeredGrid along spatial dimensions.")
AssertionError: Cannot slice StaggeredGrid along spatial dimensions.

The above error occurred when running the code v_trj, p_trj = iterate(multi_step, batch(time=100), v0, p0, dt=0.005, range=trange). Would you please help me solve this problem?

[holl-]

I can run the Higher-order Fluid Simulations notebook in Colab without error.

You are talking about this cell, right?

v0 = StaggeredGrid(0, **DOMAIN)
p0 = CenteredGrid(0, **DOMAIN)
multi_step = lambda *x, **kwargs: iterate(rk4_step, 25, *x, **kwargs)
v_trj, p_trj = iterate(multi_step, batch(time=2), v0, p0, dt=0.005, range=trange)
vis.plot(field.curl(v_trj.with_extrapolation(0)), animate='time')

That cell should show a progress bar. Does that show up?

Could you try running

!pip uninstall phiflow phiml
!pip install phiflow

and restarting your runtime?

[hangzhou188]

Hi, when I ran the code with default environment, i.e. jax and cpu, I got no error here. However, when I tried to replace the first line with from phi.torch.flow import * and tried to run the simulation with cuda, I got the above error.

[holl-]

I'll look into it.

The higher-order solvers are still under development and at the moment only the Jax version is fully tested. So if you can switch to Jax, I'd recommend using it.