From 74a090f25f77e796b6f471f753ebf36ef6587e61 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 30 Jun 2023 09:48:36 +0200 Subject: [PATCH 1/6] Add Euler freestream channel example --- .../elixir_euler_channel_freestream.jl | 90 +++++++++++++++++++ 1 file changed, 90 insertions(+) create mode 100644 channel-transport/elixir_euler_channel_freestream.jl diff --git a/channel-transport/elixir_euler_channel_freestream.jl b/channel-transport/elixir_euler_channel_freestream.jl new file mode 100644 index 000000000..791052666 --- /dev/null +++ b/channel-transport/elixir_euler_channel_freestream.jl @@ -0,0 +1,90 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_constant + +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D +# function mapping(xi_, eta_) +# # Transform input variables between -1 and 1 onto [0,3] +# xi = 1.5 * xi_ + 1.5 +# eta = 1.5 * eta_ + 1.5 +# +# y = eta + 3/8 * (cos(1.5 * pi * (2 * xi - 3)/3) * +# cos(0.5 * pi * (2 * eta - 3)/3)) +# +# x = xi + 3/8 * (cos(0.5 * pi * (2 * xi - 3)/3) * +# cos(2 * pi * (2 * y - 3)/3)) +# +# return SVector(x, y) +# end + +############################################################################### +# Get the uncurved mesh from a file (downloads the file if not available locally) + +# Unstructured mesh with 48 cells of the square domain [-1, 1]^n +# mesh_file = joinpath(@__DIR__, "square_unstructured_1.inp") +# isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", +# mesh_file) +# +# # Map the unstructured mesh with the mapping above +# mesh = P4estMesh{2}(mesh_file, polydeg=3, mapping=mapping, initial_refinement_level=1) + +coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y)) +coordinates_max = (6.0, 2.0) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh +trees_per_dimension = (6, 2) +mesh = P4estMesh(trees_per_dimension, + polydeg=1, initial_refinement_level=1, + coordinates_min=coordinates_min, coordinates_max=coordinates_max, + periodicity=(true, true)) + # periodicity=(false, false)) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + # boundary_conditions=Dict( + # :all => BoundaryConditionDirichlet(initial_condition) + # ) + ) + + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +# save_solution = SaveSolutionCallback(interval=100, +# save_initial_solution=true, +# save_final_solution=true, +# solution_variables=cons2prim) + +stepsize_callback = StepsizeCallback(cfl=1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + # save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); +summary_callback() # print the timer summary From 10657ff0008656bd886f887bc6be9fb71acdc778 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 30 Jun 2023 09:58:32 +0200 Subject: [PATCH 2/6] Add horizontal freestream --- .../elixir_euler_channel_freestream.jl | 35 ++++++------------- 1 file changed, 10 insertions(+), 25 deletions(-) diff --git a/channel-transport/elixir_euler_channel_freestream.jl b/channel-transport/elixir_euler_channel_freestream.jl index 791052666..af755fa23 100644 --- a/channel-transport/elixir_euler_channel_freestream.jl +++ b/channel-transport/elixir_euler_channel_freestream.jl @@ -8,35 +8,20 @@ using Trixi equations = CompressibleEulerEquations2D(1.4) -initial_condition = initial_condition_constant +function initial_condition_channel(x, t, equations::CompressibleEulerEquations2D) + rho = 1.0 + rho_v1 = 0.2 + rho_v2 = 0.0 + rho_e = 10.0 + return SVector(rho, rho_v1, rho_v2, rho_e) +end -solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) +initial_condition = initial_condition_channel -# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D -# function mapping(xi_, eta_) -# # Transform input variables between -1 and 1 onto [0,3] -# xi = 1.5 * xi_ + 1.5 -# eta = 1.5 * eta_ + 1.5 -# -# y = eta + 3/8 * (cos(1.5 * pi * (2 * xi - 3)/3) * -# cos(0.5 * pi * (2 * eta - 3)/3)) -# -# x = xi + 3/8 * (cos(0.5 * pi * (2 * xi - 3)/3) * -# cos(2 * pi * (2 * y - 3)/3)) -# -# return SVector(x, y) -# end +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) ############################################################################### -# Get the uncurved mesh from a file (downloads the file if not available locally) - -# Unstructured mesh with 48 cells of the square domain [-1, 1]^n -# mesh_file = joinpath(@__DIR__, "square_unstructured_1.inp") -# isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", -# mesh_file) -# -# # Map the unstructured mesh with the mapping above -# mesh = P4estMesh{2}(mesh_file, polydeg=3, mapping=mapping, initial_refinement_level=1) +# Create mesh coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y)) coordinates_max = (6.0, 2.0) # maximum coordinates (max(x), max(y)) From b8b58b92827d37c1d964e7d2f92b1a5e616a4ce7 Mon Sep 17 00:00:00 2001 From: Benjamin Uekermann Date: Fri, 30 Jun 2023 11:01:39 +0200 Subject: [PATCH 3/6] Add preCICE calls to Trixi elixir --- .../fluid.jl} | 87 ++++++++++++++++++- 1 file changed, 84 insertions(+), 3 deletions(-) rename channel-transport/{elixir_euler_channel_freestream.jl => fluid-trixi/fluid.jl} (67%) diff --git a/channel-transport/elixir_euler_channel_freestream.jl b/channel-transport/fluid-trixi/fluid.jl similarity index 67% rename from channel-transport/elixir_euler_channel_freestream.jl rename to channel-transport/fluid-trixi/fluid.jl index af755fa23..58e8de00b 100644 --- a/channel-transport/elixir_euler_channel_freestream.jl +++ b/channel-transport/fluid-trixi/fluid.jl @@ -2,6 +2,7 @@ using Downloads: download using OrdinaryDiffEq using Trixi +using PreCICE ############################################################################### # semidiscretization of the compressible Euler equations @@ -10,15 +11,15 @@ equations = CompressibleEulerEquations2D(1.4) function initial_condition_channel(x, t, equations::CompressibleEulerEquations2D) rho = 1.0 - rho_v1 = 0.2 - rho_v2 = 0.0 + rho_v1 = 1.0 + rho_v2 = 1.0 rho_e = 10.0 return SVector(rho, rho_v1, rho_v2, rho_e) end initial_condition = initial_condition_channel -solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) +solver = DGSEM(polydeg=2, surface_flux=flux_lax_friedrichs) ############################################################################### # Create mesh @@ -73,3 +74,83 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks); summary_callback() # print the timer summary + + +################################################################################## + + +commRank = 0 +commSize = 1 + + +PreCICE.createSolverInterface("Fluid", "../precice-config.xml", commRank, commSize) + +dimensions = PreCICE.getDimensions() + + +meshID = PreCICE.getMeshID("Fluid-Mesh") + +dataID = PreCICE.getDataID("Velocity", meshID) + + +numberOfVertices = nelements(solver, semi.cache) + +vertices = transpose(semi.cache.elements.node_coordinates[1:2,2,2,1:numberOfVertices]) + +vertexIDs = PreCICE.setMeshVertices(meshID, vertices) + + +writeData = zeros(numberOfVertices, dimensions) + +u = Trixi.wrap_array(sol.u[end], semi) + +let # setting local scope for dt outside of the while loop + + dt = PreCICE.initialize() + + while PreCICE.isCouplingOngoing() + + for i = 1:numberOfVertices + node_vars = Trixi.get_node_vars(u, equations, solver, 2, 2, i) + _, v1, v2, _ = cons2prim(node_vars, equations) + writeData[i, 1] = v1 + writeData[i, 2] = v2 + end + + PreCICE.writeBlockVectorData(dataID, vertexIDs, writeData) + + dt = PreCICE.advance(dt) + + end # while + +end # let + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From a20f27af13d59760e36868f2daa3d41502a9809b Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 30 Jun 2023 11:09:56 +0200 Subject: [PATCH 4/6] Use external time stepping loop --- channel-transport/fluid-trixi/fluid.jl | 58 +++++++++----------------- 1 file changed, 20 insertions(+), 38 deletions(-) diff --git a/channel-transport/fluid-trixi/fluid.jl b/channel-transport/fluid-trixi/fluid.jl index 58e8de00b..0148de743 100644 --- a/channel-transport/fluid-trixi/fluid.jl +++ b/channel-transport/fluid-trixi/fluid.jl @@ -63,18 +63,22 @@ alive_callback = AliveCallback(analysis_interval=analysis_interval) stepsize_callback = StepsizeCallback(cfl=1.0) callbacks = CallbackSet(summary_callback, - analysis_callback, alive_callback, + analysis_callback, alive_callback,) # save_solution, - stepsize_callback) + # stepsize_callback) ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); -summary_callback() # print the timer summary +# sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), +# dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback +# save_everystep=false, callback=callbacks); +# summary_callback() # print the timer summary + +integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); ################################################################################## @@ -102,14 +106,20 @@ vertexIDs = PreCICE.setMeshVertices(meshID, vertices) writeData = zeros(numberOfVertices, dimensions) -u = Trixi.wrap_array(sol.u[end], semi) + +solver_dt = 0.005 + +u = Trixi.wrap_array(integrator.u, semi) let # setting local scope for dt outside of the while loop - dt = PreCICE.initialize() + precice_dt = PreCICE.initialize() while PreCICE.isCouplingOngoing() + dt = min(solver_dt, precice_dt) + step!(integrator, dt, true) + for i = 1:numberOfVertices node_vars = Trixi.get_node_vars(u, equations, solver, 2, 2, i) _, v1, v2, _ = cons2prim(node_vars, equations) @@ -119,38 +129,10 @@ let # setting local scope for dt outside of the while loop PreCICE.writeBlockVectorData(dataID, vertexIDs, writeData) - dt = PreCICE.advance(dt) + precice_dt = PreCICE.advance(dt) end # while end # let - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +summary_callback() # print the timer summary From 4aff57d455b5050fd98c9db67bcd0160cbf12ed3 Mon Sep 17 00:00:00 2001 From: Benjamin Uekermann Date: Fri, 30 Jun 2023 11:37:43 +0200 Subject: [PATCH 5/6] Extend README --- channel-transport/README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/channel-transport/README.md b/channel-transport/README.md index 7a7d8e5ff..9ad8f942b 100644 --- a/channel-transport/README.md +++ b/channel-transport/README.md @@ -28,6 +28,8 @@ Fluid participant: * Nutils. For more information, have a look at the [Nutils adapter documentation](https://www.precice.org/adapter-nutils.html). This Nutils solver requires at least Nutils v7.0. +* Trixi. Required dependencies: [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) (tested with v0.5.30), [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) (tested with v6.53.2), [PreCICE.jl](https://github.com/precice/PreCICE.jl) (tested with v2.5.0). Currently, Trixi uses Euler equations (not incompressible Navier-Stokes) and periodic boundary conditions in all directions. + Transport participant: * Nutils. For more information, have a look at the [Nutils adapter documentation](https://www.precice.org/adapter-nutils.html). This Nutils solver requires at least Nutils v7.0. From d087b71248023d62aff4ed6a90e04f6616a03c46 Mon Sep 17 00:00:00 2001 From: Benjamin Uekermann Date: Fri, 30 Jun 2023 11:38:07 +0200 Subject: [PATCH 6/6] Change initial velocity to more channel-like setup --- channel-transport/fluid-trixi/fluid.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/channel-transport/fluid-trixi/fluid.jl b/channel-transport/fluid-trixi/fluid.jl index 0148de743..bfe58a26a 100644 --- a/channel-transport/fluid-trixi/fluid.jl +++ b/channel-transport/fluid-trixi/fluid.jl @@ -12,7 +12,7 @@ equations = CompressibleEulerEquations2D(1.4) function initial_condition_channel(x, t, equations::CompressibleEulerEquations2D) rho = 1.0 rho_v1 = 1.0 - rho_v2 = 1.0 + rho_v2 = 0.0 rho_e = 10.0 return SVector(rho, rho_v1, rho_v2, rho_e) end