tum-pbs / PhiFlow

A differentiable PDE solving framework for machine learning
MIT License
1.39k stars 189 forks source link

How to Insert center values into staggerGrid class? #151

Closed XueruiSu closed 1 month ago

XueruiSu commented 6 months ago

When I try the karman vortex case, I want to calculate from a given ndarray. such as: velocity_field = np.random.rand(128, 128, 2) velocity_repeat = np.repeat(velocity_field[np.newaxis, :, :, :], repeats=10, axis=0)

I take the velocity_repeat array as the center grid values. How could I insert the velocity_repeat array into StaggeredGrid class? I tried the CenteredGrid class but it turned out terrible numerical error.

holl- commented 6 months ago

Hi @XueruiSu , You can convert a CenteredGrid into a StaggeredGrid, like in the following example

centered = CenteredGrid(Noise(), PERIODIC, Box(x=1, y=1), x=32, y=32)
staggered = StaggeredGrid(centered, centered.boundary, centered.bounds, centered.resolution)
XueruiSu commented 6 months ago

I used your example and It raised an error: "ValueError: <class 'phiml.math.magic.BoundDim'> is not supported. Only (Tensor, tuple, list, np.ndarray, native tensors) are allowed. Current backends: [numpy, Python]" Can the first parameter of StaggeredGrid not be the CenteredGrid object?

holl- commented 6 months ago

Sorry, it should be

centered = CenteredGrid(Noise(), PERIODIC, Box(x=1, y=1), x=32, y=32)
staggered = centered.at(StaggeredGrid(0, centered.extrapolation, centered.bounds, centered.resolution))
XueruiSu commented 6 months ago

Thanks for your patient reply. Then I wrote these code to check the value equivalence between CenteredGrid and StaggeredGrid. But I found that the value would change a lot when I extracted the center values from the StaggeredGrid which been inserted by CenteredGrid. The mse loss is about 0.35 for value from the uniform(-1, 1) distribution. Is it because the interpolation method I chose when doing the conversion is incorrect?

from phi.torch.flow import *
import numpy as np
res_x = 128
res_y = 128
res_x_real, res_y_real = 96, 64
rotation_speed = 0
speed = 2.
radius = 5
SPEED = vis.control(speed)
velocity_repeat = np.random.uniform(-1, 1, (10, res_x, res_y, 2))
velocity_CenteredGrid = CenteredGrid(
    values=math.tensor(velocity_repeat, batch('inflow_loc'), spatial('x,y'), channel(vector='x,y')),
    extrapolation=ZERO_GRADIENT,
    x=res_x, y=res_y,
    bounds=Box(x=res_x_real, y=res_y_real),
)
velocity_StaggeredGrid = velocity_CenteredGrid.at(StaggeredGrid(0, velocity_CenteredGrid.extrapolation, velocity_CenteredGrid.bounds, velocity_CenteredGrid.resolution))
velocity_StaggeredGrid_extracted = velocity_StaggeredGrid.at_centers().data.numpy('inflow_loc,y,x,vector')
velocity_CenteredGrid_extracted = velocity_CenteredGrid.values.numpy('inflow_loc,y,x,vector')
print(velocity_StaggeredGrid_extracted[0, 0, 0, 0], velocity_CenteredGrid_extracted[0, 0, 0, 0])
# -0.43487194 -0.53899693
print(np.sqrt(((velocity_StaggeredGrid_extracted-velocity_CenteredGrid_extracted)**2).mean()))
# 0.35093784
holl- commented 6 months ago

Hi, the interpolation always has a smoothing effect. You can try 6th order implicit interpolation.

for interp in [
    dict(order=2, implicit=None),
    dict(order=6, implicit=Solve('biCG-stab(2)')),
    ]:
  centered = CenteredGrid(Noise(), PERIODIC, Box(x=1, y=1), x=32, y=32)
  staggered = centered.at(StaggeredGrid(0, centered.extrapolation, centered.bounds, centered.resolution), **interp)
  centered2 = staggered.at(centered, **interp)
  print(interp, "err:", field.l2_loss(centered - centered2))
  plot(centered, centered2['x'])