weymouth / waterlily.jl Goto Github PK

View Code? Open in Web Editor NEW

554.0 554.0 68.0 8.25 MB

Fast and simple fluid simulator in Julia

License: Other

Julia 97.24% Shell 2.76%

waterlily.jl's People

Contributors

Stargazers

Watchers

waterlily.jl's Issues

How to use separate map functions for different bodies?

Respected sir,

I am facing an issue while using this code. Actually, I want to use separate map functions for each of the separate bodies. Could you please help me with how to do it or refer me to any other resources for the same?

Thank you for your time and consideration.

Yours sincerely,
Kunal Ghosh

Courant number (CFL)

Hi,
I have an issue with having too much noise for the force signal acting on a body. I want to know how can I change/optimize the CFL number here [https://github.com/weymouth/WaterLily.jl/blob/master/src/Flow.jl#L165] ? ? what are these 0.5 showing in these lines:

    BDIM!(a); scale_u!(a,0.5); BC!(a.u,U,a.exitBC,a.perdir)
    project!(a,b,0.5); BC!(a.u,U,a.exitBC,a.perdir)

Should the package have a .jl name?

I found the package from mentions on twitter. Really cool. Wondering if it needs a .jl extension in the package name?

Error on running the example in Youtube video.

I am running the example code in Youtube video: (https://www.youtube.com/watch?v=OO0Nz3uWDc8), and my code looks like:

And after I run this code , the error info is:

What is the reason?

Docs are pointing to wrong/old repo URL

The docs are pointing to an old repo URL.
This can be seen when trying to click on "source code" for any of the methods linked in the docs, and it'll point to a 404 page.
This can also be seen when trying to click on "Edit on Github" link at the top of the page of the docs

Add support for objects specified by NURBS?

I've worked through the dogfish example, which produces an sdf from a B-spline and is very nice. However, it's my understanding that the sdf generating function doesn't measure closest distance but are just conformal maps in the y-direction? Will this cause some numerical error when calculating the gradients? The projection (and closest distance) of any point onto the curve should be orthogonal to the curve.

If the order and knot vector of the NURBS or B-spline is such that second derivatives exist, then equation 6.3 in the NURBS Book 2nd ed by Piegl and Tiller allows calculation of the sdf. The downside is that you need to check every span. Additionally, if the curve is not closed or C2 continuous, it will make everything so much more complicated. That said, the continuity can be simply read off the knot vector. Additionally, it makes physical sense that the head and tail of the dogfish are C0 continuous.

Anyway, I think this more of a topic for discussion rather than an action item.

I've searched for Python, Julia and C++ tools to calculate the distance of a point from a NURBS curve and haven't found anything on github yet. Of the NURBS packages in Julia, the BasicBSpline package is in active development and is pretty well organized. It's just one person though and while off to a good start, there's still a lot left to do. Python-NURBS is really good for curve fitting.

Change default time in force integration

In the pressure force integration the default is to use t=0 for the body position, but this doesn't correspond to the pressure, which is stored at time(sim).

I think the default behavior should be to measure the body at time(sim), where the pressure is also defined. This would require passing sim or flow directly to the function.

ERROR: LoadError: importing WaterLily into Main conflicts with an existing identifier

Hello,
this is my first Julia Repo I am trying to run, so maybe it is a beginners Julia problem. However I want to ask this here as well, if I may.

When trying to run a Script from folder examples I get the mentioned Error. Does someone know the solution to this?
Thank you!

Heaving Cylinder in still water

Hello,
Thanks a lot for the work done for this package, it is very intuitive and easy to compute three-dimensional flows on my laptop.

For my work, I need to assess the inertial forces on a heaving cylinder in still water for a large number of conditions (hence my interest in fast simulations!). I could easily realise simulations of a moving cylinder with a constant inflow velocity.
However, I encountered two main problems when attempting to do it in still water:

1. I get confused with the non-dimensional time as it seems to rely on the inlet fluid velocity. When the inflow velocity is set to 0, the simulation gets stuck at the first time-step. How are we meant to progress the simulation forward in time with no inflow velocity?
1. Additionally, is it possible to set all the domain boundaries as outlets?

Many thanks for your help!

Integration with Flux

Thank you for this nice package!

I am trying to use Flux.jl together with a simulation, but Zygote complains with
Mutating arrays is not supported -- called copyto!(Array{Float32, 3}, ...)

Thank you for your help

A toy example is:

using WaterLily, Flux

sim = Simulation((32, 32), (0, 0), 5; U=1, ν=1e-4, T=Float32)
model = Conv((3,3), 2 => 2; pad=SamePad() )

function comploss(model, sim) 
    t₀ = round(sim_time(sim))
    loss = 0.0
    for tᵢ in range(t₀, t₀+1; step=0.1)
        sim_step!(sim,tᵢ; remeasure=false)
        sim.flow.u .+= model(sim.flow.u[:,:,:,:])[:,:,:]
        loss += sum(sim.flow.u.^2)
    end
    return loss
end

grad = gradient(m -> comploss(m, sim), model)

examples/TwoD_circle.jl failt with julia 1.4.1

example/TwoD_circle.jl fails with julia 1.4.1:

/home/lisergey/WaterLily.jl
$ JULIA_LOAD_PATH=.: julia examples/TwoD_circle.jl
ERROR: LoadError: LoadError: syntax: invalid keyword argument syntax "ϵ"
Stacktrace:
 [1] top-level scope at /home/lisergey/WaterLily.jl/src/Body.jl:27
 [2] include(::Module, ::String) at ./Base.jl:377
 [3] include(::String) at /home/lisergey/WaterLily.jl/src/WaterLily.jl:1
 [4] top-level scope at /home/lisergey/WaterLily.jl/src/WaterLily.jl:15
 [5] include(::Module, ::String) at ./Base.jl:377
 [6] top-level scope at none:2
 [7] eval at ./boot.jl:331 [inlined]
 [8] eval(::Expr) at ./client.jl:449
 [9] top-level scope at ./none:3
in expression starting at /home/lisergey/WaterLily.jl/src/Body.jl:27
in expression starting at /home/lisergey/WaterLily.jl/src/WaterLily.jl:15
ERROR: LoadError: Failed to precompile WaterLily [ed894a53-35f9-47f1-b17f-85db9237eebd] to /home/lisergey/.julia/compiled/v1.4/WaterLily/9Odjf_w5bnh.ji.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] compilecache(::Base.PkgId, ::String) at ./loading.jl:1272
 [3] _require(::Base.PkgId) at ./loading.jl:1029
 [4] require(::Base.PkgId) at ./loading.jl:927
 [5] require(::Module, ::Symbol) at ./loading.jl:922
 [6] include(::Module, ::String) at ./Base.jl:377
 [7] exec_options(::Base.JLOptions) at ./client.jl:288
 [8] _start() at ./client.jl:484
in expression starting at /home/lisergey/WaterLily.jl/examples/TwoD_circle.jl:1

The example works with 1.6.0, but 1.4.1 is a version from current Ubuntu.

$ julia --version
julia version 1.4.1
$ cat /etc/lsb-release
DISTRIB_ID=Ubuntu
DISTRIB_RELEASE=20.04
DISTRIB_CODENAME=focal
DISTRIB_DESCRIPTION="Ubuntu 20.04.2 LTS"

examples/TwoD_block.jl fails

examples/TwoD_block.jl does not run for me. Here is what I do

$ julia -e 'using Pkg; Pkg.add("Plots")'
$ JULIA_LOAD_PATH=src: julia -e 'include("examples/TwoD_block.jl"); TwoD_block_video(6)'
(L, ν) = (8.0, 0.032)
ERROR: MethodError: no method matching CartesianIndex{2}(::CartesianIndex{2})
Closest candidates are:
  CartesianIndex{2}() where N at multidimensional.jl:73
  CartesianIndex{2}(!Matched::Integer...) where N at multidimensional.jl:69
  CartesianIndex{2}(!Matched::Tuple{Vararg{Integer,N}}) where N at multidimensional.jl:64
  ...
Stacktrace:
 [1] oneunit(::CartesianIndex{2}) at ./number.jl:320
 [2] near(::CartesianIndex{2}, ::Int64) at /home/lisergey/src/WaterLily/src/GMG.jl:1
 [3] restrictPS(::Array{Float64,3}) at /home/lisergey/src/WaterLily/src/GMG.jl:7
 [4] MultiLevelPS(::Array{Float64,3}) at /home/lisergey/src/WaterLily/src/GMG.jl:29
 [5] TwoD_block_video(::Int64, ::Int64) at /home/lisergey/src/WaterLily/examples/TwoD_block.jl:18
 [6] TwoD_block_video(::Int64) at /home/lisergey/src/WaterLily/examples/TwoD_block.jl:5
 [7] top-level scope at none:0

$ julia --version
julia version 1.0.3
$ lsb_release -da
No LSB modules are available.
Distributor ID:	Debian
Description:	Debian GNU/Linux 10 (buster)
Release:	10
Codename:	buster

StackOverflowError

I am running the ThreeD_donut.jl example without using GPU and meet the error of StackOverflowError.

The code is as follows:

using WaterLily
using StaticArrays
function donut(p=6;Re=1e3,U=1)
    # Define simulation size, geometry dimensions, viscosity
    n = 2^p
    center,R,r = SA[n/2,n/2,n/2], n/4, n/16
    ν = U*R/Re

    # Apply signed distance function for a torus
    norm2(x) = √sum(abs2,x)
    body = AutoBody() do xyz,t
        x,y,z = xyz - center
        norm2(SA[x,norm2(SA[y,z])-R])-r
    end

    # Initialize simulation and return center for flow viz
    Simulation((2n,n,n),(U,0,0),R;ν,body),center
end

# import CUDA
# @assert CUDA.functional()
# sim,center = donut(mem=CUDA.CuArray);
sim,center = donut(); # if you don't have a CUDA GPU
sim_step!(sim, 0.2)  
# GLMakie.mesh(body_mesh(sim), color=:green)

dat = sim.flow.σ[inside(sim.flow.σ)] |> Array;
function ω_θ!(dat,sim,center=center)
    dt, a = sim.L/sim.U, sim.flow.σ
    @inside a[I] = WaterLily.ω_θ(I,(1,0,0),center,sim.flow.u)*dt
    copyto!(dat,a[inside(a)]) 
end

include("ThreeD_Plots.jl")
@time makie_video!(sim,dat,ω_θ!,name="donut.mp4",duration=1,step=0.25) do obs
    contour(obs, levels=[-5,5], colormap=:balance)
end

================================

After I run this code, the result in terminal as as follows:

simulation 0% complete
ERROR: StackOverflowError:
Stacktrace:
 [1] setindex!(A::GLMakie.GLAbstraction.Texture{ColorTypes.RGBA{Float32}, 1}, value::Vector{ColorTypes.RGBA{Float32}}, indexes::Base.Slice{Base.OneTo{Int64}})
   @ GLMakie.GLAbstraction ~/.julia/packages/GLMakie/9nid7/src/GLAbstraction/AbstractGPUArray.jl:38
 [2] setindex!(A::GLMakie.GLAbstraction.Texture{ColorTypes.RGBA{Float32}, 1}, value::Vector{ColorTypes.RGBA{Float32}}, indexes::Base.Slice{Base.OneTo{Int64}}) (repeats 79983 times)
   @ GLMakie.GLAbstraction ~/.julia/packages/GLMakie/9nid7/src/GLAbstraction/AbstractGPUArray.jl:41

But the torus could be displayed if I only run the following line instead of plotting the simulation results:

GLMakie.mesh(body_mesh(sim), color=:green)
PS:
My laptop is mac M2. Julia version is 1.9.2. And Version of packages is as follows:
(WaterLily.jl) pkg> status

Status `~/Work/research/WaterLily.jl/Project.toml`
  [e9467ef8] GLMakie v0.8.7
  [5c1252a2] GeometryBasics v0.4.9
  [e6723b4c] Meshing v0.6.0
⌃ [dde4c033] Metal v0.3.0
  [91a5bcdd] Plots v1.38.17
  [90137ffa] StaticArrays v1.6.2
  [ed894a53] WaterLily v1.0.2

MultiLevelPoisson solver tolerance not independent of domain size

I just realized (while looking at pressure forces) that the pressure field oscillates when the domain is large (384,256), see attached script.

I have managed to track the issue to the MultiLevelPoisson/solver! function.
This is linked to the V-cycle criterion that checks if the residual on the upper level (r₂) is below a certain tolerance (tol) before performing the Vcycle

while r₂>tol && nᵖ<itmx
      Vcycle!(ml)
      smooth!(p); r₂ = L₂(p)
      log && push!(res,r₂)
      nᵖ+=1
end

The issue is that (r₂) depends on the grid size through p.r as L₂(p::Poisson) = p.r ⋅ p.r, while tol does not. This means that in some cases, the solver! doesn't perform any V-cycle, which doesn't change the pressure field.

The proposed solution is to scale the tol by the grid size, compute it once and store in in the MultiLevelPoisson structure, see #97
TwoD_circle_pressure_check.zip

Solving a passive scalar conservation equation

Hoping to add functionality to an (initially) passive scalar calculation of N scalars to the WaterLily flow solver. I'd label this as a want but not clear how to do it through the website.

Question about usage

I understand pretty much how to define a body using just the distance function. Makes a lot of sense. I have two major questions, however.

First, what is the impact of small errors in this function as long as the distance at the surface of the body is zero and the errors are relatively small? For instance, if I have an elliptical body and calculate distance to the surface using the distance to a surface point on a line between the point in question and the center of the ellipse, I will give a distance that is a bit longer than the actual shortest distance. What would be the effect of that? The same thing happens in the simulation of the swimming fish. If I define the fish body as the distance to the nearest point on the spine less the orthogonal thickness at that point, I will give a slightly higher distance than is geometrically correct.

Second, can I use the map function in a body to project a circle to a wing section and then just compute the distance in the projected space?

Modeling flow over a full 3D normal facing disk

Following along with your YouTube example for various 3D shapes (https://youtu.be/OO0Nz3uWDc8?si=FGseNtBGkDdYK05p), when changing the "map" function for the corner disk to place it in the center of the domain and change it to a full disk, the model does not seem to converge. I have tried increasing the resolution, the size of the domain, the size and thickness of the disk, and it's position relative to the inlet and outlet and the same issues persist. More specifically, the errors given show that the lambda2 data is a non-hermitian matrix.

original working corner disk

function make_sim(;n=2^6, U=1, Re=1000, mem=Array)
    D = n/2
    R = D/2

    norm2(x) = √sum(abs2,x)
    function sdf(xyz,t)
        x,y,z = xyz
        r = norm2(SA[y,z]);       
        norm2(SA[x,r-min(r,R)]) - 1.5  
    end
    map(xyz,t) = xyz - SA[2n/3,0,0]

    return Simulation((2n,n,n), (U,0,0), D; 
        ν=U*D/Re, body=AutoBody(sdf,map), mem)
end

change to full centralized disk

map(xyz,t) = xyz - SA[2n/3,n/2,n/2]

Similarly, the same issues show to persist if the corner disk is moved to

map(xyz,t) = xyz - SA[n/2,0,0]

Circle example not working

Default circle function does not work:

function circle(n,m;Re=250,U=1)
    radius, center = m/8, m/2
    body = AutoBody((x,t)->√sum(abs2, x .- center) - radius)
    Simulation((n,m), (U,0), radius; ν=U*radius/Re, body)
end

circ = circle(3*2^6,2^7)
t_end = 10
sim_step!(circ,t_end)

Yields the following error:

MethodError: no method matching WaterLily.Simulation(::Tuple{Int64, Int64}, ::Tuple{Int64, Int64}, ::Float64; ν=0.064, body=WaterLily.AutoBody{WaterLily.var"#comp#45"{Main.var"workspace#37".var"#2#3"{Float64, Float64}, WaterLily.var"#43#44"}, WaterLily.var"#43#44"}(WaterLily.var"#comp#45"{Main.var"workspace#37".var"#2#3"{Float64, Float64}, WaterLily.var"#43#44"}(Main.var"workspace#37".var"#2#3"{Float64, Float64}(64.0, 16.0), WaterLily.var"#43#44"()), WaterLily.var"#43#44"()))

Closest candidates are:

WaterLily.Simulation(::Tuple, !Matched::Vector, ::Number; Δt, ν, U, ϵ, uλ, body, T) at ~/.julia/packages/WaterLily/1t0mb/src/WaterLily.jl:60

var"#circle#1"(::Int64, ::Int64, ::typeof(Main.var"workspace#37".circle), ::Int64, ::Int64)@[Other: 4](http://localhost:1234/edit?id=0cb9c7e4-e058-11ee-2ff4-b5a8d4e4ecb5#)
circle(::Int64, ::Int64)@[Other: 1](http://localhost:1234/edit?id=0cb9c7e4-e058-11ee-2ff4-b5a8d4e4ecb5#)
top-level scope@[Local: 2](http://localhost:1234/edit?id=0cb9c7e4-e058-11ee-2ff4-b5a8d4e4ecb5#)[inlined]

examples ThreeD_Plots.jl

Delete line 2 and line 3. There is no need to import AbstractPlotting and do AbstractPlotting.inline!(true) because AbstractPlotting.jl has already been fully merged into Makie.jl. Please read Readme.md in https://github.com/JuliaPlots/AbstractPlotting.jl for more details.

Simulating two different fluids

Hi, I found this repo and love the library. I think that adding moving objects this easily is great. I have been using WaterLily for some time.

My new project involves a smoke simulation with a moving object. Therefore I was wondering, can WaterLily simulate 2 different fluids? especially because densities become variable on the domain and buoyancy forces become important.

The details of the simulation: Imagine a pipe full of smoke, then a piston accelerates the smoke out of the pipe into the air forming a smoke ring. Like in the donut example. I want to study how the ring changes with the simulation parameters. Smoke is much denser than air, that's my biggest problem.

Another approach I thought of was to add a lot of very small particles and update their position using the get force function and then add the buoyancy force but I don't know how good this approach is.

What do you recommend?

I could of course just study the vorticity of the air and don't mess with two different fluids, but I was wondering if WaterLily could do the job.

Thanks.

Docs "source" buttons do not go to source code for relevant function

As mentioned in #21 (comment)

In the deployed docs, the "source" buttons that appear when you hover over functions' docstrings do not properly link to the docstring in the source code. They simply only link to the main repo URL.

Some thoughts can be found the following link to my draft PR.

Test fail on MacOS in testing enviroment

I can replicate the error mentioned in #98 and #95 on MacOS, if I run the test suite using
] test
but not if I run the tests from the REPL. I guess that means there is something in a normal build which is not being replicated in the testing environment.

In the test environment, I've tracked this down to a NaN cropping up in the Vcycle.

Diverging pressure for rotational motion

Hi Gabriel,

I know you are busy working on another project, but I've been struggling with this issue for the past few weeks and Im starting to run out of ideas on what could solve this.

I've come across this issue while working on this kind of simulation:

Everything works fine when I try to display the vorticity, but when it comes to the pressure, the mean skyrockets:

Substracting it at every new iteration makes the pressure seem like there was no issue in the first place, but I'm sure there could be a way to not use this trick but the raw pressure instead.

I tried finding what was causing this issue and I thus used your example of TwoD_hover.jl which is also using a rotational motion. I indeed find similar results with a pressure that is not physically consistent. I then tried to figure out if it was only due to motions cumputed with rotation matrix, but it seems that sin(omega*t) displacement also cause this issue. Same for sin(kx-omega*t) if k is too small. Here you can see a motion of this kind with k=1.5 which is around the smallest value I can make it work with:

I also tried to play on the amplitude of the rotation and it indeed have an impact. When the rotation is really small the pressure gets back to normal, but what is weird is that the smallest amplitude working is way less than the amplitude of the tail in the last figure.

Since I've seen that the foil you are using in your new project is having a rotational motion, I was therefore wondering if you had ever encoutered this issue before or if you had any idea on what could be the reason behind this phenomenon.

I'll keep on trying to find a solution on my own if you don't have time to look at this in depth.

Thanks,
Bastien

Simulations of wing-like structures always have detached flow

In trying to simulate a crude wing in 2-D, I always seem to have flow that detaches from the top of the wing. This happens with low and high reynolds numbers and with low and high viscosities and very low angle of attack.

Could you help me build a simple example that does this correctly?

Here is an example:

Pluto_test.jl

Line 18: I needed to import 'Plots' as well as 'WaterLily, PlutoUI'. Would it be better to create a pull request on my fork?

Saving each time step flow field values

I don't know much about Julia. What I want to have from the simulation of a flow over a cylinder is extracting the full field of velocity and vorticity for a specific time duration. In other words, I want to save some arrays like [T,n,m], T as time steps, M and N as spatial domains. Could you please help me how can I do that?

using WaterLily
using LinearAlgebra: norm2
Re=250
n=3*2^7
m=2^8
radius=m/6
center, ν = m/2, radius/Re
body = AutoBody((x,t)->norm2(x .- center) - radius)
sim=Simulation((n+2,m+2), [1.,0.], radius; ν, body)
sim_step!(sim,10)
using MATLAB
write_matfile("Desktop\test_u.mat";sim.flow.u)
vrtx=sim.flow.σ
write_matfile("Desktop\test_vrtx.mat";vrtx)

Also, is "sim.flow.σ" vorticity of flow?

Any help would be appreciated.

Thank you!

Requires > 1.4

The code looks interesting! Just a commen treally - does it really require >1.5 version?
Current stable is 1.4.2 and I have that installed.

Taylor Vortex file

Sorry, I am a beginner with this. I downloaded everything including CUDA, but when I attempt to run this Taylor Vortex thing I get this error message. What am I doing wrong.

ERROR: AssertionError: CUDA.functional()
Stacktrace:
[1] top-level scope
@ ~/Documents/Julia-March-24/Tiina-WL-2.jl:18

Plotting gif of examples seems broken

Hello again,

when trying to plot the 2D examples (such as the TwoD_circle) with the functions provided in TwoD_plots.jl,
the sim_gif! function returns the following error:

MethodError: no method matching _as_gradient(::Array{RGBA{Float64},1}) Closest candidates are: _as_gradient(!Matched::ColorGradient) at C:\Users\Carlo\.julia\packages\Plots\8GUYs\src\utils.jl:136 _as_gradient(!Matched::Colorant) at C:\Users\Carlo\.julia\packages\Plots\8GUYs\src\utils.jl:137 in top-level scope at lily_plots.jl:58 in at lily_plots.jl:28 in #sim_gif!#58 at lily_plots.jl:29 in macro expansion at base\timing.jl:174 in macro expansion at Plots\8GUYs\src\animation.jl:181 in macro expansion at lily_plots.jl:34 in macro expansion at Plots\8GUYs\src\animation.jl:170 in frame at Plots\8GUYs\src\animation.jl:18 in frame at Plots\8GUYs\src\animation.jl:20 in png at Plots\8GUYs\src\output.jl:7 in show at Plots\8GUYs\src\output.jl:215 in _show at Plots\8GUYs\src\backends\gr.jl:1955 in gr_display at Plots\8GUYs\src\backends\gr.jl:723 in gr_display at Plots\8GUYs\src\backends\gr.jl:1550 in gr_set_gradient at Plots\8GUYs\src\backends\gr.jl:661 in gr_set_gradient at Plots\8GUYs\src\backends\gr.jl:650

I have played around with the whole thing a bit, and found that using the mu-body_plot! function also returns the same error.
Is this maybe due to some updated dependencies, since it seems to originate somewhere in the GR backend?
I am using [email protected] with [email protected]

Thanks for your help in advance!

Write VTK size/shape of field in writer

WaterLily.jl/ext/WaterLilyWriteVTKExt.jl

Line 58 in 52a66f3

k = w.count[1]; N=size(a.flow.p)

Is it possible to change this line to:
k = w.count[1]; N=size(w.output_attrib"p")

I think that should work generally and gives cleaner vti files (in case you decided to crop the fields before saving for example)

documentation/commnts are missing

What problem WaterLilly solves? Which method it uses?

Cannot run the benchmarks: `ERROR: Case arguments are arrays of different sizes.`

$ WATERLILY_ROOT=$(julia --project=. --startup-file=no -e 'using WaterLily; print(pkgdir(WaterLily))')
$ bash "${WATERLILY_ROOT}/benchmark/benchmark.sh"  -v "1.10.0" -t "12 24 36 48 60 72" -c "tgv jelly" -p "5,6,7,8"
ERROR: Case arguments are arrays of different sizes.

CC @b-fg

Modelling porous region using extended Darcy-Forchheimer-Brinkman Model

I would love to contribute in developing solvers for simulating convective heat transfer in packed bed using the extended Darcy-Brinkman-Forchheimer model.

Confusion with the Circle Tutorial and Unicode Characters

Hello,
first up thanks for this package, it looks really interesting.
When I was following the first code example, I stumbled upon some weird errors in Julia even though my code seemed literally the same. Upon using a code diff checker I noted that the "v" used in the example is not the latin "v" (U+0076) but the greek nu (U+03BD), which tripped up my code which was half typed by me and half copied from the example.
I guess the usage of the greek nu is intended to denote the correct formula symbol, but I would suggest to add a clarification in the tutorial to prevent confusion with these letters without optical differences.
Cheers, Carlo

Particle flows visualization

I have adapted the TwoBody_circle example to use a different shape. The flow seems to penetrate into the body.

Is this a problem?

client code is missing

Is there a code for the simulation from a twit?
https://twitter.com/gabrielweymouth/status/1281979466319896576

Prescribing boundary conditions

I've been playing around with the package a bit for 2D flows and know I can prescribe constant valued Dirichlet boundary conditions on the x and y boundaries using an array of length 2 (for example [1.0,0.0]) as an argument to the Simulation constructor. Is there a way to assign non constant boundary conditions and Neumann conditions? In particular I am interested in prescribing periodic conditions.

Issue when saving large number of file

This messes up the file name when we reach a large number of files (>100)

https://github.com/weymouth/WaterLily.jl/blob/888458479f62eb17282aedf822a46234b942a6d0/src/vtkWriter.jl#L39C41-L39C62

3D simulations for several shape: flow u is NaN or 0

I discovered this simulator and I am really appreciating it.

I have been playing with it and I found out that in several case, when dealing with 3D shape, like a cylinder, the sim.flow.u array is composed by only 0 and NaN.

Test it trying the donut example changing its sdf.

body = AutoBody() do xyz, t
    x, y, z = xyz - center
    # donut
    # norm2(SA[x, norm2(SA[y, z])-R]) - r 
    
    # cylinder
    d = SA[norm2(SA[y, z])-R, abs(x)-r]
    min(maximum(d), 0.0) + norm2(max.(d, 0.0))
end

Now, maybe I misunderstood something or it is just an unlucky choice of parameters, or it is a problem.

Edit:
The same problem arises with a 2D box.

LoadError: MethodError: no method matching record in ThreeD_donut although another package (ThreeD_TaylorGreenVortex) with “record” works

First package

using WaterLily
using LinearAlgebra: norm2
using Makie

function donut_sim(p=6,Re=1e3)
    # Define simulation size, velocity, viscosity
    n,U = 2^p, [1, 0, 0]
    center,R,r = [n/2,n/2,n/2], n/4, n/16
    ν = norm2(U)*R/Re

    # Apply signed distance function for a torus
    c = BDIM_coef(2n+2,n+2,n+2,3) do xyz  #
        x,y,z = xyz - center
        norm2([x,norm2([y,z])-R])-r
    end

    # Initialize Flow, Poisson and make struct
    u = zeros(2n+2,n+2,n+2,3)
    a = Flow(u,c,U,ν=ν)
    b = MultiLevelPoisson(c)
    Simulation(norm2(U),R,a,b),center
end

function flowdata(sim)
    @inside sim.flow.σ[I] = WaterLily.ω_θ(I,[1,0,0],center,sim.flow.u)*sim.L/sim.U
    @view sim.flow.σ[2:end-1,2:end-1,2:end-1]
end
function geomdata(sim)
    @inside sim.flow.σ[I] = sum(sim.flow.μ₀[I,i]+sim.flow.μ₀[I+δ(i,I),i] for i=1:3)
    @view sim.flow.σ[2:end-1,2:end-1,2:end-1]
end

function make_video!(sim::Simulation;name="file.mp4",verbose=true,t₀=0.0,Δprint=0.1,nprint=24*3)
    # plot the geometry and flow
    scene = contour(geomdata(sim),levels=[0.5])
    scene = contour!(scene,flowdata(sim),levels=[-7,7],
                     colormap=:balance,alpha=0.2,colorrange=[-7,7])
    scene_data = scene[end]

    # Plot flow evolution
    tprint = t₀+WaterLily.sim_time(sim)
    record(scene,name,1:nprint,compression=5) do i
        tprint+=Δprint
        sim_step!(sim,tprint;verbose)
        println("video ",round(Int,i*100/nprint),"% complete")
        scene_data[1] = flowdata(sim)
    end
    return scene
end

donut,center = donut_sim();
scene = make_video!(donut);

Second Package:

using WaterLily
using LinearAlgebra: norm2
using Makie

function flowdata(sim)
    @inside sim.flow.σ[I] = WaterLily.ω_mag(I,sim.flow.u)*sim.L/sim.U
    return @view sim.flow.σ[2:end-1,2:end-1,2:end-1]
end
function TGV_video(p=6,Re=1e5,Δprint=0.1,nprint=100)
    # Define vortex size, velocity, viscosity
    L = 2^p; U = 1; ν = U*L/Re

    # Taylor-Green-Vortex initial velocity field
    u = apply(L+2,L+2,L+2,3) do i,vx
        x,y,z = @. (vx-1.5)*π/L                # scaled coordinates
        i==1 && return -U*sin(x)*cos(y)*cos(z) # u_x
        i==2 && return  U*cos(x)*sin(y)*cos(z) # u_y
        return 0.                              # u_z
    end

    # Initialize simulation
    c = ones(L+2,L+2,L+2,3)  # no immersed solids
    a = Flow(u,c,zeros(3),ν=ν)
    b = MultiLevelPoisson(c)
    sim = Simulation(U,L,a,b)

    # plot the vorticity modulus
    scene = Scene(backgroundcolor = :black)
    scene = volume!(scene,flowdata(sim),colorrange=(π,4π),algorithm = :absorption)
    vol_plot = scene[end]

    # Plot flow evolution
    tprint = 0.0
    record(scene,"file.mp4",1:nprint,framerate=24,compression=5) do i
        tprint += Δprint
        sim_step!(sim,tprint)
        println("video ",round(Int,i/nprint*100),"% complete")
        vol_plot[1] = flowdata(sim)
    end
    return sim,scene
end

The first package throws the error concerning the record method, the second package does fine and has the record method as well.

ERROR: LoadError: MethodError: no method matching record(::var"#134#135"{Bool,Float64,Int64,Simulation,Contour{...}}, ::Scene, ::String, ::UnitRange{Int64}; compression=5)
Closest candidates are:
  record(::Any, ::Any, ::Any, ::Any; framerate) at /home/meck/.julia/packages/AbstractPlotting/rWoon/src/display.jl:463 got unsupported keyword argument "compression"
  record(::Any, ::Any, ::Any; framerate) at /home/meck/.julia/packages/AbstractPlotting/rWoon/src/display.jl:443 got unsupported keyword argument "compression"
Stacktrace:
 [1] kwerr(::NamedTuple{(:compression,),Tuple{Int64}}, ::Function, ::Function, ::Scene, ::String, ::UnitRange{Int64}) at ./error.jl:157
 [2] make_video!(::Simulation; name::String, verbose::Bool, t₀::Float64, Δprint::Float64, nprint::Int64) at /home/meck/Documents/Bachelor Arbeit/WaterLily/examples/ThreeD_donut.jl:42
 [3] make_video!(::Simulation) at /home/meck/Documents/Bachelor Arbeit/WaterLily/examples/ThreeD_donut.jl:35
 [4] top-level scope at /home/meck/Documents/Bachelor Arbeit/WaterLily/examples/ThreeD_donut.jl:52
in expression starting at /home/meck/Documents/Bachelor Arbeit/WaterLily/examples/ThreeD_donut.jl:52

I am confused because those functions reside in the same folder.

circle position

Hello again,

I am wondering how to define the position of the circle as:

x= a y=b, where a and be are not the same

n=32^7
m=22^7
radius=m/6
center, ν = m/2, radius/Re
body = AutoBody((x,t)->norm2(x .- center) - radius)
sim=Simulation((n+2,m+2), [1.,0.], radius; ν, body)

with this code, x=y=m/6

Thank you

Adjusting for Surface Mesh

Hi, cool package and very nice work on WaterLily.jl!

In the JuliaCon 2021 talk it is mentioned that importing a surface mesh is a fairly easy adjustment. However, I am having a go at it and could use more insight. My goal is to perform an apropos Santa's Sleigh CFD visualization. Figured out how to import the geometry, but the next step is confusing me.

Can AutoBody be used in conjunction with the surface mesh? If so then how does sdf and map come into play? Or should I consider making a new struct (something like MeshBody) when dealing with imported geometries?

Part of my confusion comes down to not fully understanding the following code from the 3D Donut example below. Maybe someone can elaborate.

body = AutoBody() do xyz,t
    x,y,z = xyz - center
    norm2([x,norm2([y,z])-R])-r
end

Repeated AutoBody additions cause StackOverflow

Respected sir,

I am facing an issue while using this code. Actually, when I am using more than a certain number (6) of members as "sdf" in a code, it gives me "StackOverflowError". Can you please help me figure out any possible solution for it?

Thank you for your time and consideration.

Yours sincerely,
Kunal Ghosh

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Jellyfish example: contour not exported by GLMakie

I think here it should be changed to GLMakie.contour.

Use KernelAbstractions for loops?

Using a combination of KernelAbstractions and broadcasting should make the code GPU compatible and multithreaded without much additional effort.

https://github.com/JuliaGPU/KernelAbstractions.jl

KernelAbstractions can also help catch a few other performance gotchas, like a failure to use fused multiply-add instructions by default:

https://www.juliabloggers.com/optimizing-julia-code-improving-the-performance-of-entropy-conservative-dg-methods-in-trixi-jl/

Fluid Structure Interaction

I've been experimenting in WaterLily with heaving airfoils recently and I'm interested in implementing leading-edge flexibility via a torsional spring model. Is fluid structure interaction feasible with WaterLily? Has any work been done on coupling body displacement/dynamics with fluid forcing? I'm newer to Julia and WaterLily any advice would be appreciated.

Thanks for this great package!

Error in test and in example

ThreeD_TaylorGreenVortex.jl

ulia> TGV_video()
1%, tU/L=0.0039, Δt=0.25
1%, tU/L=0.0117, Δt=0.5
1%, tU/L=0.0195, Δt=0.5
1%, tU/L=0.0273, Δt=0.5
1%, tU/L=0.0351, Δt=0.5
1%, tU/L=0.043, Δt=0.5
1%, tU/L=0.0508, Δt=0.5
1%, tU/L=0.0586, Δt=0.5
1%, tU/L=0.0664, Δt=0.5
1%, tU/L=0.0742, Δt=0.501
1%, tU/L=0.0821, Δt=0.501
1%, tU/L=0.0899, Δt=0.501
1%, tU/L=0.0977, Δt=0.501
ERROR: ModernGL.ContextNotAvailable("glTexSubImage3D, not available for your driver, or no valid OpenGL context available")
Stacktrace:
 [1] getprocaddress_e at C:\Users\vsarn\.julia\packages\ModernGL\rVuW2\src\ModernGL.jl:45 [inlined]
 [2] glTexSubImage3D(::UInt32, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::Int64, ::UInt32, ::UInt32, ::Array{Float32,3}) at C:\Users\vsarn\.julia\packages\ModernGL\rVuW2\src\functionloading.jl:66
 [3] texsubimage at C:\Users\vsarn\.julia\packages\GLMakie\Y3YHT\src\GLAbstraction\GLTexture.jl:374 [inlined] (repeats 2 times)
 [4] gpu_setindex!(::GLMakie.GLAbstraction.Texture{Float32,3}, ::Array{Float32,3}, ::UnitRange{Int64}, ::UnitRange{Int64}, ::Vararg{UnitRange{Int64},N} where N) at C:\Users\vsarn\.julia\packages\GLMakie\Y3YHT\src\GLAbstraction\GLTexture.jl:275
 [5] setindex!(::GLMakie.GLAbstraction.Texture{Float32,3}, ::Array{Float32,3}, ::UnitRange{Int64}, ::UnitRange{Int64}, ::UnitRange{Int64}) at C:\Users\vsarn\.julia\packages\GLMakie\Y3YHT\src\GLAbstraction\AbstractGPUArray.jl:60
 [6] update!(::GLMakie.GLAbstraction.Texture{Float32,3}, ::Array{Float32,3}) at C:\Users\vsarn\.julia\packages\GLMakie\Y3YHT\src\GLAbstraction\AbstractGPUArray.jl:80
 [7] (::GLMakie.GLAbstraction.var"#10#11"{GLMakie.GLAbstraction.Texture{Float32,3}})(::Array{Float32,3}) at C:\Users\vsarn\.julia\packages\GLMakie\Y3YHT\src\GLAbstraction\AbstractGPUArray.jl:211
 [8] #invokelatest#1 at .\essentials.jl:712 [inlined]
 [9] invokelatest at .\essentials.jl:711 [inlined]
 [10] setindex!(::Observable{Array{Float32,3}}, ::Array{Float32,3}; notify::Observables.var"#6#8") at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:132
 [11] setindex! at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:126 [inlined]
 [12] MapUpdater at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:241 [inlined]
 [13] (::Observables.OnUpdate{Observables.MapUpdater{AbstractPlotting.var"#183#185"{Int64},Array{Float32,3}},Tuple{Observable{Tuple{IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},Array{Float32,3}}}}})(::Tuple{IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},Array{Float32,3}}) at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:218
 [14] setindex!(::Observable{Tuple{IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},Array{Float32,3}}}, ::Tuple{IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},Array{Float32,3}}; notify::Observables.var"#6#8") at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:130
 [15] setindex!(::Observable{Tuple{IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},Array{Float32,3}}}, ::Tuple{IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},Array{Float32,3}}) at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:126
 [16] (::AbstractPlotting.var"#201#203"{DataType,Observable{Tuple{IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},Array{Float32,3}}}})(::Tuple{}, ::Tuple{Array{Float64,3}}) at C:\Users\vsarn\.julia\packages\AbstractPlotting\tHSnG\src\interfaces.jl:615
 [17] (::Observables.OnUpdate{AbstractPlotting.var"#201#203"{DataType,Observable{Tuple{IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},IntervalSets.Interval{:closed,:closed,Float32},Array{Float32,3}}}},Tuple{Observable{Tuple{}},Observable{Tuple{Array{Float64,3}}}}})(::Tuple{Array{Float64,3}}) at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:218
 [18] setindex!(::Observable{Tuple{Array{Float64,3}}}, ::Tuple{Array{Float64,3}}; notify::Observables.var"#6#8") at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:130
 [19] setindex! at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:126 [inlined]
 [20] MapUpdater at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:241 [inlined]
 [21] (::Observables.OnUpdate{Observables.MapUpdater{typeof(tuple),Tuple{Array{Float64,3}}},Tuple{Observable{Array{Float64,3}}}})(::Array{Float64,3}) at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:218
 [22] setindex!(::Observable{Array{Float64,3}}, ::Array{Float64,3}; notify::Observables.var"#6#8") at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:130
 [23] setindex!(::Observable{Array{Float64,3}}, ::Array{Float64,3}) at C:\Users\vsarn\.julia\packages\Observables\0wrF6\src\Observables.jl:126
 [24] setindex!(::Volume{...}, ::Array{Float64,3}, ::Int64) at C:\Users\vsarn\.julia\packages\AbstractPlotting\tHSnG\src\dictlike.jl:190
 [25] (::var"#11#15"{Int64,Int64,Flow{4,3},MultiLevelPoisson{4,3},Volume{...},Float64,Int64})(::Int64) at C:\Users\vsarn\Documents\FrozenCat\Julia\Test\water.jl:42
 [26] record(::var"#11#15"{Int64,Int64,Flow{4,3},MultiLevelPoisson{4,3},Volume{...},Float64,Int64}, ::Scene, ::String, ::UnitRange{Int64}; framerate::Int64, compression::Int64, sleep::Bool) at C:\Users\vsarn\.julia\packages\AbstractPlotting\tHSnG\src\display.jl:537
 [27] TGV_video(::Int64, ::Float64) at C:\Users\vsarn\Documents\FrozenCat\Julia\Test\water.jl:30
 [28] TGV_video() at C:\Users\vsarn\Documents\FrozenCat\Julia\Test\water.jl:7
 [29] top-level scope at none:0

Test:

julia> Pkg.test("WaterLily")
    Testing WaterLily
Status `C:\Users\vsarn\AppData\Local\Temp\jl_ycgAph\Manifest.toml`
  [ed894a53] WaterLily v0.1.0 #master (https://github.com/weymouth/WaterLily.git)
  [2a0f44e3] Base64
  [8ba89e20] Distributed
  [b77e0a4c] InteractiveUtils
  [8f399da3] Libdl
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging
  [d6f4376e] Markdown
  [9a3f8284] Random
  [9e88b42a] Serialization
  [6462fe0b] Sockets
  [8dfed614] Test
Test Summary: | Pass  Total
util.jl       |    2      2
Test Summary: | Pass  Total
Poisson.jl    |    2      2
Test Summary:        | Pass  Total
MultiLevelPoisson.jl |    2      2
Test Summary: | Pass  Total
Body.jl       |    5      5
Flow.jl: Test Failed at C:\Users\vsarn\.julia\packages\WaterLily\4wkpL\test\runtests.jl:66
  Expression: L₂(a.u[:, :, 2] .- U[2]) < 1.0e-6
   Evaluated: 1.4910472215055462e-6 < 1.0e-6
Stacktrace:
 [1] top-level scope at C:\Users\vsarn\.julia\packages\WaterLily\4wkpL\test\runtests.jl:66
 [2] top-level scope at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Test\src\Test.jl:1113
 [3] top-level scope at C:\Users\vsarn\.julia\packages\WaterLily\4wkpL\test\runtests.jl:59
Test Summary: | Pass  Fail  Total
Flow.jl       |    1     1      2
ERROR: LoadError: Some tests did not pass: 1 passed, 1 failed, 0 errored, 0 broken.
in expression starting at C:\Users\vsarn\.julia\packages\WaterLily\4wkpL\test\runtests.jl:57
ERROR: Package WaterLily errored during testing
Stacktrace:
 [1] pkgerror(::String, ::Vararg{String,N} where N) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Pkg\src\Types.jl:53
 [2] test(::Pkg.Types.Context, ::Array{Pkg.Types.PackageSpec,1}; coverage::Bool, julia_args::Cmd, test_args::Cmd, test_fn::Nothing) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Pkg\src\Operations.jl:1503
 [3] test(::Pkg.Types.Context, ::Array{Pkg.Types.PackageSpec,1}; coverage::Bool, test_fn::Nothing, julia_args::Cmd, test_args::Cmd, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Pkg\src\API.jl:316
 [4] test(::Pkg.Types.Context, ::Array{Pkg.Types.PackageSpec,1}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Pkg\src\API.jl:303
 [5] #test#68 at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Pkg\src\API.jl:297 [inlined]
 [6] test at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Pkg\src\API.jl:297 [inlined]
 [7] #test#67 at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Pkg\src\API.jl:296 [inlined]
 [8] test at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Pkg\src\API.jl:296 [inlined]
 [9] test(::String; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Pkg\src\API.jl:295
 [10] test(::String) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Pkg\src\API.jl:295
 [11] top-level scope at none:0

julia>

README.md needs update? MethodError: no method matching sim_step!(::Simulation; t_end::Int64)

An attempt to run the exact code from the example 'Flow over a circle' in README.md produces the following error:

julia> sim_step!(circ,t_end=10)
ERROR: MethodError: no method matching sim_step!(::Simulation; t_end::Int64)

Closest candidates are:
  sim_step!(::Simulation, ::Any; verbose, remeasure) got unsupported keyword argument "t_end"
   @ WaterLily ~/.julia/packages/WaterLily/rv3IZ/src/WaterLily.jl:83

It seems README.md needs updating

Vtk writer, writing always to vtk_data folder

Currently the vtk writer always writes to the vtk_data folder. In case multiple simulations are done after each other, this will fill up the same folder with all the results. In the filename you can still specify which case it is, but it would be cleaner to have the option to specificy the folder where you would like the vti & pvd files to be stored

https://github.com/weymouth/WaterLily.jl/blob/888458479f62eb17282aedf822a46234b942a6d0/src/vtkWriter.jl#L16C66-L16C80

weymouth / waterlily.jl Goto Github PK

waterlily.jl's People

Contributors

Stargazers

Watchers

Forkers

waterlily.jl's Issues

The code is as follows:

Recommend Projects

Recommend Topics

Recommend Org