diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 6f7a235..38789fb 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -1,3 +1,6 @@ +# Tests on GPU back-ends. Draft pull requests skip them, unless selected by name in the last commit message: `[only cuda]` +# or `[only tests]` for all of them. `[skip amdgpu]` or `[skip tests]` does the opposite. + steps: - label: "CUDA - Julia {{matrix.version}}" matrix: @@ -7,11 +10,15 @@ steps: plugins: - JuliaCI/julia#v1: version: "{{matrix.version}}" - - JuliaCI/julia-test#v1: ~ + #- JuliaCI/julia-test#v1: ~ + command: julia --eval='println("Skipping CUDA tests - pipeline configured to do nothing")' agents: queue: "cuda" - if: build.message !~ /\[skip tests\]/ + env: + WITH_CUDA_TESTS: "true" + WITH_CPU_TESTS: "false" + # if: | + # build.message !~ /\[skip [^\]]*(tests|cuda)/ && + # (build.message =~ /\[only [^\]]*(tests|cuda)/ || + # build.message !~ /\[only / && !build.pull_request.draft) timeout_in_minutes: 60 - -env: - JULIA_PKG_SERVER_REGISTRY_PREFERENCE: "eager" diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..f64b162 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,62 @@ +name: CI +on: + push: + branches: + - main + pull_request: + types: [unlabeled, opened, synchronize, reopened] + workflow_dispatch: + + +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} + +# Only trigger the job when `draft` label is not assigned to the PR +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} + runs-on: ${{ matrix.os }} + if: contains(github.event.pull_request.labels.*.name, 'draft') == false + strategy: + matrix: + version: + - '1.10' + os: + - ubuntu-latest + arch: + - x64 + env: + GITHUB_CI: "true" + steps: + - uses: actions/checkout@v7 + with: + lfs: "true" + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: actions/cache@v6 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - uses: julia-actions/julia-buildpkg@latest + - uses: julia-actions/julia-runtest@latest + env: + WITH_CUDA_TESTS: "false" + WITH_CPU_TESTS: "true" + - uses: julia-actions/julia-processcoverage@latest + - name: Upload coverage reports to Codecov + uses: codecov/codecov-action@v7 + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + with: + file: lcov.info \ No newline at end of file diff --git a/.github/workflows/format_check.yml b/.github/workflows/format_check.yml new file mode 100644 index 0000000..ba48d9d --- /dev/null +++ b/.github/workflows/format_check.yml @@ -0,0 +1,15 @@ +## Check that the code in a PR is correctly formatted +## from https://github.com/julia-actions/julia-format +name: Format check +on: + pull_request: + types: [opened, reopened, synchronize, labeled, unlabeled] + workflow_dispatch: +jobs: + code-style: + runs-on: ubuntu-latest + steps: + - uses: julia-actions/julia-format@v3 + with: + version: "2" + suggestion-label: "format-suggest" \ No newline at end of file diff --git a/JuliaFormat.toml b/JuliaFormat.toml new file mode 100644 index 0000000..d14d08d --- /dev/null +++ b/JuliaFormat.toml @@ -0,0 +1,59 @@ +# Configuration file for code formatting +# see here for default options : https://github.com/domluna/JuliaFormatter.jl#formatting-options +style = "default" + +[options] +always_for_in = true +whitespace_typedefs = true +whitespace_ops_in_indices = true +remove_extra_newlines = true +short_to_long_function_def = true +long_to_short_function_def = false +# always_use_return = true +whitespace_in_kwargs = true +format_docstrings = false +conditional_to_if = true +trailing_comma = true +normalize_line_endings = "unix" +join_lines_based_on_source = false +separate_kwargs_with_semicolon = true +align_assignment = true +align_struct_field = true +align_conditional = true +# align_pair_arrow = true +# align_matrix = true +# format_markdown = true + + +# possible style : sciml, blue, yas + +# ### SciML options +# whitespace_ops_in_indices = true +# remove_extra_newlines = true +# always_for_in = true +# whitespace_typedefs = true +# normalize_line_endings = "unix" + +# ### YAS options +# always_for_in = true +# whitespace_ops_in_indices = true +# remove_extra_newlines = true +# import_to_using = true +# pipe_to_function_call = true +# short_to_long_function_def = true +# always_use_return = true +# whitespace_in_kwargs = false +# join_lines_based_on_source = true +# separate_kwargs_with_semicolon = true + +# ### Blue option +# always_use_return = true +# short_to_long_function_def = true +# whitespace_ops_in_indices = true +# remove_extra_newlines = true +# always_for_in = true +# import_to_using = true +# pipe_to_function_call = true +# whitespace_in_kwargs = false +# annotate_untyped_fields_with_any = false +# separate_kwargs_with_semicolon = true \ No newline at end of file diff --git a/src/BcubeAccelerated.jl b/src/BcubeAccelerated.jl index 3c6b469..6e76672 100644 --- a/src/BcubeAccelerated.jl +++ b/src/BcubeAccelerated.jl @@ -17,4 +17,4 @@ include("fespace.jl") include("domain.jl") include("assembler.jl") -end \ No newline at end of file +end diff --git a/src/adapt.jl b/src/adapt.jl index c3216bc..4dbedcb 100644 --- a/src/adapt.jl +++ b/src/adapt.jl @@ -197,4 +197,4 @@ end Adapt.adapt_structure(to::AbstractBcubeBackendAcc, x) = adapt_structure(get_backend(to), x) Adapt.adapt(to::Bcube.AbstractBcubeBackend, x) = x -Adapt.adapt(to::AbstractBcubeBackendAcc, x) = adapt(get_backend(to), x) \ No newline at end of file +Adapt.adapt(to::AbstractBcubeBackendAcc, x) = adapt(get_backend(to), x) diff --git a/src/domain.jl b/src/domain.jl index 7aee3d7..ddbc51e 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -119,4 +119,4 @@ function Bcube.get_return_type_and_codim( _T = Base.return_types(i -> _f(i, valSize), (Int,)) T = eltype(first(_T)) return T, N -end \ No newline at end of file +end diff --git a/test/covo_gpu.jl b/test/covo_gpu.jl index 669b3ec..31d6734 100644 --- a/test/covo_gpu.jl +++ b/test/covo_gpu.jl @@ -1,515 +1,516 @@ -module CovoGpu #hide -println("Running covo example...") #hide - -const dir = string(@__DIR__, "/") -using Bcube -using BcubeAccelerated -using KernelAbstractions -using CUDA, CUDA.CUSPARSE, CUDA.CUSOLVER -using BcubeGmsh -using Adapt -using BcubeVTK -using LinearAlgebra -using StaticArrays -using StaticArrays -using BenchmarkTools -using UnPack - -const WITH_VTK_OUTPUT = false - -function compute_residual(_u, V, params, cache) - u = get_fe_functions(_u) - - # alias on measures - @unpack dΩ, dΓ, dΓ_perio_x, dΓ_perio_y = params - - # face normals for each face domain (lazy, no computation at this step) - nΓ = get_face_normals(dΓ) - nΓ_perio_x = get_face_normals(dΓ_perio_x) - nΓ_perio_y = get_face_normals(dΓ_perio_y) - - # flux residuals from faces for all variables - function l(v) - ∫(flux_Ω(u, v))dΩ + - -∫(flux_Γ(u, v, nΓ))dΓ + - -∫(flux_Γ(u, v, nΓ_perio_x))dΓ_perio_x + - -∫(flux_Γ(u, v, nΓ_perio_y))dΓ_perio_y - end - # l1(v) = ∫(flux_Ω(u, v))dΩ - # l2(v) = ∫(flux_Γ(u, v, nΓ))dΓ - # l3(v) = ∫(flux_Γ(u, v, nΓ_perio_x))dΓ_perio_x - # l4(v) = ∫(flux_Γ(u, v, nΓ_perio_y))dΓ_perio_y - # @show "l1" - # assemble_linear(l1, V) - # @show "l2" - # assemble_linear(l2, V) - # @show "l3" - # assemble_linear(l3, V) - # @show "l4" - # assemble_linear(l4, V) - # error("ici") - _rhs = assemble_linear(l, V) - rhs = similar(_rhs) - Bcube.__solve!(rhs, cache.mass, _rhs, Bcube.get_bcube_backend(dΩ)) - return rhs -end - -""" - flux_Ω(u, v) - -Compute volume residual using the lazy-operators approach -""" -flux_Ω(u, v) = _flux_Ω ∘ cellvar(u, v) -cellvar(u, v) = (u, map(∇, v)) -function _flux_Ω(u, ∇v) - ρ, ρu, ρE, ϕ = u - ∇λ_ρ, ∇λ_ρu, ∇λ_ρE, ∇λ_ϕ = ∇v - - vel = ρu ./ ρ - ρuu = ρu * transpose(vel) - p = pressure(ρ, ρu, ρE, γ) - - flux_ρ = ρu - flux_ρu = ρuu + p * I - flux_ρE = (ρE + p) .* vel - flux_ϕ = ϕ .* vel - - return ∇λ_ρ ⋅ flux_ρ + ∇λ_ρu ⊡ flux_ρu + ∇λ_ρE ⋅ flux_ρE + ∇λ_ϕ ⋅ flux_ϕ -end - -""" - flux_Γ(u,v,n) - -Flux at the interface is defined by a composition of two functions: -* facevar(u,v,n) defines the input states which are needed for - the riemann flux using operator notations -* flux_roe(w) defines the Riemann flux (as usual) -""" -flux_Γ(u, v, n) = flux_roe ∘ (side⁻(u), side⁺(u), jump(v), side⁻(n)) - -""" - flux_roe(w) -""" -function flux_roe(ui, uj, δv, nij) - # destructuring inputs given by `facevar` function - - nx, ny = nij - ρ1, ρu1, ρE1, ϕ1 = ui - ρ2, ρu2, ρE2, ϕ2 = uj - δλ_ρ1, δλ_ρu1, δλ_ρE1, δλ_ϕ1 = δv - ρux1, ρuy1 = ρu1 - ρux2, ρuy2 = ρu2 - - # Closure - u1 = ρux1 / ρ1 - v1 = ρuy1 / ρ1 - u2 = ρux2 / ρ2 - v2 = ρuy2 / ρ2 - p1 = pressure(ρ1, ρu1, ρE1, γ) - p2 = pressure(ρ2, ρu2, ρE2, γ) - - H2 = (γ / (γ - 1)) * p2 / ρ2 + (u2 * u2 + v2 * v2) / 2.0 - H1 = (γ / (γ - 1)) * p1 / ρ1 + (u1 * u1 + v1 * v1) / 2.0 - - R = √(ρ1 / ρ2) - invR1 = 1.0 / (R + 1) - uAv = (R * u1 + u2) * invR1 - vAv = (R * v1 + v2) * invR1 - Hav = (R * H1 + H2) * invR1 - cAv = √(abs((γ - 1) * (Hav - (uAv * uAv + vAv * vAv) / 2.0))) - ecAv = (uAv * uAv + vAv * vAv) / 2.0 - - λ1 = nx * uAv + ny * vAv - λ3 = λ1 + cAv - λ4 = λ1 - cAv - - d1 = ρ1 - ρ2 - d2 = ρ1 * u1 - ρ2 * u2 - d3 = ρ1 * v1 - ρ2 * v2 - d4 = ρE1 - ρE2 - - # computation of the centered part of the flux - flu11 = nx * ρ2 * u2 + ny * ρ2 * v2 - flu21 = nx * p2 + flu11 * u2 - flu31 = ny * p2 + flu11 * v2 - flu41 = H2 * flu11 - - # Temp variables - rc1 = (γ - 1) / cAv - rc2 = (γ - 1) / cAv / cAv - uq41 = ecAv / cAv + cAv / (γ - 1) - uq42 = nx * uAv + ny * vAv - - fdc1 = max(λ1, 0.0) * (d1 + rc2 * (-ecAv * d1 + uAv * d2 + vAv * d3 - d4)) - fdc2 = max(λ1, 0.0) * ((nx * vAv - ny * uAv) * d1 + ny * d2 - nx * d3) - fdc3 = - max(λ3, 0.0) * ( - (-uq42 * d1 + nx * d2 + ny * d3) / 2.0 + - rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0 - ) - fdc4 = - max(λ4, 0.0) * ( - (uq42 * d1 - nx * d2 - ny * d3) / 2.0 + - rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0 - ) - - duv1 = fdc1 + (fdc3 + fdc4) / cAv - duv2 = uAv * fdc1 + ny * fdc2 + (uAv / cAv + nx) * fdc3 + (uAv / cAv - nx) * fdc4 - duv3 = vAv * fdc1 - nx * fdc2 + (vAv / cAv + ny) * fdc3 + (vAv / cAv - ny) * fdc4 - duv4 = - ecAv * fdc1 + - (ny * uAv - nx * vAv) * fdc2 + - (uq41 + uq42) * fdc3 + - (uq41 - uq42) * fdc4 - - v₁₂ = 0.5 * ((u1 + u2) * nx + (v1 + v2) * ny) - fluxϕ = max(0.0, v₁₂) * ϕ1 + min(0.0, v₁₂) * ϕ2 - - return ( - δλ_ρ1 * (flu11 + duv1) + - δλ_ρu1 ⋅ (SA[flu21 + duv2, flu31 + duv3]) + - δλ_ρE1 * (flu41 + duv4) + - δλ_ϕ1 * (fluxϕ) - ) -end - -""" -Time integration of `f(q, t)` over a timestep `Δt`. -""" -forward_euler(q, f::Function, t, Δt) = get_dof_values(q) .+ Δt .* f(q, t) - -""" - rk3_ssp(q, f::Function, t, Δt) - -`f(q, t)` is the function to integrate. -""" -function rk3_ssp(q, f::Function, t, Δt) - stepper(q, t) = forward_euler(q, f, t, Δt) - _q0 = get_dof_values(q) - - _q1 = stepper(q, Δt) - - set_dof_values!(q, _q1) - _q2 = (3 / 4) .* _q0 .+ (1 / 4) .* stepper(q, t + Δt) - - set_dof_values!(q, _q2) - _q1 .= (1 / 3) * _q0 .+ (2 / 3) .* stepper(q, t + Δt / 2) - - return _q1 -end - -""" - pressure(ρ, ρu, ρE, γ) - -Computes pressure from perfect gaz law -""" -function pressure(ρ::Number, ρu::AbstractVector, ρE::Number, γ) - vel = ρu ./ ρ - ρuu = ρu * transpose(vel) - p = (γ - 1) * (ρE - tr(ρuu) / 2) - return p -end - -""" - Init field with a vortex (for the COVO test case) -""" -function covo!(q, dΩ) - - # Intermediate vars - Cₚ = γ * r / (γ - 1) - - r²(x) = ((x[1] .- xvc) .^ 2 + (x[2] .- yvc) .^ 2) ./ Rc^2 - # Temperature - T(x) = T₀ .- β^2 * U₀^2 / (2 * Cₚ) .* exp.(-r²(x)) - # Velocity - ux(x) = U₀ .- β * U₀ / Rc .* (x[2] .- yvc) .* exp.(-r²(x) ./ 2) - uy(x) = V₀ .+ β * U₀ / Rc .* (x[1] .- xvc) .* exp.(-r²(x) ./ 2) - # Density - ρ(x) = ρ₀ .* (T(x) ./ T₀) .^ (1.0 / (γ - 1)) - # momentum - ρu(x) = SA[ρ(x) * ux(x), ρ(x) * uy(x)] - # Energy - ρE(x) = ρ(x) * ((Cₚ / γ) .* T(x) + (ux(x) .^ 2 + uy(x) .^ 2) ./ 2) - # Passive scalar - ϕ(x) = Rc^2 * r²(x) < 0.01 ? exp(-r²(x) ./ 2) : 0.0 - - f = map(PhysicalFunction, (ρ, ρu, ρE, ϕ)) - projection_l2!(q, f, dΩ) - return nothing -end - -""" - Tiny struct to ease the VTK output -""" -mutable struct VtkHandler - basename::Any - mesh::Any - ite::Any - VtkHandler(basename, mesh) = new(basename, mesh, 0) -end - -""" - Write solution to vtk - Wrapper for `write_vtk` -""" -function append_vtk(vtk, vars, t, params) - ρ, ρu, ρE, ϕ = vars - - mesh = vtk.mesh - mesh_degree = 1 - vtk_degree = maximum(x -> get_degree(Bcube.get_function_space(get_fespace(x))), vars) - vtk_degree = max(1, mesh_degree, vtk_degree) - - # p = pressure(ρ, ρu, ρE, γ) - p = pressure ∘ (ρ, ρu, ρE, γ) - dict_vars_dg = Dict("rho" => ρ, "rhou" => ρu, "rhoE" => ρE, "phi" => ϕ, "p" => p) - write_file( - vtk.basename * "_DG.pvd", - mesh, - dict_vars_dg, - vtk.ite, - t; - discontinuous = true, - mesh_degree = vtk_degree, - collection_append = vtk.ite > 0, - ) - - # Update counter - vtk.ite += 1 -end - -function Adapt.adapt_structure(to, mesh::BcubeGmsh.GmshMetaData) - nothing -end -function Bcube.__solve!( - x, - A, - b, - backend::BcubeAccelerated.BcubeBackendAcc{K, CUDA.CUDABackend}, -) where {K} - F = CUSOLVER.SparseCholesky(A) - CUSOLVER.spcholesky_factorise(F, A, 1.e-12) - CUSOLVER.spcholesky_solve(F, b, x) - return nothing -end -function Bcube.__solve!( - x, - A::CUSOLVER.SparseCholesky, - b, - backend::BcubeAccelerated.BcubeBackendAcc{K, CUDA.CUDABackend}, -) where {K} - CUSOLVER.spcholesky_solve(A, b, x) - return nothing -end -# function Bcube.get_bcube_backend(a::CuArray) -# BcubeAccelerated.BcubeBackendAcc(KernelAbstractions.get_backend(a)) -# end - -# Settings -if get(ENV, "BenchmarkMode", "false") == "false" #hide - const cellfactor = 1 - const nx = 32 * cellfactor + 1 - const ny = 32 * cellfactor + 1 - const fspace = :Lagrange - const timeScheme = :ForwardEuler -else #hide - const nx = 128 + 1 #hide - const ny = 128 + 1 #hide - const fspace = :Lagrange - const timeScheme = :ForwardEuler -end #hide -const nperiod = 1 # number of turn -const CFL = 0.1 -const degree = 1 # FunctionSpace degree -const degquad = 2 * degree + 1 -const γ = 1.4 -const β = 0.2 # vortex intensity -const r = 287.15 # Perfect gaz constant -const T₀ = 300 # mean-flow temperature -const P₀ = 1e5 # mean-flow pressure -const M₀ = 0.5 # mean-flow mach number -const ρ₀ = 1.0 # mean-flow density -const xvc = 0.0 # x-center of vortex -const yvc = 0.0 # y-center of vortex -const Rc = 0.005 # Charasteristic vortex radius -const c₀ = √(γ * r * T₀) # Sound velocity -const U₀ = M₀ * c₀ # mean-flow velocity -const V₀ = 0.0 # mean-flow velocity -const ϕ₀ = 1.0 -const l = 0.05 # half-width of the domain -const Δt = CFL * 2 * l / (nx - 1) / ((1 + β) * U₀ + c₀) / (2 * degree + 1) -#const Δt = 5.e-7 -const nout = 100 # Number of time steps to save -const outputpath = "./covo/" -const output = joinpath(@__DIR__, outputpath, "covo_deg$degree") -const nite = Int(floor(nperiod * 2 * l / (U₀ * Δt))) + 1 - -function build_mass_matrix_test(U, V, dΩ::D) where {D <: AbstractMeasure} - _m(u, v) = ∫(u ⋅ v)dΩ - # @code_warntype _m(1,2) - @show typeof(_m(1, 2)) - BcubeAccelerated.test_arg_AK_foreachindex(get_backend(Bcube.get_bcube_backend(dΩ)), dΩ) - BcubeAccelerated.test_arg_AK_foreachindex( - get_backend(Bcube.get_bcube_backend(dΩ)), - _m(1, 2), - ) - return assemble_bilinear(_m, U, V) -end -function setup(backend, tmp_path, fs) - mesh_cpu = read_mesh(tmp_path) - rm(tmp_path) - U_sca_cpu = TrialFESpace(fs, mesh_cpu, :discontinuous; size = 1) # DG, scalar - U_vec_cpu = TrialFESpace(fs, mesh_cpu, :discontinuous; size = 2) # DG, vectoriel - U_cpu = MultiFESpace(U_sca_cpu, U_vec_cpu, U_sca_cpu, U_sca_cpu) - V_sca_cpu = TestFESpace(U_sca_cpu) - V_vec_cpu = TestFESpace(U_vec_cpu) - V_cpu = MultiFESpace(V_sca_cpu, V_vec_cpu, V_sca_cpu, V_sca_cpu) - - mesh = adapt(backend, mesh_cpu) - BcubeAccelerated.test_arg(backend, mesh) - U = adapt(backend, U_cpu) - V = adapt(backend, V_cpu) - BcubeAccelerated.test_arg(backend, U) - BcubeAccelerated.test_arg(backend, V) - return mesh, U, V, mesh_cpu -end - -_factorize(backend::Bcube.AbstractBcubeBackend, M) = __factorize(get_backend(backend), M) -__factorize(backend, M) = factorize(M) -function __factorize(backend::CUDA.CUDABackend, M) - F = CUSOLVER.SparseCholesky(M) - CUSOLVER.spcholesky_factorise(F, M, 1.e-12) - return F -end - -function run_covo(backend) - println("Starting run_covo...") - - # Build mesh - meshParam = (nx = nx, ny = ny, lx = 2l, ly = 2l, xc = 0.0, yc = 0.0) - tmp_path = "tmp.msh" - if get(ENV, "BenchmarkMode", "false") == "false" #hide - BcubeGmsh.gen_rectangle_mesh(tmp_path, :quad; meshParam...) - else #hide - if get(ENV, "MeshConfig", "quad") == "triquad" #hide - BcubeGmsh.gen_rectangle_mesh_with_tri_and_quad(tmp_path; meshParam...) #hide - else #hide - BcubeGmsh.gen_rectangle_mesh(tmp_path, :quad; meshParam...) #hide - end #hide - end #hide - - # Define variables and test functions - fs = FunctionSpace(fspace, degree) - mesh, U, V, mesh_cpu = setup(backend, tmp_path, fs) - u = FEFunction(U, KernelAbstractions.zeros(backend, Float64, get_ndofs(U))) - - @show Bcube.get_ndofs(U) - - # Define measures for cell and interior face integrations - dΩ = Measure(CellDomain(mesh), degquad) - dΓ = Measure(InteriorFaceDomain(mesh), degquad) - - # Declare periodic boundary conditions and - # create associated domains and measures - periodicBCType_x = PeriodicBCType(Translation(SA[-2l, 0.0]), ("East",), ("West",)) - periodicBCType_y = PeriodicBCType(Translation(SA[0.0, 2l]), ("South",), ("North",)) - Γ_perio_x_cpu = BoundaryFaceDomain(mesh_cpu, periodicBCType_x) - Γ_perio_y_cpu = BoundaryFaceDomain(mesh_cpu, periodicBCType_y) - Γ_perio_x = adapt(backend, Γ_perio_x_cpu) - Γ_perio_y = adapt(backend, Γ_perio_y_cpu) - dΓ_perio_x = Measure(Γ_perio_x, degquad) - dΓ_perio_y = Measure(Γ_perio_y, degquad) - - params = (dΩ = dΩ, dΓ = dΓ, dΓ_perio_x = dΓ_perio_x, dΓ_perio_y = dΓ_perio_y) - - # Init vtk - mkpath(joinpath(@__DIR__, outputpath)) - vtk = VtkHandler(output, mesh_cpu) - - # Init solution - t = 0.0 - - covo!(u, dΩ) - - # cache mass matrices - cache = (mass = _factorize(backend, Bcube.build_mass_matrix(U, V, dΩ)),) - - if get(ENV, "BenchmarkMode", "false") == "true" #hide - return u, U, V, params, cache - end - - # Write initial solution - WITH_VTK_OUTPUT && append_vtk(vtk, map(x -> adapt(get_backend(zeros(1)), x), (u...,)), t, params) - - # Time loop - for i in 1:nite - println("") - println("") - println("Iteration ", i, " / ", nite) - - ## Step forward in time - rhs(u, t) = compute_residual(u, V, params, cache) - if timeScheme == :ForwardEuler - unew = forward_euler(u, rhs, t, Δt) - elseif timeScheme == :RK3 - unew = rk3_ssp(u, rhs, t, Δt) - else - error("Unknown time scheme: $timeScheme") - end - - set_dof_values!(u, unew) - - t += Δt - - # Write solution to file - if (i % Int(max(floor(nite / nout), 1)) == 0) && WITH_VTK_OUTPUT - println("--> VTK export") - append_vtk(vtk, map(x -> adapt(get_backend(zeros(1)), x), (u...,)), t, params) - end - end - - # Summary and benchmark # ndofs total = 20480 - _rhs(u, t) = compute_residual(u, V, params, cache) - @btime begin - forward_euler($u, $_rhs, $time, $Δt) # 5.639 ms (1574 allocations: 2.08 MiB) - KernelAbstractions.synchronize(get_backend($backend)) - end - println("ndofs total = ", Bcube.get_ndofs(U)) - @show Δt, U₀, U₀ * t - @show boundary_names(mesh) - return nothing -end - -#const backendDevice = get_backend(ones(2)) ## CPU -const backendDevice = get_backend(CUDA.ones(2)) ## GPU -const backend = BcubeAccelerated.BcubeBackendAcc( - backendDevice; - kernel = BcubeAccelerated.KernelKA(), - sync = true, -) - -#const backend = Bcube.get_bcube_backend() -@show backend - -if get(ENV, "BenchmarkMode", "false") == "false" - mkpath(outputpath) - run_covo(backend) -end - -end #hide - -# CPU serial -# 330.144 ms (52 allocations: 31.00 MiB) -# ndofs total = 327680 - -# CPU (1threads) nx=128 -# 336.991 ms (397 allocations: 31.03 MiB) -# ndofs total = 327680 - -# CPU (80threads) nx=128 -# 21.236 ms (2741 allocations: 32.05 MiB) -# ndofs total = 327680 - -# GPU nx=128 -# 8.990 ms (13675 allocations: 1.22 MiB) -# ndofs total = 327680 \ No newline at end of file +module CovoGpu #hide +println("Running covo example...") #hide + +const dir = string(@__DIR__, "/") +using Bcube +using BcubeAccelerated +using KernelAbstractions +using CUDA, CUDA.CUSPARSE, CUDA.CUSOLVER +using BcubeGmsh +using Adapt +using BcubeVTK +using LinearAlgebra +using StaticArrays +using StaticArrays +using BenchmarkTools +using UnPack + +const WITH_VTK_OUTPUT = false + +function compute_residual(_u, V, params, cache) + u = get_fe_functions(_u) + + # alias on measures + @unpack dΩ, dΓ, dΓ_perio_x, dΓ_perio_y = params + + # face normals for each face domain (lazy, no computation at this step) + nΓ = get_face_normals(dΓ) + nΓ_perio_x = get_face_normals(dΓ_perio_x) + nΓ_perio_y = get_face_normals(dΓ_perio_y) + + # flux residuals from faces for all variables + function l(v) + ∫(flux_Ω(u, v))dΩ + + -∫(flux_Γ(u, v, nΓ))dΓ + + -∫(flux_Γ(u, v, nΓ_perio_x))dΓ_perio_x + + -∫(flux_Γ(u, v, nΓ_perio_y))dΓ_perio_y + end + # l1(v) = ∫(flux_Ω(u, v))dΩ + # l2(v) = ∫(flux_Γ(u, v, nΓ))dΓ + # l3(v) = ∫(flux_Γ(u, v, nΓ_perio_x))dΓ_perio_x + # l4(v) = ∫(flux_Γ(u, v, nΓ_perio_y))dΓ_perio_y + # @show "l1" + # assemble_linear(l1, V) + # @show "l2" + # assemble_linear(l2, V) + # @show "l3" + # assemble_linear(l3, V) + # @show "l4" + # assemble_linear(l4, V) + # error("ici") + _rhs = assemble_linear(l, V) + rhs = similar(_rhs) + Bcube.__solve!(rhs, cache.mass, _rhs, Bcube.get_bcube_backend(dΩ)) + return rhs +end + +""" + flux_Ω(u, v) + +Compute volume residual using the lazy-operators approach +""" +flux_Ω(u, v) = _flux_Ω ∘ cellvar(u, v) +cellvar(u, v) = (u, map(∇, v)) +function _flux_Ω(u, ∇v) + ρ, ρu, ρE, ϕ = u + ∇λ_ρ, ∇λ_ρu, ∇λ_ρE, ∇λ_ϕ = ∇v + + vel = ρu ./ ρ + ρuu = ρu * transpose(vel) + p = pressure(ρ, ρu, ρE, γ) + + flux_ρ = ρu + flux_ρu = ρuu + p * I + flux_ρE = (ρE + p) .* vel + flux_ϕ = ϕ .* vel + + return ∇λ_ρ ⋅ flux_ρ + ∇λ_ρu ⊡ flux_ρu + ∇λ_ρE ⋅ flux_ρE + ∇λ_ϕ ⋅ flux_ϕ +end + +""" + flux_Γ(u,v,n) + +Flux at the interface is defined by a composition of two functions: +* facevar(u,v,n) defines the input states which are needed for + the riemann flux using operator notations +* flux_roe(w) defines the Riemann flux (as usual) +""" +flux_Γ(u, v, n) = flux_roe ∘ (side⁻(u), side⁺(u), jump(v), side⁻(n)) + +""" + flux_roe(w) +""" +function flux_roe(ui, uj, δv, nij) + # destructuring inputs given by `facevar` function + + nx, ny = nij + ρ1, ρu1, ρE1, ϕ1 = ui + ρ2, ρu2, ρE2, ϕ2 = uj + δλ_ρ1, δλ_ρu1, δλ_ρE1, δλ_ϕ1 = δv + ρux1, ρuy1 = ρu1 + ρux2, ρuy2 = ρu2 + + # Closure + u1 = ρux1 / ρ1 + v1 = ρuy1 / ρ1 + u2 = ρux2 / ρ2 + v2 = ρuy2 / ρ2 + p1 = pressure(ρ1, ρu1, ρE1, γ) + p2 = pressure(ρ2, ρu2, ρE2, γ) + + H2 = (γ / (γ - 1)) * p2 / ρ2 + (u2 * u2 + v2 * v2) / 2.0 + H1 = (γ / (γ - 1)) * p1 / ρ1 + (u1 * u1 + v1 * v1) / 2.0 + + R = √(ρ1 / ρ2) + invR1 = 1.0 / (R + 1) + uAv = (R * u1 + u2) * invR1 + vAv = (R * v1 + v2) * invR1 + Hav = (R * H1 + H2) * invR1 + cAv = √(abs((γ - 1) * (Hav - (uAv * uAv + vAv * vAv) / 2.0))) + ecAv = (uAv * uAv + vAv * vAv) / 2.0 + + λ1 = nx * uAv + ny * vAv + λ3 = λ1 + cAv + λ4 = λ1 - cAv + + d1 = ρ1 - ρ2 + d2 = ρ1 * u1 - ρ2 * u2 + d3 = ρ1 * v1 - ρ2 * v2 + d4 = ρE1 - ρE2 + + # computation of the centered part of the flux + flu11 = nx * ρ2 * u2 + ny * ρ2 * v2 + flu21 = nx * p2 + flu11 * u2 + flu31 = ny * p2 + flu11 * v2 + flu41 = H2 * flu11 + + # Temp variables + rc1 = (γ - 1) / cAv + rc2 = (γ - 1) / cAv / cAv + uq41 = ecAv / cAv + cAv / (γ - 1) + uq42 = nx * uAv + ny * vAv + + fdc1 = max(λ1, 0.0) * (d1 + rc2 * (-ecAv * d1 + uAv * d2 + vAv * d3 - d4)) + fdc2 = max(λ1, 0.0) * ((nx * vAv - ny * uAv) * d1 + ny * d2 - nx * d3) + fdc3 = + max(λ3, 0.0) * ( + (-uq42 * d1 + nx * d2 + ny * d3) / 2.0 + + rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0 + ) + fdc4 = + max(λ4, 0.0) * ( + (uq42 * d1 - nx * d2 - ny * d3) / 2.0 + + rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0 + ) + + duv1 = fdc1 + (fdc3 + fdc4) / cAv + duv2 = uAv * fdc1 + ny * fdc2 + (uAv / cAv + nx) * fdc3 + (uAv / cAv - nx) * fdc4 + duv3 = vAv * fdc1 - nx * fdc2 + (vAv / cAv + ny) * fdc3 + (vAv / cAv - ny) * fdc4 + duv4 = + ecAv * fdc1 + + (ny * uAv - nx * vAv) * fdc2 + + (uq41 + uq42) * fdc3 + + (uq41 - uq42) * fdc4 + + v₁₂ = 0.5 * ((u1 + u2) * nx + (v1 + v2) * ny) + fluxϕ = max(0.0, v₁₂) * ϕ1 + min(0.0, v₁₂) * ϕ2 + + return ( + δλ_ρ1 * (flu11 + duv1) + + δλ_ρu1 ⋅ (SA[flu21 + duv2, flu31 + duv3]) + + δλ_ρE1 * (flu41 + duv4) + + δλ_ϕ1 * (fluxϕ) + ) +end + +""" +Time integration of `f(q, t)` over a timestep `Δt`. +""" +forward_euler(q, f::Function, t, Δt) = get_dof_values(q) .+ Δt .* f(q, t) + +""" + rk3_ssp(q, f::Function, t, Δt) + +`f(q, t)` is the function to integrate. +""" +function rk3_ssp(q, f::Function, t, Δt) + stepper(q, t) = forward_euler(q, f, t, Δt) + _q0 = get_dof_values(q) + + _q1 = stepper(q, Δt) + + set_dof_values!(q, _q1) + _q2 = (3 / 4) .* _q0 .+ (1 / 4) .* stepper(q, t + Δt) + + set_dof_values!(q, _q2) + _q1 .= (1 / 3) * _q0 .+ (2 / 3) .* stepper(q, t + Δt / 2) + + return _q1 +end + +""" + pressure(ρ, ρu, ρE, γ) + +Computes pressure from perfect gaz law +""" +function pressure(ρ::Number, ρu::AbstractVector, ρE::Number, γ) + vel = ρu ./ ρ + ρuu = ρu * transpose(vel) + p = (γ - 1) * (ρE - tr(ρuu) / 2) + return p +end + +""" + Init field with a vortex (for the COVO test case) +""" +function covo!(q, dΩ) + + # Intermediate vars + Cₚ = γ * r / (γ - 1) + + r²(x) = ((x[1] .- xvc) .^ 2 + (x[2] .- yvc) .^ 2) ./ Rc^2 + # Temperature + T(x) = T₀ .- β^2 * U₀^2 / (2 * Cₚ) .* exp.(-r²(x)) + # Velocity + ux(x) = U₀ .- β * U₀ / Rc .* (x[2] .- yvc) .* exp.(-r²(x) ./ 2) + uy(x) = V₀ .+ β * U₀ / Rc .* (x[1] .- xvc) .* exp.(-r²(x) ./ 2) + # Density + ρ(x) = ρ₀ .* (T(x) ./ T₀) .^ (1.0 / (γ - 1)) + # momentum + ρu(x) = SA[ρ(x) * ux(x), ρ(x) * uy(x)] + # Energy + ρE(x) = ρ(x) * ((Cₚ / γ) .* T(x) + (ux(x) .^ 2 + uy(x) .^ 2) ./ 2) + # Passive scalar + ϕ(x) = Rc^2 * r²(x) < 0.01 ? exp(-r²(x) ./ 2) : 0.0 + + f = map(PhysicalFunction, (ρ, ρu, ρE, ϕ)) + projection_l2!(q, f, dΩ) + return nothing +end + +""" + Tiny struct to ease the VTK output +""" +mutable struct VtkHandler + basename::Any + mesh::Any + ite::Any + VtkHandler(basename, mesh) = new(basename, mesh, 0) +end + +""" + Write solution to vtk + Wrapper for `write_vtk` +""" +function append_vtk(vtk, vars, t, params) + ρ, ρu, ρE, ϕ = vars + + mesh = vtk.mesh + mesh_degree = 1 + vtk_degree = maximum(x -> get_degree(Bcube.get_function_space(get_fespace(x))), vars) + vtk_degree = max(1, mesh_degree, vtk_degree) + + # p = pressure(ρ, ρu, ρE, γ) + p = pressure ∘ (ρ, ρu, ρE, γ) + dict_vars_dg = Dict("rho" => ρ, "rhou" => ρu, "rhoE" => ρE, "phi" => ϕ, "p" => p) + write_file( + vtk.basename * "_DG.pvd", + mesh, + dict_vars_dg, + vtk.ite, + t; + discontinuous = true, + mesh_degree = vtk_degree, + collection_append = vtk.ite > 0, + ) + + # Update counter + vtk.ite += 1 +end + +function Adapt.adapt_structure(to, mesh::BcubeGmsh.GmshMetaData) + nothing +end +function Bcube.__solve!( + x, + A, + b, + backend::BcubeAccelerated.BcubeBackendAcc{K, CUDA.CUDABackend}, +) where {K} + F = CUSOLVER.SparseCholesky(A) + CUSOLVER.spcholesky_factorise(F, A, 1.e-12) + CUSOLVER.spcholesky_solve(F, b, x) + return nothing +end +function Bcube.__solve!( + x, + A::CUSOLVER.SparseCholesky, + b, + backend::BcubeAccelerated.BcubeBackendAcc{K, CUDA.CUDABackend}, +) where {K} + CUSOLVER.spcholesky_solve(A, b, x) + return nothing +end +# function Bcube.get_bcube_backend(a::CuArray) +# BcubeAccelerated.BcubeBackendAcc(KernelAbstractions.get_backend(a)) +# end + +# Settings +if get(ENV, "BenchmarkMode", "false") == "false" #hide + const cellfactor = 1 + const nx = 32 * cellfactor + 1 + const ny = 32 * cellfactor + 1 + const fspace = :Lagrange + const timeScheme = :ForwardEuler +else #hide + const nx = 128 + 1 #hide + const ny = 128 + 1 #hide + const fspace = :Lagrange + const timeScheme = :ForwardEuler +end #hide +const nperiod = 1 # number of turn +const CFL = 0.1 +const degree = 1 # FunctionSpace degree +const degquad = 2 * degree + 1 +const γ = 1.4 +const β = 0.2 # vortex intensity +const r = 287.15 # Perfect gaz constant +const T₀ = 300 # mean-flow temperature +const P₀ = 1e5 # mean-flow pressure +const M₀ = 0.5 # mean-flow mach number +const ρ₀ = 1.0 # mean-flow density +const xvc = 0.0 # x-center of vortex +const yvc = 0.0 # y-center of vortex +const Rc = 0.005 # Charasteristic vortex radius +const c₀ = √(γ * r * T₀) # Sound velocity +const U₀ = M₀ * c₀ # mean-flow velocity +const V₀ = 0.0 # mean-flow velocity +const ϕ₀ = 1.0 +const l = 0.05 # half-width of the domain +const Δt = CFL * 2 * l / (nx - 1) / ((1 + β) * U₀ + c₀) / (2 * degree + 1) +#const Δt = 5.e-7 +const nout = 100 # Number of time steps to save +const outputpath = "./covo/" +const output = joinpath(@__DIR__, outputpath, "covo_deg$degree") +const nite = Int(floor(nperiod * 2 * l / (U₀ * Δt))) + 1 + +function build_mass_matrix_test(U, V, dΩ::D) where {D <: AbstractMeasure} + _m(u, v) = ∫(u ⋅ v)dΩ + # @code_warntype _m(1,2) + @show typeof(_m(1, 2)) + BcubeAccelerated.test_arg_AK_foreachindex(get_backend(Bcube.get_bcube_backend(dΩ)), dΩ) + BcubeAccelerated.test_arg_AK_foreachindex( + get_backend(Bcube.get_bcube_backend(dΩ)), + _m(1, 2), + ) + return assemble_bilinear(_m, U, V) +end +function setup(backend, tmp_path, fs) + mesh_cpu = read_mesh(tmp_path) + rm(tmp_path) + U_sca_cpu = TrialFESpace(fs, mesh_cpu, :discontinuous; size = 1) # DG, scalar + U_vec_cpu = TrialFESpace(fs, mesh_cpu, :discontinuous; size = 2) # DG, vectoriel + U_cpu = MultiFESpace(U_sca_cpu, U_vec_cpu, U_sca_cpu, U_sca_cpu) + V_sca_cpu = TestFESpace(U_sca_cpu) + V_vec_cpu = TestFESpace(U_vec_cpu) + V_cpu = MultiFESpace(V_sca_cpu, V_vec_cpu, V_sca_cpu, V_sca_cpu) + + mesh = adapt(backend, mesh_cpu) + BcubeAccelerated.test_arg(backend, mesh) + U = adapt(backend, U_cpu) + V = adapt(backend, V_cpu) + BcubeAccelerated.test_arg(backend, U) + BcubeAccelerated.test_arg(backend, V) + return mesh, U, V, mesh_cpu +end + +_factorize(backend::Bcube.AbstractBcubeBackend, M) = __factorize(get_backend(backend), M) +__factorize(backend, M) = factorize(M) +function __factorize(backend::CUDA.CUDABackend, M) + F = CUSOLVER.SparseCholesky(M) + CUSOLVER.spcholesky_factorise(F, M, 1.e-12) + return F +end + +function run_covo(backend) + println("Starting run_covo...") + + # Build mesh + meshParam = (nx = nx, ny = ny, lx = 2l, ly = 2l, xc = 0.0, yc = 0.0) + tmp_path = "tmp.msh" + if get(ENV, "BenchmarkMode", "false") == "false" #hide + BcubeGmsh.gen_rectangle_mesh(tmp_path, :quad; meshParam...) + else #hide + if get(ENV, "MeshConfig", "quad") == "triquad" #hide + BcubeGmsh.gen_rectangle_mesh_with_tri_and_quad(tmp_path; meshParam...) #hide + else #hide + BcubeGmsh.gen_rectangle_mesh(tmp_path, :quad; meshParam...) #hide + end #hide + end #hide + + # Define variables and test functions + fs = FunctionSpace(fspace, degree) + mesh, U, V, mesh_cpu = setup(backend, tmp_path, fs) + u = FEFunction(U, KernelAbstractions.zeros(backend, Float64, get_ndofs(U))) + + @show Bcube.get_ndofs(U) + + # Define measures for cell and interior face integrations + dΩ = Measure(CellDomain(mesh), degquad) + dΓ = Measure(InteriorFaceDomain(mesh), degquad) + + # Declare periodic boundary conditions and + # create associated domains and measures + periodicBCType_x = PeriodicBCType(Translation(SA[-2l, 0.0]), ("East",), ("West",)) + periodicBCType_y = PeriodicBCType(Translation(SA[0.0, 2l]), ("South",), ("North",)) + Γ_perio_x_cpu = BoundaryFaceDomain(mesh_cpu, periodicBCType_x) + Γ_perio_y_cpu = BoundaryFaceDomain(mesh_cpu, periodicBCType_y) + Γ_perio_x = adapt(backend, Γ_perio_x_cpu) + Γ_perio_y = adapt(backend, Γ_perio_y_cpu) + dΓ_perio_x = Measure(Γ_perio_x, degquad) + dΓ_perio_y = Measure(Γ_perio_y, degquad) + + params = (dΩ = dΩ, dΓ = dΓ, dΓ_perio_x = dΓ_perio_x, dΓ_perio_y = dΓ_perio_y) + + # Init vtk + mkpath(joinpath(@__DIR__, outputpath)) + vtk = VtkHandler(output, mesh_cpu) + + # Init solution + t = 0.0 + + covo!(u, dΩ) + + # cache mass matrices + cache = (mass = _factorize(backend, Bcube.build_mass_matrix(U, V, dΩ)),) + + if get(ENV, "BenchmarkMode", "false") == "true" #hide + return u, U, V, params, cache + end + + # Write initial solution + WITH_VTK_OUTPUT && + append_vtk(vtk, map(x -> adapt(get_backend(zeros(1)), x), (u...,)), t, params) + + # Time loop + for i in 1:nite + println("") + println("") + println("Iteration ", i, " / ", nite) + + ## Step forward in time + rhs(u, t) = compute_residual(u, V, params, cache) + if timeScheme == :ForwardEuler + unew = forward_euler(u, rhs, t, Δt) + elseif timeScheme == :RK3 + unew = rk3_ssp(u, rhs, t, Δt) + else + error("Unknown time scheme: $timeScheme") + end + + set_dof_values!(u, unew) + + t += Δt + + # Write solution to file + if (i % Int(max(floor(nite / nout), 1)) == 0) && WITH_VTK_OUTPUT + println("--> VTK export") + append_vtk(vtk, map(x -> adapt(get_backend(zeros(1)), x), (u...,)), t, params) + end + end + + # Summary and benchmark # ndofs total = 20480 + _rhs(u, t) = compute_residual(u, V, params, cache) + @btime begin + forward_euler($u, $_rhs, $time, $Δt) # 5.639 ms (1574 allocations: 2.08 MiB) + KernelAbstractions.synchronize(get_backend($backend)) + end + println("ndofs total = ", Bcube.get_ndofs(U)) + @show Δt, U₀, U₀ * t + @show boundary_names(mesh) + return nothing +end + +#const backendDevice = get_backend(ones(2)) ## CPU +const backendDevice = get_backend(CUDA.ones(2)) ## GPU +const backend = BcubeAccelerated.BcubeBackendAcc( + backendDevice; + kernel = BcubeAccelerated.KernelKA(), + sync = true, +) + +#const backend = Bcube.get_bcube_backend() +@show backend + +if get(ENV, "BenchmarkMode", "false") == "false" + mkpath(outputpath) + run_covo(backend) +end + +end #hide + +# CPU serial +# 330.144 ms (52 allocations: 31.00 MiB) +# ndofs total = 327680 + +# CPU (1threads) nx=128 +# 336.991 ms (397 allocations: 31.03 MiB) +# ndofs total = 327680 + +# CPU (80threads) nx=128 +# 21.236 ms (2741 allocations: 32.05 MiB) +# ndofs total = 327680 + +# GPU nx=128 +# 8.990 ms (13675 allocations: 1.22 MiB) +# ndofs total = 327680 diff --git a/test/linear_transport_gpu.jl b/test/linear_transport_gpu.jl index b1b8092..af46295 100644 --- a/test/linear_transport_gpu.jl +++ b/test/linear_transport_gpu.jl @@ -214,4 +214,4 @@ enable_timer!(to) main(nx, ny, nite, degree, backend) show(to) -end \ No newline at end of file +end diff --git a/test/runtests.jl b/test/runtests.jl index c626e0e..6313719 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,14 +8,23 @@ using CUDA, CUDA.CUSPARSE, CUDA.CUSOLVER using SparseArrays, LinearAlgebra using TimerOutputs +WITH_CPU_TESTS = get(ENV, "WITH_CPU_TESTS", "true") == "true" +WITH_CUDA_TESTS = get(ENV, "WITH_CUDA_TESTS", "true") == "true" + +@show WITH_CPU_TESTS, WITH_CUDA_TESTS + @testset "BcubeAccelerated.jl" begin - include("test_domain.jl") + if WITH_CPU_TESTS + include("test_domain.jl") + end - include("linear_transport_gpu.jl") - using .LinearTransportGpu + if WITH_CUDA_TESTS + include("linear_transport_gpu.jl") + using .LinearTransportGpu - include("covo_gpu.jl") - using .CovoGpu + include("covo_gpu.jl") + using .CovoGpu + end @test true # Checkpoint end diff --git a/test/test_domain.jl b/test/test_domain.jl index 7972d5f..7b630b9 100644 --- a/test/test_domain.jl +++ b/test/test_domain.jl @@ -14,4 +14,4 @@ @test Bcube.get_return_type_and_codim(f1, Ω, backend) == (Float64, (1,)) f2 = PhysicalFunction(x -> x) @test Bcube.get_return_type_and_codim(f2, Ω, backend) == (Float64, (2,)) -end \ No newline at end of file +end