General plan and layout

[ranocha]

In GitLab by @ranocha on May 7, 2020, 17:07

Hi! I'm totally new to this project and haven't been involved in the discussions and the work resulting in the current version. Hence, I hesitate a bit to start discussions that can possibly necessitate huge changes to your workflow and/or the code. On the other hand, I couldn't find such discussions online (since they probably happened in person) and the resulting comments may be of interest for future users of Trixi.

My first impression of Trixi is that it's heavily inspired by monolithic Fortran codes like Flexi/Fluxo, which is totally natural, of course. In my opinion, a different approach might be useful in Julia: Instead of providing a single executable running a parameter file, an approach to create a library of tools might be of interest.
For example, if a student shall run some experiments with a slightly different initial condition, Trixi itself has to be modified, e.g. at https://gitlab.mi.uni-koeln.de/numsim/code/Trixi.jl/-/blob/master/src/equations/linearscalaradvection.jl#L46.
An alternative to parameter files could be to supply the parameters as pure Julia code using the Trixi library. In that case, students could write the initial condition as a Julia function without the need to modify Trixi itself. Since functions are first-class citizens in Julia, they can be passed around as arguments (or as part of structs).

The structure of Trixi using many modules and global variables reminds me of Fortran modules etc. While using constant global variables is okay (type stable) in Julia, passing parameters is more encouraged.

Another reminder of Fortran and co are the type annotations for function arguments. Unless they direct multiple dispatch, they shouldn't be used in Julia.

function calc_surface_flux!(surface_flux, neighbor_ids,
                            u_surfaces, dg::D, ::Val{true},
                            orientations)
   ...
end

is compiled to exactly the same machine code as

function calc_surface_flux!(surface_flux::Array{Float64, 4}, neighbor_ids::Matrix{Int},
                            u_surfaces::Array{Float64, 4}, dg::Dg, ::Val{true},
                            orientations::Vector{Int})
   ...
end

if the arguments are the same. The first version is more general, since the array type can be exchanged. For example, it is possible to use the same code to run a simulation with Float32 (or possibly some floating point type with higher accuracy), which can be useful sometimes. I prefer to include the additional kind of documentation provided by the type annotations in docstrings for the functions.

The multiple dispatch feature of Julia seems to be abused sometimes in the current version, e.g. in https://gitlab.mi.uni-koeln.de/numsim/code/Trixi.jl/-/blob/master/src/solvers/dg.jl#L1171. Such type-unstable code results in a possibly huge performance impact.
For example, I started from the current master branch with this git diff.

diff --git a/parameters_ec.toml b/parameters_ec.toml
index 8cd35d1..3bf5612 100644
--- a/parameters_ec.toml
+++ b/parameters_ec.toml
@@ -18,7 +18,7 @@ volume_flux_type = "chandrashekar_ec"
 # volume_flux_type = "central"
 # sources = 
 t_start = 0.0
-t_end   = 0.4
+t_end   = 4.0 # 0.4
 
 # restart = true
 # restart_filename = "out/restart_000100.h5"
diff --git a/src/equations/euler.jl b/src/equations/euler.jl
index c3978b9..5fdc582 100644
--- a/src/equations/euler.jl
+++ b/src/equations/euler.jl
@@ -670,62 +670,67 @@ function Equations.riemann!(surface_flux::AbstractArray{Float64, 1},
   calcflux1D!(f_ll, equation, rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll, orientation)
   calcflux1D!(f_rr, equation, rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr, orientation)
 
-  if equation.surface_flux_type == :laxfriedrichs
-    λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
-    surface_flux[1] = 1/2 * (f_ll[1] + f_rr[1]) - 1/2 * λ_max * (rho_rr    - rho_ll)
-    surface_flux[2] = 1/2 * (f_ll[2] + f_rr[2]) - 1/2 * λ_max * (rho_v1_rr - rho_v1_ll)
-    surface_flux[3] = 1/2 * (f_ll[3] + f_rr[3]) - 1/2 * λ_max * (rho_v2_rr - rho_v2_ll)
-    surface_flux[4] = 1/2 * (f_ll[4] + f_rr[4]) - 1/2 * λ_max * (rho_e_rr  - rho_e_ll)
-  elseif equation.surface_flux_type in (:central, :kennedygruber, :chandrashekar_ec, :yuichi)
-    symmetric_twopoint_flux!(surface_flux, Val(equation.surface_flux_type),
-                             equation, orientation,
-                             rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll,
-                             rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr)
-
-  elseif equation.surface_flux_type == :hllc
-    error("not yet implemented or tested")
-    v_tilde = (sqrt(rho_ll) * v_ll + sqrt(rho_rr) * v_rr) / (sqrt(rho_ll) + sqrt(rho_rr))
-    h_ll = (rho_e_ll + p_ll) / rho_ll
-    h_rr = (rho_e_rr + p_rr) / rho_rr
-    h_tilde = (sqrt(rho_ll) * h_ll + sqrt(rho_rr) * h_rr) / (sqrt(rho_ll) + sqrt(rho_rr))
-    c_tilde = sqrt((equation.gamma - 1) * (h_tilde - 1/2 * v_tilde^2))
-    s_ll = v_tilde - c_tilde
-    s_rr = v_tilde + c_tilde
-
-    if s_ll > 0
-      surface_flux[1, surface_id] = f_ll[1]
-      surface_flux[2, surface_id] = f_ll[2]
-      surface_flux[3, surface_id] = f_ll[3]
-    elseif s_rr < 0
-      surface_flux[1, surface_id] = f_rr[1]
-      surface_flux[2, surface_id] = f_rr[2]
-      surface_flux[3, surface_id] = f_rr[3]
-    else
-      s_star = ((p_rr - p_ll + rho_ll * v_ll * (s_ll - v_ll) - rho_rr * v_rr * (s_rr - v_rr))
-                / (rho_ll * (s_ll - v_ll) - rho_rr * (s_rr - v_rr)))
-      if s_ll <= 0 && 0 <= s_star
-        surface_flux[1, surface_id] = (f_ll[1] + s_ll *
-            (rho_ll * (s_ll - v_ll)/(s_ll - s_star) - rho_ll))
-        surface_flux[2, surface_id] = (f_ll[2] + s_ll *
-            (rho_ll * (s_ll - v_ll)/(s_ll - s_star) * s_star - rho_v_ll))
-        surface_flux[3, surface_id] = (f_ll[3] + s_ll *
-            (rho_ll * (s_ll - v_ll)/(s_ll - s_star) *
-            (rho_e_ll/rho_ll + (s_star - v_ll) * (s_star + rho_ll/(rho_ll * (s_ll - v_ll))))
-            - rho_e_ll))
-      else
-        surface_flux[1, surface_id] = (f_rr[1] + s_rr *
-            (rho_rr * (s_rr - v_rr)/(s_rr - s_star) - rho_rr))
-        surface_flux[2, surface_id] = (f_rr[2] + s_rr *
-            (rho_rr * (s_rr - v_rr)/(s_rr - s_star) * s_star - rho_v_rr))
-        surface_flux[3, surface_id] = (f_rr[3] + s_rr *
-            (rho_rr * (s_rr - v_rr)/(s_rr - s_star) *
-            (rho_e_rr/rho_rr + (s_star - v_rr) * (s_star + rho_rr/(rho_rr * (s_rr - v_rr))))
-            - rho_e_rr))
-      end
-    end
-  else
-    error("unknown Riemann solver '$(string(equation.surface_flux_type))'")
-  end
+  symmetric_twopoint_flux!(surface_flux, Val(:chandrashekar_ec),
+                           equation, orientation,
+                           rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll,
+                           rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr)
+
 end
 
 # Original riemann! implementation, non-optimized but easier to understand

If I understood the code correctly, the result should be the same since I use surface_flux_type = "chandrashekar_ec" in parameters_ec.toml. (Of course, this code is less general but it serves the purpose to demonstrate type stability in a simple way.)
Running the current master with this parameters_ec.toml, I get

              trixi                     Time                   Allocations      
                                ----------------------   -----------------------
        Tot / % measured:            18.1s / 99.5%            589MiB / 83.5%    

 Section                ncalls     time   %tot     avg     alloc   %tot      avg
 -------------------------------------------------------------------------------
 main loop                   1    18.0s   100%   18.0s    480MiB  97.5%   480MiB
   rhs                   1.75k    17.3s  95.8%  9.88ms   99.2MiB  20.2%  58.1KiB
     surface flux        1.75k    13.6s  75.5%  7.78ms   10.2MiB  2.06%  5.95KiB
     volume integral     1.75k    2.89s  16.0%  1.65ms   76.3MiB  15.5%  44.6KiB
     ...

With the type stable version from git diff, I get

        Tot / % measured:            5.22s / 90.0%            675MiB / 73.0%    

 Section                ncalls     time   %tot     avg     alloc   %tot      avg
 -------------------------------------------------------------------------------
 main loop                   1    4.70s   100%   4.70s    480MiB  97.5%   480MiB
   rhs                   1.75k    3.99s  84.8%  2.28ms   99.3MiB  20.2%  58.1KiB
     volume integral     1.75k    2.92s  62.2%  1.67ms   76.3MiB  15.5%  44.6KiB
     ...
     surface flux        1.75k    235ms  4.99%   134μs   10.2MiB  2.07%  5.97KiB

That impact is way more than I expected - I hope I didn't code bullshit here...
The additional git diff

diff --git a/src/equations/euler.jl b/src/equations/euler.jl
index c3978b9..b1f1b73 100644
--- a/src/equations/euler.jl
+++ b/src/equations/euler.jl
@@ -358,7 +358,9 @@ end
                                               equation::Euler,
                                               u::AbstractArray{Float64},
                                               element_id::Int, n_nodes::Int)
-  calcflux_twopoint!(f1, f2, f1_diag, f2_diag, Val(equation.volume_flux_type),
+  # calcflux_twopoint!(f1, f2, f1_diag, f2_diag, Val(equation.volume_flux_type),
+  #                    equation, u, element_id, n_nodes)
+  calcflux_twopoint!(f1, f2, f1_diag, f2_diag, Val(:chandrashekar_ec),
                      equation, u, element_id, n_nodes)
 end

improves the performance for the volume terms to

        Tot / % measured:            3.13s / 97.4%            576MiB / 83.2%    

 Section                ncalls     time   %tot     avg     alloc   %tot      avg
 -------------------------------------------------------------------------------
 main loop                   1    3.04s   100%   3.04s    467MiB  97.5%   467MiB
   rhs                   1.75k    2.25s  73.9%  1.29ms   85.6MiB  17.9%  50.1KiB
     volume integral     1.75k    1.20s  39.5%   687μs   62.6MiB  13.1%  36.6KiB
     ...
     surface flux        1.75k    237ms  7.79%   136μs   10.2MiB  2.13%  5.97KiB

A general solution to get type stable code in a library like approach could be to use parametric types and supply the fluxes as Julia functions instead of strings or symbols.

I saw the additional repositories Abraxas and trixi-tests but haven't looked at them yet. In Julia, it is pretty common to include some (unit) tests in every package which are triggered for every PR on GitHub (merge request on GitLab). In this way, some additional example code is provided and some basic functionality is guaranteed to work with new changes. While I didn't like that approach at the beginning (since I wasn't used to it), I really started to like it while developing more complex code, in particular if some time has passed between development cycles. I think it would be really nice if some small/basic tests are included here in test/runtests.jl etc. They could be triggered for every merge request as GitLab pipeline, I think. If Trixi shall be registered as Julia package (as hinted to in https://gitlab.mi.uni-koeln.de/numsim/code/Trixi.jl/-/issues/20#note_276), running tests and reporting coverage results (e.g. via Travis and Codecov/Coveralls) is usually seen as a requirement and sign of good coding standards.

Okay, this grew into a lot of text. I hope we can start a constructive discussion and improve Trixi, building on your nice work :)

[ranocha]

In GitLab by @gregorgassner on May 7, 2020, 17:13

Woah!! So somehow your modifications made the code 6 times faster??? That is cool and shocking at the same time 🆒

[ranocha]

In GitLab by @gregorgassner on May 7, 2020, 17:17

I am very much open to a discussion on structure, layout and changes to Trixi. Also, I personally do not mind to restructure the code a lot and also make it more Julia to use the Julia features as good as possible (and hopefully simplify its usage for us, and also students!)

[ranocha]

In GitLab by @sloede on May 7, 2020, 20:33

I am very much open to a discussion on structure, layout and changes to Trixi. Also, I personally do not mind to restructure the code a lot and also make it more Julia to use the Julia features as good as possible (and hopefully simplify its usage for us, and also students!)

I fully support this statement! The main goals of Trixi - which are often in conflict - are that it should be as simple and flexible as possible, such that it remains easy to learn & use (especially for students) but also powerful enough to be used for research. Other than that, I do not feel the need to stick to any fixed paradigm in terms of structure or style, but I am very much in favor of making it as Julian as reasonable.

[ranocha]

In GitLab by @andrewwinters5000 on May 7, 2020, 22:34

I agree with the above statements. My old Fortran habits die hard, though I am adaptable (and willing) to learn and exploit Julia's strengths. I am also shocked that some simple type changes made the code six times faster.

[ranocha]

In GitLab by @ranocha on May 11, 2020, 15:34

mentioned in merge request !43

[ranocha]

In GitLab by @ranocha on May 14, 2020, 07:34

I'll post a possible structure and approach here to have something that we can discuss later (and update here).

• Reduce the usage of modules. I don't really see a reason to use so many modules. Additionally, it makes re-using some code more complicated since all modules have to be imported explicitly.
• Reorganize some structs.
  • From my point of view, the numerical fluxes do not belong to the equation/model but the DG method. For example, stating that a lax_friedrichs_flux is used specifies the DG method and is not related directly to the equation. Additionally, the equation/model should be re-usable for other purposes where no volume and surface numerical fluxes are necessary.
  • Maybe it might make sense to compress the basic structs as much as possible and introduce additional caches for the computations. For example, it might be nice to be able to specify parameters of a DG method (say, polynomial degree, central volume fluxes, Lax-Friedrichs surface flux) and combining that into a discretization for linear advection or the Euler equations. Additionally, the required caches might even depend on more settings, e.g. the time integration scheme (cf. Implement paired Runge-Kutta scheme for non-uniform meshes #21). Then, a function like cache = create_cache(time_algorithm, space_algorithm, mesh, model) could be called at the beginning of the simulation. In this way, we could specialize such a function for the combination we want to support. For example, if a new space_algorithm (say, WENO FD) instead of dg is implemented and a new time_algorithm (say, paired Runge-Kutta), we code only a cache for the combination we want to support, i.e. paried RK with AMR-DG. A method error will be thrown for all other combinations that are not supported and we do not need to blow-up the structs such that they can support all possible combinations of parameters.