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.traces[key](natives) # 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.