Adds a pertubation advection open boundary matching scheme #3977

jagoosw · 2024-12-06T19:01:09Z

This PR adds a PertubationAdvection open boundary condition which assumes that there is some mean flow (i.e. prescribed in a channel or from a coarser resolution model) which is homogeneous on the boundary but can vary in time.

On the boundary, we approximate (for an x-boundary) to:

$$\frac{\partial u}{\partial t} \approx -U\frac{\partial u\prime}{\partial x},$$

where $u=u\prime+U$
I have chosen to take an explicit/backwards Euler step:

$$\frac{u\prime(t^{n+1}) - u\prime(t^n)}{\Delta t} + \frac{U\prime(t^{n+1}) - U\prime(t^n)}{\Delta t} \approx -U\frac{u\prime(x_b, t^{n+1}) - u\prime(x_{b-1}, t^{n+1})}{\Delta x},$$

because then we don't have to store any other information (i.e. if we did a forward step we would need $u\prime(x_{b-1}, t^n)$ ). This results in:

$$u(t^{n+1}) = \frac{u^n+\tilde{U}u(x_{b-1}, t^{n+1})}{1+\tilde{U}},$$

where $\tilde{U}=\min\left(1, U\Delta t / \Delta x\right)$. I have also added restoring to $U$ (i.e. damping $u\prime\to0$) in some timescale $\tau$ which can be different for inflow and outflow to allow the scheme to work when there is inflow, which results in:

$$u(t^{n+1}) = \frac{u^n+\tilde{U}u(x_{b-1}, t^{n+1}) + U\tilde{\tau}}{1+\tilde{U}+\tilde{\tau}},$$

where $\tilde{\tau}=\Delta t / \tau$.

This has the limitation that if we want to extend to have wall tangent mean velocity ($\partial_t u \approx -(\vec{U}\cdot\nabla) u$) then we either have to solve the whole boundary at once since every point is going to depend on its boundary neighbours, or we need to take a forward Euler step in which case we need the first interior point at the previous time step.

I think we probably will need to address this for real nesting cases at some point.

u.mp4

For now, it seems to be working okay (left plot vorticity, right plot x-velocity):

(old plot)

https://github.com/user-attachments/assets/21ae3eb2-5d3b-4c08-86fb-bb9c222e8e6a

and the drag coefficient is ~1.4 at Re = 100 and I'm running a case at Re=1000 now.

I have also checked the Strouhal number which matches exactly to the expected values (~0.17).

jagoosw · 2024-12-06T19:03:32Z

@tomchor, do you think something like this:

Oceananigans.jl/validation/open_boundaries/cylinder.jl

Lines 18 to 76 in 32cd354

    
           """ 
        
               drag(modell; bounding_box = (-1, 3, -2, 2), ν = 1e-3) 
        
           Returns the drag within the `bounding_box` computed by: 
        
           ```math 
        
           \\frac{\\partial \\vec{u}}{\\partial t} + (\\vec{u}\\cdot\\nabla)\\vec{u}=-\\nabla P + \\nabla\\cdot\\vec{\\tau} + \\vec{F},\\newline 
        
           \\vec{F}_T=\\int_\\Omega\\vec{F}dV = \\int_\\Omega\\left(\\frac{\\partial \\vec{u}}{\\partial t} + (\\vec{u}\\cdot\\nabla)\\vec{u}+\\nabla P - \\nabla\\cdot\\vec{\\tau}\\right)dV,\\newline 
        
           \\vec{F}_T=\\int_\\Omega\\left(\\frac{\\partial \\vec{u}}{\\partial t}\\right)dV + \\oint_{\\partial\\Omega}\\left(\\vec{u}(\\vec{u}\\cdot\\hat{n}) + P\\hat{n} - \\vec{\\tau}\\cdot\\hat{n}\\right)dS,\\newline 
        
           \\F_u=\\int_\\Omega\\left(\\frac{\\partial u}{\\partial t}\\right)dV + \\oint_{\\partial\\Omega}\\left(u(\\vec{u}\\cdot\\hat{n}) - \\tau_{xx}\\right)dS + \\int_{\\partial\\Omega}P\\hat{x}\\cdot d\\vec{S},\\newline 
        
           F_u=\\int_\\Omega\\left(\\frac{\\partial u}{\\partial t}\\right)dV - \\int_{\\partial\\Omega_1} \\left(u^2 - 2\\nu\\frac{\\partial u}{\\partial x} + P\\right)dS + \\int_{\\partial\\Omega_2}\\left(u^2 - 2\\nu\\frac{\\partial u}{\\partial x}+P\\right)dS -  \\int_{\\partial\\Omega_2} uvdS + \\int_{\\partial\\Omega_4} uvdS, 
        
           ``` 
        
           where the bounding box is ``\\Omega`` which is formed from the boundaries ``\\partial\\Omega_{1}``, ``\\partial\\Omega_{2}``, ``\\partial\\Omega_{3}``, and ``\\partial\\Omega_{4}``  
        
           which have outward directed normals ``-\\hat{x}``, ``\\hat{x}``, ``-\\hat{y}``, and ``\\hat{y}`` 
        
           """ 
        
           function drag(model; 
        
                         bounding_box = (-1, 3, -2, 2), 
        
                         ν = 1e-3) 
        
               u, v, _ = model.velocities 
        
               uᶜ = Field(@at (Center, Center, Center) u) 
        
               vᶜ = Field(@at (Center, Center, Center) v) 
        
               xc, yc, _ = nodes(uᶜ) 
        
               i₁ = findfirst(xc .> bounding_box[1]) 
        
               i₂ = findlast(xc .< bounding_box[2]) 
        
               j₁ = findfirst(yc .> bounding_box[3]) 
        
               j₂ = findlast(yc .< bounding_box[4]) 
        
               uₗ² = Field(uᶜ^2, indices = (i₁, j₁:j₂, 1)) 
        
               uᵣ² = Field(uᶜ^2, indices = (i₂, j₁:j₂, 1)) 
        
               uvₗ = Field(uᶜ*vᶜ, indices = (i₁:i₂, j₁, 1)) 
        
               uvᵣ = Field(uᶜ*vᶜ, indices = (i₁:i₂, j₂, 1)) 
        
               ∂₁uₗ = Field(∂x(uᶜ), indices = (i₁, j₁:j₂, 1)) 
        
               ∂₁uᵣ = Field(∂x(uᶜ), indices = (i₂, j₁:j₂, 1)) 
        
               ∂ₜuᶜ = Field(@at (Center, Center, Center) model.timestepper.Gⁿ.u) 
        
               ∂ₜu = Field(∂ₜuᶜ, indices = (i₁:i₂, j₁:j₂, 1)) 
        
               p = model.pressures.pNHS 
        
               ∫∂ₓp = Field(∂x(p), indices = (i₁:i₂, j₁:j₂, 1)) 
        
               a_local = Field(Integral(∂ₜu)) 
        
               a_flux = Field(Integral(uᵣ²)) - Field(Integral(uₗ²)) + Field(Integral(uvᵣ)) - Field(Integral(uvₗ)) 
        
               a_viscous_stress = 2ν * (Field(Integral(∂₁uᵣ)) - Field(Integral(∂₁uₗ))) 
        
               a_pressure = Field(Integral(∫∂ₓp)) 
        
               return a_local + a_flux + a_pressure - a_viscous_stress 
        
           end

could belong in Oceanostics?

jagoosw · 2024-12-06T19:04:33Z

I have also run the same case with different domain lengths (6m, 12m - shown, 18m and 30m) and they all produce very similar Cd and St

tomchor · 2024-12-06T19:09:13Z

@tomchor, do you think something like this:

Oceananigans.jl/validation/open_boundaries/cylinder.jl

Lines 18 to 76 in 32cd354

"""

drag(modell; bounding_box = (-1, 3, -2, 2), ν = 1e-3)

Returns the drag within the `bounding_box` computed by:

```math

\\frac{\\partial \\vec{u}}{\\partial t} + (\\vec{u}\\cdot\\nabla)\\vec{u}=-\\nabla P + \\nabla\\cdot\\vec{\\tau} + \\vec{F},\\newline

\\vec{F}_T=\\int_\\Omega\\vec{F}dV = \\int_\\Omega\\left(\\frac{\\partial \\vec{u}}{\\partial t} + (\\vec{u}\\cdot\\nabla)\\vec{u}+\\nabla P - \\nabla\\cdot\\vec{\\tau}\\right)dV,\\newline

\\vec{F}_T=\\int_\\Omega\\left(\\frac{\\partial \\vec{u}}{\\partial t}\\right)dV + \\oint_{\\partial\\Omega}\\left(\\vec{u}(\\vec{u}\\cdot\\hat{n}) + P\\hat{n} - \\vec{\\tau}\\cdot\\hat{n}\\right)dS,\\newline

\\F_u=\\int_\\Omega\\left(\\frac{\\partial u}{\\partial t}\\right)dV + \\oint_{\\partial\\Omega}\\left(u(\\vec{u}\\cdot\\hat{n}) - \\tau_{xx}\\right)dS + \\int_{\\partial\\Omega}P\\hat{x}\\cdot d\\vec{S},\\newline

F_u=\\int_\\Omega\\left(\\frac{\\partial u}{\\partial t}\\right)dV - \\int_{\\partial\\Omega_1} \\left(u^2 - 2\\nu\\frac{\\partial u}{\\partial x} + P\\right)dS + \\int_{\\partial\\Omega_2}\\left(u^2 - 2\\nu\\frac{\\partial u}{\\partial x}+P\\right)dS - \\int_{\\partial\\Omega_2} uvdS + \\int_{\\partial\\Omega_4} uvdS,

```

where the bounding box is ``\\Omega`` which is formed from the boundaries ``\\partial\\Omega_{1}``, ``\\partial\\Omega_{2}``, ``\\partial\\Omega_{3}``, and ``\\partial\\Omega_{4}``

which have outward directed normals ``-\\hat{x}``, ``\\hat{x}``, ``-\\hat{y}``, and ``\\hat{y}``

"""

function drag(model;

bounding_box = (-1, 3, -2, 2),

ν = 1e-3)

u, v, _ = model.velocities

uᶜ = Field(@at (Center, Center, Center) u)

vᶜ = Field(@at (Center, Center, Center) v)

xc, yc, _ = nodes(uᶜ)

i₁ = findfirst(xc .> bounding_box[1])

i₂ = findlast(xc .< bounding_box[2])

j₁ = findfirst(yc .> bounding_box[3])

j₂ = findlast(yc .< bounding_box[4])

uₗ² = Field(uᶜ^2, indices = (i₁, j₁:j₂, 1))

uᵣ² = Field(uᶜ^2, indices = (i₂, j₁:j₂, 1))

uvₗ = Field(uᶜ*vᶜ, indices = (i₁:i₂, j₁, 1))

uvᵣ = Field(uᶜ*vᶜ, indices = (i₁:i₂, j₂, 1))

∂₁uₗ = Field(∂x(uᶜ), indices = (i₁, j₁:j₂, 1))

∂₁uᵣ = Field(∂x(uᶜ), indices = (i₂, j₁:j₂, 1))

∂ₜuᶜ = Field(@at (Center, Center, Center) model.timestepper.Gⁿ.u)

∂ₜu = Field(∂ₜuᶜ, indices = (i₁:i₂, j₁:j₂, 1))

p = model.pressures.pNHS

∫∂ₓp = Field(∂x(p), indices = (i₁:i₂, j₁:j₂, 1))

a_local = Field(Integral(∂ₜu))

a_flux = Field(Integral(uᵣ²)) - Field(Integral(uₗ²)) + Field(Integral(uvᵣ)) - Field(Integral(uvₗ))

a_viscous_stress = 2ν * (Field(Integral(∂₁uᵣ)) - Field(Integral(∂₁uₗ)))

a_pressure = Field(Integral(∫∂ₓp))

return a_local + a_flux + a_pressure - a_viscous_stress

end

could belong in Oceanostics?

Sure! Why not? :)

tomchor · 2024-12-06T19:14:35Z

src/BoundaryConditions/perturbation_advection_open_boundary_matching_scheme.jl

+@inline function _fill_east_halo!(j, k, grid, u, bc::PAOBC, loc::Tuple{Face, Any, Any}, clock, model_fields)
+    i = grid.Nx + 1
+
+    Δt = clock.last_stage_Δt
+
+    Δt = ifelse(isinf(Δt), 0, Δt)
+
+    Δx = xspacing(i, j, k, grid, loc...)
+
+    ūⁿ⁺¹ = getbc(bc, j, k, grid, clock, model_fields)
+
+    uᵢⁿ     = @inbounds u[i, j, k]
+    uᵢ₋₁ⁿ⁺¹ = @inbounds u[i - 1, j, k]
+
+    U = max(0, min(1, Δt / Δx * ūⁿ⁺¹))
+
+    τ = ifelse(ūⁿ⁺¹ >= 0, 
+               bc.classification.matching_scheme.outflow_timescale, 
+               bc.classification.matching_scheme.inflow_timescale)
+
+
+    τ̃ = Δt / τ
+
+    uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + U)
+
+    @inbounds u[i, j, k] = uᵢⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid))
+end


I think we could implement everything using getindex() so that we implement it once (or at max twice; one for left and other for right BC) and can use it in three dimensions.

For example this

uᵢⁿ = @inbounds u[i, j, k] uᵢ₋₁ⁿ⁺¹ = @inbounds u[i - 1, j, k]

would become

uᵢⁿ = @inbounds getindex(u, boundary_indices...) uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, one_off_boundary_indices...)

This would avoid errors when transcribing to other dimensions (I remember catching a couple for the flat extrapolation matching scheme).

yeah, I was trying to think of a clean way to do this

Attempted this now (but I've still not written out the other directions in case we want to change things first)

validation/open_boundaries/cylinder.jl

glwagner · 2024-12-06T20:45:15Z

There's a typo in your top equation right? I think it should be $\partial_t (u + U) = U \partial_x u$.

glwagner · 2024-12-06T20:46:30Z

A problem I think this scheme might have is that as the flow changes direction the assumption that the perturbation is small falls apart, which is maybe why it's been a bit unstable so far in an oscillating case.

What do you mean? What is the consequence of this approximation?

jagoosw · 2024-12-06T20:57:23Z

What do you mean? What is the consequence of this approximation?

I was thinking about this more and really to remain valid then in :

$$\partial_t (u\prime+U)\approx u\prime\partial_xu\prime + U\partial_xu\prime + (U - u\prime)/\tau$$

as $U\to0$ then we need $u\prime\partial_xu\prime\ll u\prime/\tau$ to still be able to throw away the first term, so I think its actually okay as long as you have the relaxation.

I found a typo in the west boundary which might be why the oscillating case wasn't working, rerunning now

jagoosw · 2024-12-06T20:57:31Z

There's a typo in your top equation right? I think it should be ∂ t ( u + U ) = U ∂ x u .
Yeah thanks for catching!

glwagner · 2024-12-06T21:57:59Z

What do you mean? What is the consequence of this approximation?

I was thinking about this more and really to remain valid then in :

as U → 0 then we need u ′ ∂ x u ′ ≪ u ′ / τ to still be able to throw away the first term, so I think its actually okay as long as you have the relaxation.

I found a typo in the west boundary which might be why the oscillating case wasn't working, rerunning now

Okay, I agree with you. You can also analyze the update formula / backward Euler step in the limit U -> 0. It doesn't look like it needs to be regularized in any way in that case though, it looks ok.

tomchor · 2024-12-09T19:14:17Z

Is there a reference for this method? If there is, please include it on the methods docstring. If not, it'd be good to include a brief explanation like the one at the top comment in the dosctring. I don't think this radiation method is very well known.

jagoosw · 2025-01-18T18:06:23Z

Also feel free to mark this as ready for review once you add tests.

Should be ready now

jagoosw · 2025-01-20T15:40:53Z

Boundary condition tests passed just some GPU tests failed because they lost Sverdrup

tomchor · 2025-01-21T02:57:41Z

src/boundary_mean.jl

@jagoosw you still need to remove this file from here, right?

Is this only used for the perturbation OBC? If so, it can go in the same directory.

We never resolved this, the problem is that AbstractOperations are defined after BoundaryConditions so this can't go in the boundary module

Hmm, I see. I don't have a good answer here, but I don't think it can be in /src. @glwagner @simone-silvestri can you advise here?

Also, is there a simple way to define these functions without AbstractOperations?

It would be possible to write it without but it would not be trivial to do it for general grids and would end up reinventing the wheel to get the flux.

Perhaps we could move it to Utils?

I'm okay with that.

src/boundary_mean.jl

tomchor · 2025-01-21T03:01:00Z

test/test_boundary_conditions_integration.jl

+
+        u = normal_velocity(Val(orientation), model)
+
+        CUDA.@allowscalar u.data .= 1


I think we're trying to avoid introducing new @alloscalar instances to testing. @glwagner can you confirm?

I thought the problem was that tests currently have a global @allowscalar, but some tests need them or will be more complicated to write?

I think you can write most of them simply without allowscalar, but I'm not an expert on this. For example you may be able to write something like interior(u) .= 1 without it (I haven't tested it though), and it's a pretty simple change.

I say try to avoid it with simple statements like the one I just mentioned. If we can't do it simply and @glwagner doesn't manifest any opinions on it, we move forward. How does that sound?

I sorted this now, it wasn't as straight forward as just using interior because I needed to set a peripheral node, but I just used view instead

test/test_boundary_conditions_integration.jl

validation/open_boundaries/cylinder.jl

tomchor

Made a few comments, but it looks like there are old unresolved comments and some files here and there which do not seem to be ready for prime yet. let me know if anything needs clarification.

Co-authored-by: Tomas Chor <[email protected]>

jagoosw · 2025-01-21T15:40:35Z

@tomchor I think I've addressed your comments now

validation/open_boundaries/cylinder.jl

tomchor · 2025-01-21T15:49:08Z

validation/open_boundaries/cylinder.jl


-ax = Axis(fig[1, 1], aspect = DataAspect())
+    simulation.output_writers[:drag] = JLD2OutputWriter(model, (; drag_force),
+                                                        schedule = IterationInterval(Int(0.1/Δt)),#TimeInterval(0.1),


Suggested change

schedule = IterationInterval(Int(0.1/Δt)),#TimeInterval(0.1),

schedule = IterationInterval(Int(0.1/Δt)),

tomchor · 2025-01-21T15:49:23Z

validation/open_boundaries/cylinder.jl

+    outputs = (; u, v, p, ζ)
+
+    simulation.output_writers[:jld2] = JLD2OutputWriter(model, outputs,
+                                                        schedule = IterationInterval(Int(2/Δt)),#TimeInterval(0.1),


Suggested change

schedule = IterationInterval(Int(2/Δt)),#TimeInterval(0.1),

schedule = IterationInterval(Int(2/Δt)),

tomchor · 2025-01-21T15:49:38Z

validation/open_boundaries/cylinder.jl


-progress(sim) = @info "$(time(sim)) with Δt = $(prettytime(sim.Δt)) in $(prettytime(sim.run_wall_time))"
+    simulation = Simulation(model; Δt, stop_time)
+    #conjure_time_step_wizard!(simulation, cfl=1.0, IterationInterval(3); max_Δt)


Suggested change

#conjure_time_step_wizard!(simulation, cfl=1.0, IterationInterval(3); max_Δt)

tomchor · 2025-01-21T15:56:00Z

validation/open_boundaries/cylinder.jl

+# there is some problem using ConjugateGradientPoissonSolver with TimeInterval because the timestep can go really small
+# so while I identify the issue I'm using IterationInterval and a fixed timestep


Have you tried setting minimum_relative_step? (see #3606)

This seems to fix the issue! I hadn't tried this since early november/december so I hadn't seen you'd made this fix so this is great, thank you

I'll update the validation script

Co-authored-by: Tomas Chor <[email protected]>

glwagner · 2025-01-22T20:52:51Z

validation/open_boundaries/cylinder.jl

+F⃗ₜ = ∫ᵥF⃗dV = ∫ᵥ(∂ₜ u⃗ + (u⃗⋅∇)u⃗ + ∇P − ∇⋅τ⃗)dV
+F⃗ₜ = ∫ᵥ(∂ₜu⃗)dV + ∮ₛ(u⃗(u⃗⋅n̂) + Pn̂ − τ⃗⋅n̂)dS
+Fᵤ = ∫ᵥ ∂ₜu dV + ∮ₛ(u(u⃗⋅n̂) − τₓₓ)dS + ∮ₛPx̂⋅dS⃗
+Fᵤ = ∫ᵥ ∂ₜ u dV − ∫ₛ₁(u² − 2ν∂ₓ u + P)dS + ∫ₛ₂(u² − 2ν∂ₓ u+P)dS − ∫ₛ₃uvdS + ∫ₛ₄ uvdS


FYI what this looks like on my screen / browser

I looks similarly messed up on my browser, but the actual document looks fine on my laptop:

Probably it's best to prioritize rendering on github, since a lot of people use it. I guess the easiest way to do that is to copy-paste stuff on a github text box and see what displays correctly and what doesn't. I just tested it out right now and it seems to work. The \oint character seems to not work for example.

Right. Arguably everyone uses github btw. Eg we are using github now.

Probably it's best to prioritize rendering on github, since a lot of people use it. I guess the easiest way to do that is to copy-paste stuff on a github text box and see what displays correctly and what doesn't. I just tested it out right now and it seems to work. The \oint character seems to not work for example.

I tried this but then when I copy it over it renders wrong:
$$\frac{\partial \vec{u}}{\partial t} + (\vec{u}\cdot\nabla)\vec{u}=-\nabla P + \nabla\cdot\vec{\tau} + \vec{F},$$ $$\vec{F}_T=\int_\Omega\vec{F}dV = \int_\Omega\left(\frac{\partial \vec{u}}{\partial t} + (\vec{u}\cdot\nabla)\vec{u}+\nabla P - \nabla\cdot\vec{\tau}\right)dV,\newline$$ $$\vec{F}_T=\int_\Omega\left(\frac{\partial \vec{u}}{\partial t}\right)dV + \oint_{\partial\Omega}\left(\vec{u}(\vec{u}\cdot\hat{n}) + P\hat{n} - \vec{\tau}\cdot\hat{n}\right)dS,$$ $$F_u=\int_\Omega\left(\frac{\partial u}{\partial t}\right)dV + \oint_{\partial\Omega}\left(u(\vec{u}\cdot\hat{n}) - \tau_{xx}\right)dS + \int_{\partial\Omega}P\hat{x}\cdot d\vec{S},$$ $$F_u=\int_\Omega\left(\frac{\partial u}{\partial t}\right)dV - \int_{\partial\Omega_1} \left(u^2 - 2\nu\frac{\partial u}{\partial x} + P\right)dS + \int_{\partial\Omega_2}\left(u^2 - 2\nu\frac{\partial u}{\partial x}+P\right)dS - \int_{\partial\Omega_2} uvdS + \int_{\partial\Omega_4} uvdS$$
where the bounding box is $\Omega$ which is formed from the boundaries $\partial\Omega_1$, $\partial\Omega_2$, $\partial\Omega_3$, $\partial\Omega_4$ which have outward directed normals $-\hat{x}$, $\hat{x}$, $-\hat{y}$ and $\hat{y}$

Rendered in julia file:
∂ u → ∂ t + ( u → ⋅ ∇ ) u → = − ∇ P + ∇ ⋅ τ → + F → ,
F → T = ∫ Ω F → d V = ∫ Ω ( ∂ u → ∂ t + ( u → ⋅ ∇ ) u → + ∇ P − ∇ ⋅ τ → ) d V ,
F → T = ∫ Ω ( ∂ u → ∂ t ) d V + ∮ ∂ Ω ( u → ( u → ⋅ n ^ ) + P n ^ − τ → ⋅ n ^ ) d S ,
F u = ∫ Ω ( ∂ u ∂ t ) d V + ∮ ∂ Ω ( u ( u → ⋅ n ^ ) − τ x x ) d S + ∫ ∂ Ω P x ^ ⋅ d S → ,
F u = ∫ Ω ( ∂ u ∂ t ) d V − ∫ ∂ Ω 1 ( u 2 − 2 ν ∂ u ∂ x + P ) d S + ∫ ∂ Ω 2 ( u 2 − 2 ν ∂ u ∂ x + P ) d S − ∫ ∂ Ω 2 u v d S + ∫ ∂ Ω 4 u v d S

where the bounding box is Ω which is formed from the boundaries ∂ Ω 1 , ∂ Ω 2 , ∂ Ω 3 , ∂ Ω 4 which have outward directed normals − x ^ , x ^ , − y ^ and y

I'll play around and try and find an intermediary

How does this look:

Suggested change

Fᵤ = ∫ᵥ ∂ₜ u dV − ∫ₛ₁(u² − 2ν∂ₓ u + P)dS + ∫ₛ₂(u² − 2ν∂ₓ u+P)dS − ∫ₛ₃uvdS + ∫ₛ₄ uvdS

∂ₜU + (U⋅∇)U = −∇P + ∇⋅τ + F⃗

F = ∫ᵥFdV = ∫ᵥ(∂U/∂t + (U⋅∇)U + ∇P − ∇⋅τ)dV

F = ∫ᵥ(∂U/∂t)dV + ∫ₛ(U(U⋅n̂) + Pn̂ − τ⋅n̂)dS

Fᵤ = ∫ᵥ ∂U/∂t dV + ∫ₛ(u(U⋅n̂) − τₓₓ)dS + ∮ₛPx̂⋅dS⃗

Fᵤ = ∫ᵥ ∂U/∂t dV − ∫ₛ₁(u² − 2ν∂ₓ u + P)dS + ∫ₛ₂(u² − 2ν∂ₓ u+P)dS − ∫ₛ₃uvdS + ∫ₛ₄ uvdS

I think I got all of the unrendered characters gone (vectors and oint?).

Strange that GH renders the original fine for me, maybe because I have a systemwide juliamono installation?

not sure. But your update renders fine.

jagoosw and others added 3 commits November 18, 2024 17:32

fix FEOBC

d992d81

x-boundary condition implemented

2a374cb

changed relaxation a bit

babe922

jagoosw added feature 🌟 Something new and shiny science 🌊 Sometimes addictive, sometimes allergenic numerics 🧮 So things don't blow up and boil the lobsters alive labels Dec 6, 2024

cleaning up + validation case

32cd354

tomchor reviewed Dec 6, 2024

View reviewed changes

jagoosw mentioned this pull request Dec 6, 2024

Unexpected FlatExtrapolationBoundaryCondition behaviour in validation script cylinder.jl #3975

Closed

jagoosw and others added 4 commits December 6, 2024 14:18

changed to pointwise

5a0285f

oops

fef6185

fixed open boundary filling + tests

fa5a471

Merge branch 'main' into jsw/fix-feobc

29a3b4a

jagoosw force-pushed the jsw/pertubation-advection-obc branch from 0e645c5 to 32cd354 Compare December 6, 2024 20:11

jagoosw added 2 commits December 6, 2024 15:33

generalised and fixed bug in left boundary

b55395c

bump patch

1a9fbec

jagoosw commented Dec 6, 2024

View reviewed changes

validation/open_boundaries/cylinder.jl Outdated Show resolved Hide resolved

Formating

e66f3d1

jagoosw force-pushed the jsw/pertubation-advection-obc branch from afecfb3 to e66f3d1 Compare December 6, 2024 20:43

jagoosw and others added 2 commits December 10, 2024 15:33

Merge branch 'main' into jsw/fix-feobc

dc73eba

test problem

51c9938

Merge branch 'main' into jsw/pertubation-advection-obc

dddb91a

tomchor reviewed Jan 21, 2025

View reviewed changes

src/boundary_mean.jl Outdated Show resolved Hide resolved

tomchor reviewed Jan 21, 2025

View reviewed changes

test/test_boundary_conditions_integration.jl Outdated Show resolved Hide resolved

tomchor reviewed Jan 21, 2025

View reviewed changes

test/test_boundary_conditions_integration.jl Outdated Show resolved Hide resolved

tomchor reviewed Jan 21, 2025

View reviewed changes

validation/open_boundaries/cylinder.jl Outdated Show resolved Hide resolved

tomchor reviewed Jan 21, 2025

View reviewed changes

validation/open_boundaries/cylinder.jl Outdated Show resolved Hide resolved

tomchor requested changes Jan 21, 2025

View reviewed changes

jagoosw and others added 3 commits January 21, 2025 15:08

Apply suggestions from code review

00d56cf

Co-authored-by: Tomas Chor <[email protected]>

commented tests

0de469d

tidied up validation

270d26a

tomchor reviewed Jan 21, 2025

View reviewed changes

validation/open_boundaries/cylinder.jl Outdated Show resolved Hide resolved

tomchor reviewed Jan 21, 2025

View reviewed changes

jagoosw added 2 commits January 22, 2025 13:20

fixed adapt

f20894e

fix validation case

ef4d4d7

jagoosw changed the title ~~[WIP] Adds a pertubation advection open boundary matching scheme~~ Adds a pertubation advection open boundary matching scheme Jan 22, 2025

jagoosw and others added 6 commits January 22, 2025 15:47

Update validation/open_boundaries/cylinder.jl

2b15848

Co-authored-by: Tomas Chor <[email protected]>

added plotting to validation

d12414c

added explanation to docstring

335926e

corrected validation script

c3241ca

Merge branch 'main' into jsw/pertubation-advection-obc

fa94dd1

removed @allowscalars

35f42b2

glwagner reviewed Jan 22, 2025

View reviewed changes


		u = normal_velocity(Val(orientation), model)

		CUDA.@allowscalar u.data .= 1

	schedule = IterationInterval(Int(0.1/Δt)),#TimeInterval(0.1),
	schedule = IterationInterval(Int(0.1/Δt)),

	schedule = IterationInterval(Int(2/Δt)),#TimeInterval(0.1),
	schedule = IterationInterval(Int(2/Δt)),

		# there is some problem using ConjugateGradientPoissonSolver with TimeInterval because the timestep can go really small
		# so while I identify the issue I'm using IterationInterval and a fixed timestep

-Fᵤ = ∫ᵥ ∂ₜ u dV − ∫ₛ₁(u² − 2ν∂ₓ u + P)dS + ∫ₛ₂(u² − 2ν∂ₓ u+P)dS − ∫ₛ₃uvdS + ∫ₛ₄ uvdS
+∂ₜU + (U⋅∇)U = −∇P + ∇⋅τ + F⃗
+F = ∫ᵥFdV = ∫ᵥ(∂U/∂t + (U⋅∇)U + ∇P − ∇⋅τ)dV
+F = ∫ᵥ(∂U/∂t)dV + ∫ₛ(U(U⋅n̂) + Pn̂ − τ⋅n̂)dS
+Fᵤ = ∫ᵥ ∂U/∂t dV + ∫ₛ(u(U⋅n̂) − τₓₓ)dS + ∮ₛPx̂⋅dS⃗
+Fᵤ = ∫ᵥ ∂U/∂t dV − ∫ₛ₁(u² − 2ν∂ₓ u + P)dS + ∫ₛ₂(u² − 2ν∂ₓ u+P)dS − ∫ₛ₃uvdS + ∫ₛ₄ uvdS

Adds a pertubation advection open boundary matching scheme #3977

Are you sure you want to change the base?

Adds a pertubation advection open boundary matching scheme #3977

Conversation

jagoosw commented Dec 6, 2024 • edited Loading

jagoosw commented Dec 6, 2024

jagoosw commented Dec 6, 2024

tomchor commented Dec 6, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

glwagner commented Dec 6, 2024 • edited Loading

glwagner commented Dec 6, 2024

jagoosw commented Dec 6, 2024 • edited Loading

jagoosw commented Dec 6, 2024 • edited Loading

glwagner commented Dec 6, 2024 • edited Loading

tomchor commented Dec 9, 2024

jagoosw commented Jan 18, 2025

jagoosw commented Jan 20, 2025

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tomchor left a comment

Choose a reason for hiding this comment

jagoosw commented Jan 21, 2025

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jagoosw commented Dec 6, 2024 •

edited

Loading

glwagner commented Dec 6, 2024 •

edited

Loading

jagoosw commented Dec 6, 2024 •

edited

Loading

jagoosw commented Dec 6, 2024 •

edited

Loading

glwagner commented Dec 6, 2024 •

edited

Loading