Pressure solves on `GPU` are not ready for `Flat`

[francispoulin]

I have put together an example for my own benefit (not to be merged into master unless people want it to be) that looks at the simulation of inertial instability in 2D. See here for the code.

I believe everything is working on a CPU but I am having two difficulties with GPUs.

1. In line 39, where I define the background buoyancy of the jet, I can't use coriolis.f because CUDA seems to need a global variable and this doesn't cut it. I am presently using f instead and this works, but should this work? ERROR: LoadError: InvalidIRError: compiling kernel gpu__compute!
2. When I make a simulation from a gpu calculation, I get that the buoyancy perturbation has perturbations at the top and bottom, which are not physical. Which is to say that it has different boundary conditions. I didn't actually specify the boundary conditions differently but just wanted a solid top and bottom. If I specify the boundary conditions explicitly, should this fix the problem? Maybe that's a fix but is this expected behaviour?

https://user-images.githubusercontent.com/8239041/113899192-8d851900-979a-11eb-97a7-b7f8085864b7.mp4

[glwagner]

In line 39, where I define the background buoyancy of the jet, I can't use coriolis.f because CUDA seems to need a global variable and this doesn't cut it.

This isn't a GPU issue I don't think because initialization is done on the CPU. I don't see a coriolis object defined in your script. You may need to write model.coriolis.f.

[glwagner]

When I make a simulation from a gpu calculation, I get that the buoyancy perturbation has perturbations at the top and bottom, which are not physical. Which is to say that it has different boundary conditions. I didn't actually specify the boundary conditions differently but just wanted a solid top and bottom. If I specify the boundary conditions explicitly, should this fix the problem? Maybe that's a fix but is this expected behaviour?

Can you show both CPU and GPU results?

[francispoulin]

Thanks @glwagner .

First, what I tried and does not work is the following

coriolis = FPlane(f=f)
b̄(x, y, z) = U_max * tanh(y/L_jet)   * exp( - (z + D)^2/H_jet^2 ) * 2 * coriolis.f * L_jet / H_jet^2 * (z + D)

Should this work?

Second, the above is the gpu result and I am now including the cpu.

https://user-images.githubusercontent.com/8239041/113900529-e903d680-979b-11eb-9ef4-9d53d61bc5e8.mp4

[glwagner]

Should this work?

set! on the GPU does this:

1. Create a CPU field
2. call set!(cpu_field, f) with the function f
3. Copy CPU data to GPU

So initialization is done on the CPU, not the GPU. If set! works on the CPU then it should work on the GPU. The only place where something might go wrong is copying the data over, somehow.

You can see this here:

https://github.com/CliMA/Oceananigans.jl/blob/5c4f493ee9795613ca2cf83050ad7fe9687e2ec6/src/Fields/set!.jl#L94-L102

[ali-ramadhan]

@francispoulin I think it's always helpful to post the full error and stacktrace. From your original post

ERROR: LoadError: InvalidIRError: compiling kernel gpu__compute!

suggests that the error could be happening since b̄ is being used in a ComputedField's compute!: https://github.com/CliMA/Oceananigans.jl/blob/9eca5780658bb8f5c0debd34146a0ad5cb73c872/examples/inertially_unstable_jet.jl#L51

[ali-ramadhan]

On the CPU vs. GPU issue. Your GPU movie with the artifacts kinda reminds of https://github.com/CliMA/Oceananigans.jl/issues/1170 which had artifacts at the top-most grid point due to a bug in the pressure solver.

Is there a reason to think that the pressure solver is wrong for the case when there's a background buoyancy field? Is it forgetting to add a tendency to the source term or something maybe? Actually it seems to only be a GPU issue so can't be the source term.

[glwagner]

Oh nice, the first use I've seen that mixes a function with a concrete Field in an AbstractOperation! Sweet! I'm not sure it would work but you can try const coriolis = Coriolis(f=f). Alternatively you can build the FunctionField manually with something like

b̄(x, y, z, f=coriolis.f) = # definition
b̄_field = FunctionField{Center, Center, Center}(b̄, grid, parameters=coriolis.f)

b̃ = ComputedField(b - b̄_field) 

But there's not too much of a downside in writing const f = # whatever right?

[francispoulin]

Thansk for the quick feedback.

1. I should say this is not as much of a concern as I found a work around. In my definition of b\tilde, I changed f to model.coriolis.f and received an error. Below is the beginning and it's huge so can't copy the whole thing. I will stick to my simple solution for the moment but want to point this out, in case there was a concern.

 include("inertially_unstable_jet.jl")
ERROR: LoadError: InvalidIRError: compiling kernel gpu__compute!(Cassette.Context{nametype(CUDACtx),KernelAbstractions.CompilerMetadata{KernelAbstractions.NDIteration.StaticSize{(1, 64, 64)},KernelAbstractions.NDIteration.DynamicCheck,Nothing,Nothing,KernelAbstractions.NDIteration.NDRange{3,KernelAbstractions.NDIteration.StaticSize{(1, 1, 64)},KernelAbstractions.NDIteration.StaticSize{(1, 64, 1)},Nothing,Nothing}},Nothing,KernelAbstractions.var"##PassType#253",Nothing,Cassette.DisableHooks}, typeof(Oceananigans.Fields.gpu__compute!), OffsetArrays.OffsetArray{Float64,3,CuDeviceArray{Float64,3,1}}, Oceananigans.AbstractOperations.BinaryOperation{Center,Center,Center,typeof(-),OffsetArrays.OffsetArray{Float64,3,CuDeviceArray{Float64,3,1}},Oceananigans.Fields.FunctionField{Center,Center,Center,Nothing,Nothing,typeof(b̄),RegularRectilinearGrid{Float64,Flat,Bounded,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}},typeof(identity),typeof(identity),typeof(identity),RegularRectilinearGrid{Float64,Flat,Bounded,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}) resulted in invalid LLVM IR
Reason: unsupported dynamic function invocation (call to overdub)
Stacktrace:
 [1] b̄ at /home/fpoulin/software/Oceananigans.jl/examples/inertially_unstable_jet.jl:39
 [2] call_func at /home/fpoulin/software/Oceananigans.jl/src/Fields/function_field.jl:61
 [3] getindex at /home/fpoulin/software/Oceananigans.jl/src/Fields/function_field.jl:63
 [4] identity at /home/fpoulin/software/Oceananigans.jl/src/Operators/interpolation_utils.jl:6
 [5] - at /home/fpoulin/software/Oceananigans.jl/src/AbstractOperations/binary_operations.jl:59
 [6] identity at /home/fpoulin/software/Oceananigans.jl/src/Operators/interpolation_utils.jl:11
 [7] getindex at /home/fpoulin/software/Oceananigans.jl/src/AbstractOperations/binary_operations.jl:34
 [8] macro expansion at /home/fpoulin/software/Oceananigans.jl/src/Fields/computed_field.jl:114
 [9] gpu__compute! at /home/fpoulin/.julia/packages/KernelAbstractions/mKsXc/src/macros.jl:80
 [10] overdub at /home/fpoulin/.julia/packages/Cassette/Wjztv/src/overdub.jl:0
Reason: unsupported dynamic function invocation (call to overdub)
Stacktrace:
 [1] - at /home/fpoulin/software/Oceananigans.jl/src/AbstractOperations/binary_operations.jl:59
 [2] identity at /home/fpoulin/software/Oceananigans.jl/src/Operators/interpolation_utils.jl:11
 [3] getindex at /home/fpoulin/software/Oceananigans.jl/src/AbstractOperations/binary_operations.jl:34
 [4] macro expansion at /home/fpoulin/software/Oceananigans.jl/src/Fields/computed_field.jl:114
 [5] gpu__compute! at /home/fpoulin/.julia/packages/KernelAbstractions/mKsXc/src/macros.jl:80
 [6] overdub at /home/fpoulin/.julia/packages/Cassette/Wjztv/src/overdub.jl:0
Reason: unsupported dynamic function invocation (call to overdub)
Stacktrace:
 [1] macro expansion at /home/fpoulin/software/Oceananigans.jl/src/Fields/computed_field.jl:114
 [2] gpu__compute! at /home/fpoulin/.julia/packages/KernelAbstractions/mKsXc/src/macros.jl:80
 [3] overdub at /home/fpoulin/.julia/packages/Cassette/Wjztv/src/overdub.jl:0
Reason: unsupported dynamic function invocation (call to getproperty)
Stacktrace:
 [1] call at /home/fpoulin/.julia/packages/Cassette/Wjztv/src/context.jl:456
 [2] fallback at /home/fpoulin/.julia/packages/Cassette/Wjztv/src/context.jl:454
 [3] overdub at /home/fpoulin/.julia/packages/Cassette/Wjztv/src/context.jl:279
 [4] b̄ at /home/fpoulin/software/Oceananigans.jl/examples/inertially_unstable_jet.jl:39
 [5] call_func at /home/fpoulin/software/Oceananigans.jl/src/Fields/function_field.jl:61
 [6] getindex at /home/fpoulin/software/Oceananigans.jl/src/Fields/function_field.jl:63
 [7] identity at /home/fpoulin/software/Oceananigans.jl/src/Operators/interpolation_utils.jl:6
 [8] - at /home/fpoulin/software/Oceananigans.jl/src/AbstractOperations/binary_operations.jl:59
 [9] identity at /home/fpoulin/software/Oceananigans.jl/src/Operators/interpolation_utils.jl:11
 [10] getindex at /home/fpoulin/software/Oceananigans.jl/src/AbstractOperations/binary_operations.jl:34
 [11] macro expansion at /home/fpoulin/software/Oceananigans.jl/src/Fields/computed_field.jl:114
 [12] gpu__compute! at /home/fpoulin/.julia/packages/KernelAbstractions/mKsXc/src/macros.jl:80
 [13] overdub at /home/fpoulin/.julia/packages/Cassette/Wjztv/src/overdub.jl:0
Stacktrace:
 [1] check_ir(::GPUCompiler.CompilerJob{GPUCompiler.PTXCompilerTarget,CUDA.CUDACompilerParams}, ::LLVM.Module) at /home/fpoulin/.julia/packages/GPUCompiler/uTpNx/src/validation.jl:123
 [2] macro expansion at /home/fpoulin/.julia/packages/GPUCompiler/uTpNx/src/driver.jl:239 [inlined]
 [3] macro expansion at /home/fpoulin/.julia/packages/TimerOutputs/4QAIk/src/TimerOutput.jl:206 [inlined]
 [4] codegen(::Symbol, ::GPUCompiler.CompilerJob; libraries::Bool, deferred_codegen::Bool, optimize::Bool, strip::Bool, validate::Bool, only_entry::Bool) at /home/fpoulin/.julia/packages/GPUCompiler/uTpNx/src/driver.jl:237
 [5] compile(::Symbol, ::GPUCompiler.CompilerJob; libraries::Bool, deferred_codegen::Bool, optimize::Bool, strip::Bool, validate:
...

[francispoulin]

2. Is more of a concern as there seems to be some anomalies developing near the boundary, that should not be there.

I should also point out that the values at the boundary are increasing very rapidly, causing a numerical instability.

I tried doing a purely 1D profile simulation in the z direction and that was stable for many days. This suggests to me that it is a 2D problem.

I am at a loss as to what to do but wonder if I should be explicitly setting the boundary conditions, as they seem differenet in the two cases.

Any suggestions of what I can try to fix this?

[ali-ramadhan]

I'm not sure where to look right now but I think it's pretty worrying that the CPU and GPU simulations are producing different results. They should produce the exact same results (up to machine precision), or at least this is what Oceananigans.jl claims to do.

[francispoulin]

I did a coarse run with 32x32 and it still develops this numerial instability.

Below is a plot of the buoyancy at a constant value of y, to show the vertical profile of buoyancy. We see that we do not have Neumann (Derivative) boundary condtiions in b, which I think we should. No?

[ali-ramadhan]

Hmmm, the default boundary condition for tracers on walls is no-flux which I think does translate to zero Neumann boundary condition.

But Oceananigans stores the value of b at (Center, Center, Center) so it might not look exactly like a Neumann boundary condition since the gradient goes to zero at the wall which is located at a Face...

The plot you posted look very far from a Neumann boundary condition though.

Hmmm, I wonder if it has something to do with ComputedField boundary conditions. But the weird thing is that it's fine on the CPU...

[francispoulin]

Minor update:

I have a reduced case that is 4x4 that I'm playing with to try and determine why the results from CPUs differ from GPUs.

After one time step, using QuasiAdamsBashforth2, I get that the two sets of values in the first column are identical.

CPU
10-element OffsetArray(::Array{Float64,1}, -2:7) with eltype Float64 with indices -2:7:
 0.006762989865623905
 0.006762989865623905
 2.673025080318251e-8
 2.673025080318251e-8
 0.006762989865623905
 0.006762989865623905
 2.673025080318251e-8
 2.673025080318251e-8
 0.006762989865623905
 0.006762989865623905

GPU
10-element OffsetArray(::CuArray{Float64,1}, -2:7) with eltype Float64 with indices -2:7:
 0.006762989865623905
 0.006762989865623905
 2.673025080318251e-8
 2.673025080318251e-8
 0.006762989865623905
 0.006762989865623905
 2.673025080318251e-8
 2.673025080318251e-8
 0.006762989865623905
 0.006762989865623905

However, when I use RungeKutta3 I get that the results are different.

CPU
10-element OffsetArray(::Array{Float64,1}, -2:7) with eltype Float64 with indices -2:7:
 0.00676299502124483
 0.00676299502124483
 1.696248222222971e-8
 1.696248222222971e-8
 0.00676299502124483
 0.00676299502124483
 1.6962482222229713e-8
 1.6962482222229713e-8
 0.00676299502124483
 0.00676299502124483

GPU
10-element OffsetArray(::CuArray{Float64,1}, -2:7) with eltype Float64 with indices -2:7:
 0.006762989865623905
 0.006762989865623905
 2.673025080318251e-8
 2.673025080318251e-8
 0.006762989865623905
 0.006762989865623905
 2.673025080318251e-8
 2.673025080318251e-8
 0.006762989865623905
 0.006762989865623905

I know the differences are small, 8th decimal place or so, but on such a coarse grid should we be expecting the same answers?

I observed that QAB2 then starts to differ at the next time step. I suppose the differences in the RK3 method appear sooner because it is a multistep method, so three steps in one?

I know that the pressure solve is what differs significantly between the CPU and GPU approaches but are these reasonable errors or do these seem too big?

I can say that the errors at the boundary happen with both time stepping methods. So I cannot run this example on a GPU because of this numerical instability that develops. I know that many other examples run successfully, but not sure what I am doing differently compared other codes.

[glwagner]

The second time-step is when QuasiAdamsBashforth2 uses tendency information saved from the previous time-step. On the first time-step it's equivalent to a forward Euler scheme.

I know the differences are small, 8th decimal place or so, but on such a coarse grid should we be expecting the same answers?

Hardware differences can produce results that differ by something like sqrt(eps(Float64)) I think. This is the tolerance we use for regression tests (but even those sometimes fail to pass on the GPU for some reason).

Since the background buoyancy field has a linear gradient and your domain is bounded, can you try running the same simulation but explicitly resolving the background buoyancy field? (This could suggest that the issue has to do with BackgroundField...)

[francispoulin]

Very good point. The background buoyancy does not satisfy the boundary conditions that are imposed. I will check it out later.

If this is the case, then it is still mysterious why GPUs have a harder time with this than CPUs, but we will see and learn. Thanks for the suggestion.

Francis

---

From: Gregory L. Wagner @.> Sent: Wednesday, April 7, 2021 5:44:23 PM To: CliMA/Oceananigans.jl @.> Cc: Francis Poulin @.>; Mention @.> Subject: Re: [CliMA/Oceananigans.jl] strange results on gpus (#1554)

The second time-step is when QuasiAdamsBashforth2 uses tendency information saved from the previous time-step. On the first time-step it's equivalent to a forward Euler scheme.

I know the differences are small, 8th decimal place or so, but on such a coarse grid should we be expecting the same answers?

Hardware differences can produce results that differ by something like sqrt(eps(Float64)) I think. This is the tolerance we use for regression tests (but even those sometimes fail to pass on the GPU for some reason).

Since the background buoyancy field has a linear gradient and your domain is bounded, can you try running the same simulation but explicitly resolving the background buoyancy field? (This could suggest that the issue has to do with BackgroundField...)

- You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/CliMA/Oceananigans.jl/issues/1554#issuecomment-815284094, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AB63PQICFGOXKMU4HFN3NOLTHTG3PANCNFSM42RDHB2Q.

[francispoulin]

I thought I would also try freezing the velocity field, if possible. If this still weirdness at the boundary still happens when the velocity doesn't evolve, then it has nothing to do with the pressure solve.

Could someone point me to an example that does this?

---

From: Francis Poulin @.> Sent: Wednesday, April 7, 2021 7:43:11 PM To: CliMA/Oceananigans.jl @.>; CliMA/Oceananigans.jl @.> Cc: Mention @.> Subject: Re: [CliMA/Oceananigans.jl] strange results on gpus (#1554)

Very good point. The background buoyancy does not satisfy the boundary conditions that are imposed. I will check it out later.

If this is the case, then it is still mysterious why GPUs have a harder time with this than CPUs, but we will see and learn. Thanks for the suggestion.

Francis

---

From: Gregory L. Wagner @.> Sent: Wednesday, April 7, 2021 5:44:23 PM To: CliMA/Oceananigans.jl @.> Cc: Francis Poulin @.>; Mention @.> Subject: Re: [CliMA/Oceananigans.jl] strange results on gpus (#1554)

The second time-step is when QuasiAdamsBashforth2 uses tendency information saved from the previous time-step. On the first time-step it's equivalent to a forward Euler scheme.

I know the differences are small, 8th decimal place or so, but on such a coarse grid should we be expecting the same answers?

Hardware differences can produce results that differ by something like sqrt(eps(Float64)) I think. This is the tolerance we use for regression tests (but even those sometimes fail to pass on the GPU for some reason).

Since the background buoyancy field has a linear gradient and your domain is bounded, can you try running the same simulation but explicitly resolving the background buoyancy field? (This could suggest that the issue has to do with BackgroundField...)

- You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/CliMA/Oceananigans.jl/issues/1554#issuecomment-815284094, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AB63PQICFGOXKMU4HFN3NOLTHTG3PANCNFSM42RDHB2Q.

[francispoulin]

I am happy to report that following @glwagner 's suggestion, I set the buoyancy frequency to zero and that seemed to solve the problem. Nothing formed at the boundary and therefore no numerical instability. This is good evidence to suggest that it's because the background state is not being consistent with the top and bottom boundary conditions, which require a Neuman (Derivative) boundary condtiion. A big step forward!

A few questions come to mind.

• Where is it exactly that GPUs have difficulties dealing with this?
• Why can CPUs do this without any problem? Or maybe where is the difference between the two methods?

A simple solution to my problem is to change the boundary conditions to be consistent with the background state and that should fix things up nicely, i.e. set a non-zero Neumann (Derivative) Boundary Condition at the top and bottom.

I will give that a try tomorrow, but I am still interested in the differences between what the CPU and GPU are doing.

[glwagner]

I think what you have discovered is that simulations with BackgroundField can behave differently on the GPU than on the CPU for some reason.

What do you mean when you say the BackgroundField is not consistent with the boundary condition? The boundary condition is applied only to the prognostic variable (the deviation) only. Indeed, one of the most important applications of BackgroundField is to run simulations in which the deviations are periodic even though the BackgroundField is not. The Eady turbulence problem is an example of one such problem. In that case the background velocity and buoyancy fields are not periodic in the y-direction; but because they only vary linearly, deviations from that background state can be idealized as periodic.

There may indeed be some issue regarding BackgroundField and Bounded directions, but I believe this is a bug and not expected. Certainly there is a bug if different results occur on the CPU versus the GPU. Hopefully we can get to the bottom of it.

[francispoulin]

Ah, thanks for clarifying as clearly I was mistaken. I thought I needed to impose the boundary conditions on the total field, background + deviation. If that's not the case then the boundary conditions are set up correctly.

I believe the next thing to try is the same set up but without BackgroundField and imposing the Neumann boundary conditions associated with the buoyancy frequency. If that behaves correctly, then the problem would appear to be in BackgroundField, as you suspect.

Update: I made the following changes,

• Removed the background field from model,
• Added N^2 * (z + D) into the initial conditiosn of the buoyancy,
• Imposed buoyancy_gradient_bc at the top and bottom of the domain.

Unfortunately, the same problems occur at the boundary.

This does not seem to be imposing the correct boundary conditions on a tracer field at the top and bottom.

Periodic condition: I did try setting the vertical direction to periodic and unfortunately that gave an error, copied below.

Is this a seperate problem or do people think it's related?

ERROR: LoadError: ArgumentError: batching dims must be sequential
Stacktrace:
 [1] create_plan(::CUDA.CUFFT.cufftType_t, ::Tuple{Int64,Int64,Int64}, ::Array{Int64,1}) at /home/fpoulin/.julia/packages/CUDA/wTQsK/lib/cufft/fft.jl:140
 [2] plan_fft! at /home/fpoulin/.julia/packages/CUDA/wTQsK/lib/cufft/fft.jl:256 [inlined]
 [3] plan_forward_transform at /home/fpoulin/software/Oceananigans.jl/src/Solvers/plan_transforms.jl:42 [inlined]
 [4] plan_transforms(::GPU, ::RegularRectilinearGrid{Float64,Flat,Bounded,Periodic,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::CUDA.CuArray{Complex{Float64},3}, ::UInt32) at /home/fpoulin/software/Oceananigans.jl/src/Solvers/plan_transforms.jl:106
 [5] Oceananigans.Solvers.FFTBasedPoissonSolver(::GPU, ::RegularRectilinearGrid{Float64,Flat,Bounded,Periodic,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::UInt32) at /home/fpoulin/software/Oceananigans.jl/src/Solvers/fft_based_poisson_solver.jl:25
 [6] FFTBasedPoissonSolver at /home/fpoulin/software/Oceananigans.jl/src/Solvers/fft_based_poisson_solver.jl:11 [inlined]
 [7] PressureSolver(::GPU, ::RegularRectilinearGrid{Float64,Flat,Bounded,Periodic,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}) at /home/fpoulin/software/Oceananigans.jl/src/Solvers/Solvers.jl:44
 [8] IncompressibleModel(; grid::RegularRectilinearGrid{Float64,Flat,Bounded,Periodic,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, architecture::GPU, float_type::Type{T} where T, clock::Clock{Float64}, advection::WENO5, buoyancy::BuoyancyTracer, coriolis::FPlane{Float64}, stokes_drift::Nothing, forcing::NamedTuple{(),Tuple{}}, closure::AnisotropicDiffusivity{Float64,Float64,Float64,Float64,Float64,Float64}, boundary_conditions::NamedTuple{(),Tuple{}}, tracers::Symbol, timestepper::Symbol, background_fields::NamedTuple{(:b,),Tuple{BackgroundField{typeof(B_func),Float64}}}, particles::Nothing, velocities::Nothing, pressures::Nothing, diffusivities::Nothing, pressure_solver::Nothing, immersed_boundary::Nothing) at /home/fpoulin/software/Oceananigans.jl/src/Models/IncompressibleModels/incompressible_model.jl:145
 [9] top-level scope at /home/fpoulin/software/Oceananigans.jl/examples/inertially_unstable_jet.jl:31
 [10] include(::String) at ./client.jl:457
 [11] top-level scope at REPL[1]:1
in expression starting at /home/fpoulin/software/Oceananigans.jl/examples/inertially_unstable_jet.jl:31

[francispoulin]

I changed the topology in the x direction from Flat to Periodic and obtained the following two simulations on GPUs and CPUs. The good news is that they are both stable and are very similar. I wonder if the problem might be because of the Flat topology?

https://user-images.githubusercontent.com/8239041/114034660-625b0200-984c-11eb-8771-ec0a1945dd89.mp4

https://user-images.githubusercontent.com/8239041/114035326-fdec7280-984c-11eb-993c-e5ae2fb08d45.mp4

[francispoulin]

The problem was that the pressure solvers are not tested for Flat. It seems that the CPU case does the right thing but the GPU case does not. I will create an issue saying this needs to be addressed. In the mean time, @ali-ramadhan and I have added a warning to tell the user that the pressure solvers are not tested for flat. #1556.

[glwagner]

We definitely need this issue open until the problem is solved. Adding a warning doesn't solve the problem!

What happens when you use Periodic with 1 grid point rather than Flat?

[francispoulin]

The simulations above were with Periodic and only 1 grid point. That worked fine.

[glwagner]

Ok next question: what happens when we use Bounded with 1 grid point (I think the Poisson solver treats Flat as it treats Bounded, but maybe it should treat Flat as Periodic)

[glwagner]

Which direction was Flat? It looks like the code linked to in this PR has changed (one reason why its better to write code directly in the PR rather than providing links that can change...)

[francispoulin]

Good to know about the links and sorry that things have changed. Flat was in the x direction and I'm going to confirm previous results now and try bounded. Will share results as soon as I have them.

[ali-ramadhan]

Ideally Flat should be treated properly but I agree it's better to treat it was Periodic with 1 grid point (FFT shouldn't do anything) rather than as Bounded.

[francispoulin]

@glwagner : unfortunately, my desktop is acting up and I can't get very far in my simulations on a GPU without running out of memory, and it's not far at all. I don't understand this at all as it ran a couple of hours ago, perfectly fine.

However, as for the bounded case, I don't think it is a reasonable test case as in this problem the jet is in the x direction. If we change it from Periodic to Bounded, then u=0 and it will give rise to a very different scenario since I'm looking at a jet u(y,z). Do you see my concern, because of the no-normal flow boundary conditions?

When my computer seems better behaved I will happily run the example again but at the moment, sadly, I can't.

[francispoulin]

Finally, here is an animation with Bounded in the x direction. The good news is that it's numerically stable. However, I have been plotting the pertubations and the fields look very different because the solid walls force the total velocity in x to go to zero.

https://user-images.githubusercontent.com/8239041/114077635-04dbab00-9876-11eb-91d2-0e2981cfd883.mp4

@glwagner : how easy/difficult will it be to test the Poisson solvers for Flat?

[glwagner]

It's not too difficult.

[glwagner]

Closed by #1560