PieterjanRobbe / GaussianRandomFields.jl

A package for Gaussian random field generation in Julia
Other
66 stars 11 forks source link

Unstructured grid without elements #34

Closed gerrygralton closed 3 years ago

gerrygralton commented 3 years ago

Currently, you cannot use an unstructured grid with a non-uniform set of points unless it comes with connectivity (elements). This seems like an unnecessary complication.

My current workaround for a 3D unstructured grid is, grf = GaussianRandomField(cov_fun, KarhunenLoeve(1000), coords, Matrix{Int64}(undef, 0, 4)) but I think this should be supported.

PieterjanRobbe commented 3 years ago

Fixed by #36. Here's an example of a random field on a set of points inside a sphere:

n = 1024
x = 2*rand(n) .- 1
y = 2*rand(n) .- 1
z = 2*rand(n) .- 1
d = x.*x + y.*y + z.*z

pts = hcat(x, y, z)
pts = pts[d .< 1, :]
grf = GaussianRandomField(CovarianceFunction(3, Matern(1, 3/2)), KarhunenLoeve(128), pts)
scatter(pts[:, 1], pts[:, 2], pts[:, 3], marker_z=sample(grf), label="")

See grf_unstructured_grid.pdf.

PieterjanRobbe commented 3 years ago

If you're using three-dimensional random fields you might consider using circulant embedding instead of a KL expansion? Depending on the smoothness of the random field, it may be (much) faster to generate the random field on a structured mesh with a sufficiently small grid size, and resort to linear interpolation. That would depend on the application you have in mind, of course.

gerrygralton commented 3 years ago

Thank you for responding so quickly, it's much appreciated!

  1. The fix works perfectly for KL expansions because it passes the quadrature nodes to the.apply() function. With spectral or Cholesky methods it tries to pass the pts and the empty element array and there's no apply() function to handle that.
  2. I'm probably mistaken but I thought this kind of application was the primary use case for the KL expansion - to avoid having to assemble a full NxN matrix and being able to interpolate back to the points of interest. Why is a circulant embedding so advantageous?
PieterjanRobbe commented 3 years ago

Ah yes good point, that definition was still missing. Can you confirm that it's working now for Cholesky and Spectral in three dimensions? (See #37)

gerrygralton commented 3 years ago

Yep, that's working now! Cheers

PieterjanRobbe commented 3 years ago

As for the second point, in my experience, the KL expansion becomes expensive when you have a small correlation length and low smoothness (many KL terms) and many grid points. Circulant embedding may be (much!) cheaper in this setting. That being said, it all depends on what correlation lengths/smoothness and number of grid points you need in your application.

gerrygralton commented 3 years ago

OK, thanks for the tip. I'll have a play and see what works best.