From 77f99d7c6f8468c48a2ce6b1b25af95a5485901a Mon Sep 17 00:00:00 2001
From: CompatHelper Julia <compathelper_noreply@julialang.org>
Date: Sat, 10 Aug 2024 01:28:21 +0000
Subject: [PATCH 1/2] CompatHelper: add new compat entry for SciMLBase at
 version 2, (keep existing compat)

---
 Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Project.toml b/Project.toml
index 85723c9f..eebb294f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,6 +16,7 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 [compat]
 BenchmarkTools = "1"
 CUDA = "5"
+SciMLBase = "2"
 StaticArrays = "1"
 StrideArrays = "0.1"
 Trixi = "0.8"

From 9c5b6ad4fad2d992f60c78cd7d6b097ffb371bf2 Mon Sep 17 00:00:00 2001
From: huiyuxie <huiyuxie.sde@gmail.com>
Date: Fri, 9 Aug 2024 16:25:34 -1000
Subject: [PATCH 2/2] Simple test for 2D Linear Advection Equation

---
 .github/workflows/FormatCheck.yml |   2 +-
 examples/advection_basic_1d.jl    |  47 +++++
 examples/advection_basic_2d.jl    |  39 ++++
 out/mesh.h5                       | Bin 0 -> 36832 bytes
 out/solution_000000000.h5         | Bin 0 -> 34816 bytes
 out/solution_000000018.h5         | Bin 0 -> 34816 bytes
 out/solution_000000040.h5         | Bin 0 -> 2560 bytes
 src/TrixiGPU.jl                   |  14 +-
 src/solvers/common.jl             |   6 +-
 src/solvers/dg_1d.jl              |  69 +++---
 src/solvers/dg_2d.jl              | 338 ++++++++++++++++++++++++++++++
 src/solvers/solvers.jl            |  23 ++
 test/test_solvers.jl              |  57 +++--
 13 files changed, 527 insertions(+), 68 deletions(-)
 create mode 100644 examples/advection_basic_1d.jl
 create mode 100644 examples/advection_basic_2d.jl
 create mode 100644 out/mesh.h5
 create mode 100644 out/solution_000000000.h5
 create mode 100644 out/solution_000000018.h5
 create mode 100644 out/solution_000000040.h5
 create mode 100644 src/solvers/dg_2d.jl

diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml
index b544da6b..9d2f4424 100644
--- a/.github/workflows/FormatCheck.yml
+++ b/.github/workflows/FormatCheck.yml
@@ -24,7 +24,7 @@ jobs:
       - name: Install JuliaFormatter and format
         run: |
           julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.56"))'
-          julia  -e 'using JuliaFormatter; format(["src", "test", "utils"])'
+          julia  -e 'using JuliaFormatter; format(["examples", "src", "test", "utils"])'
       - name: Format check
         run: |
           julia -e '
diff --git a/examples/advection_basic_1d.jl b/examples/advection_basic_1d.jl
new file mode 100644
index 00000000..ba25c0f6
--- /dev/null
+++ b/examples/advection_basic_1d.jl
@@ -0,0 +1,47 @@
+# 1D Linear Advection Equation
+using Trixi, TrixiGPU
+using OrdinaryDiffEq
+
+advection_velocity = 1.0f0
+equations = LinearScalarAdvectionEquation1D(advection_velocity)
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = -1.0f0
+coordinates_max = 1.0f0
+
+mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 4,
+                n_cells_max = 30_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+tspan = (0.0f0, 1.0f0)
+
+# FIXME: Remember to export
+ode = TrixiGPU.semidiscretize_gpu(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:l2_error_primitive,
+                                                              :linf_error_primitive))
+
+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 = 0.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, save_everystep = false, callback = callbacks);
+summary_callback()
diff --git a/examples/advection_basic_2d.jl b/examples/advection_basic_2d.jl
new file mode 100644
index 00000000..bab60b39
--- /dev/null
+++ b/examples/advection_basic_2d.jl
@@ -0,0 +1,39 @@
+# 2D Linear Advection Equation
+using Trixi, TrixiGPU
+using OrdinaryDiffEq
+
+advection_velocity = (0.2f0, -0.7f0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0f0, -1.0f0)
+coordinates_max = (1.0f0, 1.0f0)
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 30_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+tspan = (0.0f0, 1.0f0)
+
+# FIXME: Remember to export
+ode = TrixiGPU.semidiscretize_gpu(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution,
+                        stepsize_callback)
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, save_everystep = false, callback = callbacks);
+summary_callback()
diff --git a/out/mesh.h5 b/out/mesh.h5
new file mode 100644
index 0000000000000000000000000000000000000000..e121eda2e52ee474a9a37baf0239763f050d8fb8
GIT binary patch
literal 36832
zcmeI)2b5k_)dt{i5(sHP=%I%mdWTRGLJK7jLg+{c$v~oMCZVZFC?cp-rC5-F(p!)s
zh)ASpK*13s3<4S>!iWJC3l<b-=6fF!FI)%L|Nq`${TbKdo%ig0_I~%-=X`hWcW*pz
z$R2ww*>jbiwYHz`-D`{0`n3P*KX2t{?uGRcoc-AR(hxj1cte-^FAcR-Yi+i^8Sj($
z{kql%_>CC8=bmjN;=@c^1pm)hV91`k4QpG1!?O<EgR`^GLfiLR+wbkp2LA_WzczQP
z2NAb-ZHxASuC<=E?ydi8n9w+N?8s@SPHwC%*QfO)wHh6X5mOo)hqhYV{TIyZR_on1
zZ{nD76Q<Vsv`(m(H0bSl|5Dk4dEGN_<mkrn<KKPeAyU^`ed+8`pEr8c<WZx?O*{2n
zQ*Lo$y~?(E15)3T@{MmCb^O0c`nK-WY5r36`4h%X95>;_2_wfhp42%0U6U^+cN;Fq
zZ~mxL{<Hb@ozk_obZz<8PoIX-jT5IePN}CKxkYW6egofL@NIh4wt1`d^QfVrzHSyC
zF09M%y=#Nof33%F?c?~yi6=}O`|snMT*vkR3*z){t+`d_<i;uECXE@_-lFAO9qW=>
z6O@R{a6_$I-8HRm$S1C8wgl|j`Z?U@y70zpleKu;o@~V~UAnYK>Yg}V8(bJ}s9!7E
zoLhgE%__Dz&w28m<7>6)ZNruSIJf8as7o1CP}Fv6+Y8J%{{i(=nC>o^^v?C|-)?NE
zFH_CB+QdcvpRe}~tJJg4*l2-@cU$kt3!DFIzgw^7?1n*w%<TPp*9}Gce}6rPwEulk
zt98!x?2>a<_7_k$=ebWq?}?4$P8fUKq$wlEjhT8%_Uk*>vv2#t|L>jumPz$azCX_Y
zLhIFY@~A1TAO3BtHGA=T4Yk%!%dAhIhHmYbq^Yy>Nm2KB=X$>Td3TL)-tn%j@!LLM
z-t9WmmZAJ}_M(r}6%1RUj`jb(^K*9W!sh?;-xpf19;3&O8$YIPIr}g1cHV{Ty+%)(
zG-b@ViKC`9PVL+4ZI4C2ELwp@E6`b2z~4d(KELWp{XJ8Lxs9*+@GL`qOJxDdkl%9g
zD?>W#pW?OJ%JDBlerv?94Ee1SzcS=kzbe#fWk~0yiMLq?$ZzZTl_9?Y@hd}qyTq>y
z=^Pxt-8(>jL*iG4{D#J_4EYU@Um4n`hsE#k;1L}lzhmN8hWtjyuMGK(jb9nkIU#-%
zJ3xL@<5!0KPKjR`@;fblWk~0l@jI&n<aciT%8=jr@hd}qm&C6O>AXCCS9E~<u8dz9
z^1CK}WytUP_?4l3dQ1Fn4c^uP^1CyBWytT|_?02Q2jW+TbUqZnhdV%iv*TBW{GNzk
z8S;B7eq~7KGx2Ng0QtQTzcS=EFMegn@8$TFA)WK%_i6{o@AdeVA-^}`SBCuFj9(er
zr}-MGmaj8vt>5>w|CAxW?(r)_eoMx$4EZe^zcQqArTF#j0Qs#RzcS>vcKpha-v;q3
zLpnE!-=-ZPzpdg|hWxgVUm5b-IeukG=Wg*E+yV02JAP%zZ~yp}A-{v-SB71(j}Hsp
z5^Tu4uHk)yn+695{r&3iQ-6Oh7rQS8@-=^LiHvs-Rt>&%bwAxZf|jgWSn@9&Y<>38
z{<BPQ*<jVs_Emb@k1thLtkc>`!K%T1i~EvwBB*~bUk0~6O|Ac|7F<18HLMu{{eo)+
z*A7+<>qfwO!S#b1R4pv|Hw^l+f8+2?I*C5_ty@RWfZ#5{&4cx)jkVgA;adf(hHWBX
z+u(M=?W-1+{5u494DJ-%IaoF98UX`?gMzyStA^bpU`TLiaCmT!U>*Y1_6pxSST*bu
z0s98`3+`XFu;f1=cwlf?@StGTaBu{S2p$qVv}$3=e|T_0@XX-lLEm5bUTSp4#|DoK
z9u+(~_<>;6Ffsy;4UP&P7pxk_L_lM(^+kpDpA&*r!?+0eVDLl1@l^{;{)xe<!Bc{#
z1@myY=52%zXM9SqYM2%QCk9Uno?NxC<Uci-hrhLthJUP+==*qdoF1GWJR?{&oD~7*
z2G0*(5<ENjiQqZGPX?=o^CI9=!A}Q2Q?;<<zaV&F@S@<w!K&f25%9U-rNPUp7MA>1
z1g{KU6TCk7`QR6VUkrXJST%e(0<H?q2wokm8om+%UkzRxysm0t$$vxe#^6oCn}b!u
ztr750@R{KJ;BCP>gZBm>2;Lr?8N4I-wP4k7R|MQ0yeIhes)Z&0eZg-8?+<=6ST%er
z0v-&`3Vu6SH9Q;vvx83rp9($_d^Gr-;CF*n!($Qfc<_6{?^i7>`JW8_AlMZAVX$iW
zQ3T8h{y6w_)xwg$Iru_wUhw7Mv%#MPp9?-8tQvkA0Y3}Q4gNe>HT)t1UJU*+_)^uv
zlD{SRtKciaUk9s(S0mu{;2Xg=gRcdD6Z~!PcfqRR_Yv@i;2(p3s#;j`|2g=V;9rA(
z3swz(kAQy!-y*2NFV&*=L%e_C{RZz3cu&Q9DBd&i9?4><WAR|O;1Zprz8=xhGuSKG
zyOZc!Dms=9_6aW2N%Soj9m@w-2(H*k^m$Lo`$yg{^8V1u(YZ=+)!=HKME@Gmv1YJe
zaIH?FZ=L8^H@IGK{Z68<e{^gZ+$gwlC(-BqH1C~xugrU6n?>j5!7YMYb`t$tN5?k7
zZG+o&5`6=rV~605!JRsZzFneY*Wke5piZLCdx_pR^gf~Y1$U3mJ%W1%_v$42heXFd
z!F_}KbrO9;qvL?!fx%&&MBnh}I5;>W=-*vd``4*|reBVSMIgK-98L)D8+<0XX>e+=
zN6_~NONE~p>=tzYyF~cn!Ja|)#ofbC4&wZm&bT<e!{zA{zF2VOpzk%7EBN;ccmFFM
zU8}_IzIe%Sx^Vu>$1cul;qt5yep-+odN<4X=|T0VcbSY^2l?n(Huf`uIRDKvF3yVK
z@@x^VkL!f%^VS7_^{7`oy4H=qby_K$E}Z{1v5T{QxIEj2cMTpMoF5z<q<4+*H-hTX
zk5x17{#U>Bea+b27q1?!9-O~=#aS!<@~jtbUmg+edyW1D|9;`_f5oHg$k_d#N30!A
z7tVjf*u^<IykBsmpnXXXy}QS*AL>!Bepv_k=-D9i^&RKGN9^Kk94^nE;rci_T%U&&
z{MDmg@#q>8f9teKI9)jZePS2q_;7jl4Zl2SKXwhjHF#yvd2-$cWc;e2^XR;8nQ`aX
zd2yb$iv8*!&VPrDi?dz0JUfQ_9MW(5VV8ox^XR;aN0)x;_x71j7tVjz*u~i^T%Lj9
z*9YmL_qL4R7*vmXcgnbRkdL07W4}3w^S?dg;tUFxXJ)uQ4h`4mI}85mQSYGG={g|X
zI_(xt7ta5#*aroN1?9OreAD2WLD%_*!s&IMoHzC8hx6#X>X*JdFV556iLV}<zk0>l
zFaGiz7;ayl74CCtc){O!bY8`y>+JYD=lh4#h4Vi+c5%)LmuEz{eMt|!567+_>QS$L
zSqJ&(IVkh>9q0c@?BW~}F3+Rk`gnf0KF==rt4F=!(RD%mt<$05bm9CTi(Q<H!sU59
zyh{Q7yDDhEPKclHU+w2LGJZ_Z_XnfH`vs2*y8k^o+;w$ixclM{gs&aM`H#uCIHSVl
zX$(IiI5FsZjj;v)W5eD5ibvO^*xeT&7fu(>e_ZV1d^r5r;0J@w6Fu}gZyUu=J?cF^
z<JLhwdQOPldBgdi9lJOmD&+Y@xIUf|uFvNd{MDmg@#s1={?=)HI9)jZ^I{k0qv7&=
zD*W(Z*P!PDmxt3kHQe_H>Nz&!lY{Pm^-JHU#qPd%O1OG({^}Lyq(YvLgxi;kh5KIP
zw1WSM;qHIMqigZl-4~x6P8ZJq<FSjgM7TVshufF*(0fJf`k@~6>X&tpkDiZZzP{u9
zKOeg|)5GQYLbyII6|T=$7W~zt-s!Q^wRE_3IwPDeod1_&pB`K$D9=^lgM;>?{pz{a
zkl39k=gswE-`Jf;=k=V7JIBt8^YqEs_YdOyFUYt!pAMJj!f>BM`fWd4Qt)>komcVb
z(og;VOy<*t^Z#t@;(Re&p3j9357I-g=Yk_*SC4vK53PfI^jw^IhX!%}o-2xTX(5m2
zjrw>^?E36^C0*)KuXuEQCG)M*Wd&V0f6o=gxi<dtc;0ws#y1VR&U@ZUuk+-*sYgGY
zN9R?)^xb)Jp1zcP>cRP|SDYE~m*=bD_T^^bKBulP_&bl&Gd?3o*XH5Q`PJcc;rwri
zeMWH0pgcE*+n4mv>$#wQs7Jl}WgX<B=eoq#cbvcHisIZ9e|bD_)W-p_>$B&T{M9of
ze&W%!L&mMs&Ea(6{5@9`XQ%kf<9VaMLt&TrLHqQfjN8Z0gsbDuaCP1rF3-$xI_?O+
zJ$P4e=it|ZYc<q=!TH}8yEylR%kz!!+ky`VcL_ca<bQYgprCkkJ(BU=f?p5Ph4cSb
zA<lQg<#{mNb833%ZI0b@P4%ev{@AU9eDr)XcF#d^{?EoP&a7~GeiE*aPlW693k84m
zShrcR)AeMyb^3NVT{!=r#y%_96qM&@;kO1|?+1s^57Ik3d}vTT`tfMShXwUZ-yh5P
zLBa0^)r0d_uQ=a}zdSz(w=Z2!hiCk$g8$>;hXlo=%k}cm*uNj73+Ml%LLBey$TKH)
z`;s1dUyWTq)T3VgvJUdm^TW*7cbxxgv5WKLaCv?cu8-bJ(&yI;{_3%A;?d<j9P9LS
z=F^4q|9$M@cyC9ZKg8~S$$q>w<6VQEyE;$K+q~F4*K{78*XLq)j-40h>G{|_2gUjS
zB6e}+hRgF}xX&T|wjW+D_&blYGCns*mwxK^&ts<x=igF@^V@KFeiiP1jvjg!i`{cg
z^{Drkv0DfE=y@r2_j@@1#bX!em2i2wh3n%R;riU8;IE#!85fVPKWE%J{W_d3oPW<k
zoWF+4(<}U$pzobr=lf=yUgybqQ;&W)kIt)p>AUmdJpC^5)r0d_uQ-2<zdU~lw=cbC
z<8$iGg1__VyoyJc_i&u^KV?2$IRC%LE{^wh<oQSJ_9Z>^t{l65s7Jl}WgX<B=Wm&>
z?>PTeVi)JFaCue@*GJzg>hl@}fAy$WJi5GxW1T$5qYLN1X5x$EdrNux#ePcx4-1F(
ztG`pvjNN{o9qxYc{BYL`*E{!<7sc*+=(_5D@8Z~(XsG|<dV%wIJrSo@{N?e!r|Z7^
zl}i%e`?dVtPhJ|kcyzfRxh(eH$wwE?-+P_ncrQ{O?|V8=^w8_Px!$Qqy{?DWK|Xqx
ziN2?!7w7LhinBue<=H%3AHDaf&)&D=uO9V^N7p)uXPvy?N*B(5>-dYaUbsBlg!{Yn
z*3@}<(BB30x{tahcJ=6o`^jr#*Drl{zjs~i-fvY8&R@Oa^ow43)(y8WZ%f?kgBuk5
z-RIpDyLfcnp7EQ5?(68n`S(wJaqh@Gc{YsQzNClV!LjRydep05)<Hgc*3W!>$NBFb
zyEq$#%d<zgKHeMs`n-3+Up?v-kFNU)>$Y*`(}nXN61zC}XTCi9#P07>`_X<qJpQf+
z&Xe=zdf|HKJUXwght9F{;yi7c{H_-`f7cUnwvE3$JB0fj(r^1==YqfU=)8(YmwxK^
zc9~BX&VQHK#o04lo?XLT_vxYch^(*coqE*kdT1Twqi3hYuLW`bM<%{F1LG&pQL*de
z{;})xF|qSkk9x(UYiQ<Mr$Gf>IRBA}FV2DSm*?2nH_7-j;jZ%&!s&IMoHzC8hx6#X
z>X*JdFV53mg?e!Q>J?|-_{(!ZxP93iy*{T7D)>8(&Z~HIJ)8N?`F@#C7tVip?BYC^
z`SKhbyM0Lyy%XcFAL>!Bepv_k=oyyy`i}FTl=$L|h`&6OW7o%d;rcu^oWFY1D;{0$
z+pN<e$wwE?e_G;;<33NG6JrMnSL>32u0i*qWmuKpBAxXHmG#B=d+pErWjz<|(|X_E
zYhT~1-}|Ni?-!km_UWR1+O`7k_54-kw`iX(+NbY_eOm8Q@xOB)tdGE!a5z64z7`H&
zDe&SrJAAPa@1?^0c?HgHK08jAxbhe0;q1jYxH!-}$Pb6ci`bNUtY>qf{^twz&M9zq
z^|Ir$3h`za=Fco}cJtYB`qVGKI_1aN)y<C6C!YM~%a600&yKgGo~A<Gbn7o&>cZLS
zWXI_gk1q4+z}d}b$L&w+Aip~0$Jy1*j?*We{N~G#vzyP3Hx<@X-B6uST{tveL^{Qz
z%RD-8c5&Eoal|c~Ck}ft4lWKf5AwsI@gmM0_^+R*P2nwB=UL(Q#msPhpIzYW_5nNI
zT!=TPF#q`iXE&c6r_Vl<-#(HbXSc7|ar(rQ-+cLTcJtZsmegZCUo6x=uTbwx1<tNs
zc6@#z-YbRquN64E`Rq7->X%==^5g94Wyk3gPk!^|$Jxzi$D2}5(}VSWJ1e|7Tz{P-
z{l(em#Gf74fARF!eEr4Q&1c8yb6#gg@9aW;oPB1&j?*We{N~G#vzyP3+n-Izr~l32
z`m0|3#o5)%j?W)h?~{1?Yrg*C?B=uM^r>Hd^~#U4tCt<8Pdxd}mmg<0pB<l@b(@>I
zTf&>d=ZD)DuZ7Qy{gnb|w@=ye7Yp%TD$Ji(;OyqJ<MZR+Jh0w3`R9bokF(oX>^Obm
z$#1^=IJ^1mcvJLQ&*pIJX&+h7IkDsH>Sf361M&3T{F#M#?B=uM^r>Hd^~#U4tCt<8
zPdxd}mmg<0pB-;WJuL(4{hlBGVz~a!3)f$q{iT8(*MIT!*L?lO+0AFi>6>5B_evo@
zPTy+<J5HZ?@|!O|&Tc+C-t=I79j#|`xcc>1z50u@tCt<ufARF!eEr4Q&1c8yQ@{M`
zl^<tUFFQ`3c=DSsKhADGJ3cpkgH7RZRydp)4rdp5ahx4C7vjw+%y(|YVK<*0r%T+j
zdE&4a<KW^z^B_MQ8ZTl?>T3$xS9IG)blEpJyL#Di`#?N>H(#G|cJtYB`qVGKI_1aN
z)y<C6C!YM~%a600&yF``9b1CxrdyqKsS9VPlO3m1Ji5%M17|m%9jDJa$gfWMadvgH
z<MfFqzxndx?B=uMErs<|H&iE77Y>aVkxudGGLH_NT^x2?9C6F$iNjuugNp;rgZyx4
Hd;$IwnX#Br

literal 0
HcmV?d00001

diff --git a/out/solution_000000000.h5 b/out/solution_000000000.h5
new file mode 100644
index 0000000000000000000000000000000000000000..1cf48619a6b668688ccab260f799cb4d13992a9c
GIT binary patch
literal 34816
zcmeI*TZ~>+6$kJ$Q!YbWLjhxIV+5Zh0f{1t56Gc9t&IvoDMX^!PU(lOwgoyPglGbU
zMDZRJ74X3aP^tzchNvXQh%rWmltj77EdwG76gyC4h!AnU@9f_Utm))TY0qr6XCAtz
z|MOpat-aSed#$t2Wj?=RaQSQIFPUF-e+CAM{fb4sU;Xn^`YEockjeii!kX0fr8d!u
zt^F0yR~%e)k1bBe7p4A>m{$So6K~bZWy`vV{<TT0ZT9gQSg~yBsog2~aGHmK)c!8T
zPx@~!*1Es<$iO{l#p~4vRGsK|Ve!@;Ft0eESXc}ceVf;9*fd-$ny{-f?W9RlUSZ+F
zVnvdhKS7+3>s;`up^**eZ5}QTTT~Q>PaJeAeQ?9(&d`?A&l%b{v}NhK3p?j@Y0DD*
z*ul!4DLV!#JI+6E<HhSbAK#@V>746IPRSjZ>|cH1(3V{ffI+X6TR^U}v9qbOd1TkJ
zMl+=^jeApPc>U^;i_h=$#@(GWYF2IObWWZ)+Ew>5J6|q$?uLz(+<`BX>-Koe`gioi
zCgk4L`(<9SfAQMxJdbSH)F}=*qMUTkaSMv$UYeK@YAZSaPT4!Jl3zD63+b_PO75U)
ze!6;xM>^;4(jl~UtvIl<B$?BfR?NP>1WW`+?<ckMCXRQH>}o_^-u~S&Pqh90{k;wb
zQh)RM(gYO6rCpw%LlfL3-?{Cl=PoRc>>cXPUoTEK4^7NpH$=N9La)n-cDNUuJz|RM
z>i0eWw|_EY$OPJ|?`|}(qKh10z83W@m<Z$zLu)s7hF8Dsiv9Q8l>h4<syQ4hhfY5H
z`1>9lEAPGatKU2F#AnK5?z(H$FAm;0iC_A{b$yF(d3qB6;sb*>zV6<g)A$c9zT^{w
zr#xG(9k}+qU-$oW8Xx}bw(nf?=B@R|=~-~szYlugs%J}b7p#8t8<(B%+@v1(&1Wt-
z>#?sqSAOINFMRaS{l=&9;ip}9`_3in#%t|sT0gsuC)x);`)4<N@{nIXSF-z|V-Nf7
zLpP33juXD^vYXfJxMRE=`Syz6JhO9r8XsN-_4Bp!)3kne8&9+k{<Zs0x$2V9@sizm
zS~wm#I3Ao{EgX*=T)%KUEgX*=91l)!-EngKN=}@_{=^BtMLv9II{$LVAN!2Mc+Csk
zyu$Ica6EEwJUG2tI378;e&KjpI377T9-Q8~<K&w1e{TKeK|II!gj?6e?NG`ukv~7>
zi7%x56#4L@Q+~^>zs^46lvmj;Psj0GlJfMLl+Pm{UR6^&rmd$jKAVoe-1DDXe{TQN
zzb2heUrOhVxY=9rcXz5Q*QU4>fAq|BT$~^M&btr&Q|IZD+`BusUipE+f0QwvZ}|G@
z_nmdzQ{`!&8U4wRE`M?wAAZ9btIqt~=;LK>{c$|eKI7~==;trq@#XC$yPsHo;dlS|
zp}$X#6F#`D@2HE;`CED7j@yp-^uPW(jSpY4_E)#x^y5cM{^#eop1_-0f6PBWM}NYb
z_CGyytv|Q_bL$tk;yK1AeES*w+uwKhW2Lx#!_wpLeDSk?nT&Jzw^rVE^4lMIq<sGm
z2bSMI{_r$D{OnuKUHpgRAFe$=n%2*5<B9gcmtT3)kq=(}NXhQOYyUj_p6f>^#|dBe
z;x|ux)0)w8@$WZZv}Ez<G(Nnl`$pURpY{RNnX#UPH&uV;UVn4zZ`c3s3cRC~=l|RC
zzq&ByfAw0d7hQj-^XmLtcbuM?@;^QD*sHDoX3GCf>t}Z^|L4}<=Krb=m*?j?|F{0H
zt^?J1b6tPtuK#lDSGTHXV|@z0f5qB+ic-BgRo_Pb1BpjYor?#jR}04@2X9J0J@nGY
zF1WfLZk>pH*9q2>$cOX07Or1%aQ(vRjeK^pJM!ViU3Z+GrshAle)C{njF(+-^9ncb
zk<U;37Wr^~&#k}CKI6PPUDvbwq_n>k$Ma;`C$OGJK3v_ATYsH>#<@K0>#=)f+TV=h
zDbhZndNT6imEXt8HvhNzzq=3YJaWD{FW}Y{xOFG;ty9*m$cI~Z>W<UXuK(4(p7}EG
z<_FGC{07H|TX(c@JaTaT!s(6t*8f$VYJSAIxQCm^y5r=`D>?HX`{O5mi+ni0=hk0m
zpK;DM|GQ5hFB+#j3YQmi>#wtq-R|qLd#3#FKA}7x`S8l`V`=|?>Ut{1Q`7O6d;W9l
z&+Y%*`s?CToPR&v*B7^*SB>#(pZK2iyg}r{?K|hzUuU0jdj6K(o=c76u}^Q`KJwx2
z`?UVA>VI(-;|bnW{N?t4ZvA!fDb79DCvLZ-=f7h-yHD%BZREpOru)9R_1D>FoSy4t
zx99ZYcswWSxxL7TSAKuK&Hw2+5Orp(kKs+#pR-+mm(}e5Y)#Lz?e6~1Z1p!)f9B4=
zx%KB>e`o9eruH9S?fkp%AM^E5={Zq#k@>}kTlaJ8ud|Qcv#tLgOz&UGi^eIB!sQ+M
z)tPaAc31!I)}M9uvAfOx)%#yDpYNUgKlT1!^uM^6ZTz+Vf93m!&NFpf^jF<+a!u8L
zx%J0*j`gWJ*LzjfdiDI|`;C#0M=tW=^lIU;Ke&Elf02(zF7n~@)*aWsx}Kakiv5Wb
zaTED)e%HeFOAfAIIK7e2PIgB=+_>wG(_>!A=|>C4BL~NW)2oH!k%Q|Oj;Dp=k%QyG
z>8(3X&b*Q{@3B9A;<w0$^LuXnb@myj>w0#3FC~sgouF=reE3_^`+vFh*V$*B@+!OK
z={O#FTHcO)c;)x@{=fR)eWCXGzv{gi_eJc7xUUFzpUD0Q+<i*-jo|q3D*o!<Pa)Ud
z|8MXAx6gl0oiENi=ZEvmdFs4{<HPOyI1indu|K%;I`$X&yQ@F;7wuzrQ~R&E_49-K
z{PK{z1aE5pwe^4X9;i5*>-d|k|Ly<V=XYN+#!FNC|Ly*7^<4pVvwhJ$-G7e#sq5{F
z(xa}2tNSB=clUp18-Go$Kh`<%9Oo0>)cyC|`uX2=S@b7-w*GgWXrDjE3*2=l{o-t{
z_1D=a&Ry5D+jV~&k9`9B29XbU-Je^3oqfh>Ux(fHed2iR6WKS4e0b&e_Wu9AUjMc0
z|GxeFSL^?M+yC;vI-|}1)pu*<<)-TodgS?-&*9DH|MvbrKgW6!-c<dW+yA-s*VU)$
zT=!Mgt?F1U9FH6v4^FQZjz<ozUpSr?jz<oT2dB60I5~AaIqO91&pKh<h<rG|YvKAO
z2iGs0-pFSsyCWZN+;zw4F|XwGqlM#<gX6*J)xz<}!SxHr)57t{!SUeq)*UBjUdfsF
z*dIUfTjay}J-7Zk`;60dJ-gi(kK?gUTel-0?!I_#{dM*kr+poEtJC6m>=UWyA|GD)
z-TQxI<xKlOaQ(*p5_Z6McmI)`am;o8&8?sR#j$u7mvHAb+<71Q)(PuI<inl!x%JoC
zXPnkmb{mHljz<oT2d7sH$0G;VFC0$`$0G;FgVS4goE*QB<M-H~I1x9I59fC+T)*Vt
z`i0XQ`RrtO<im}-?l?W>m7IRGa6EEwJUG2tI378;e&KjpI377T9-Q8~<K)aMIrAR-
z<0pQLd^o@7)?a6zaklxt`nx(YpUazYd0L)_%hT}6@ACiD^VBh(nvTEL|JC0?Qs+f~
z!keo9tONASwf?&J6t~{%7q=&*@Bhbm-kQD>uqu5Y3;(3VKQ-Z=qsXnl&OYPpPv7-o
z_v7jN|HcV-pMGijJ{o*`;>YJ(|5yFbpELQNpX=uL)%t&L>94bop6Gx6G){ize}0DJ
z!`=7Et-sDbc0U!q|Nq7rzW;yC4B!99_ng<q(tF>z_4C6__Q5?@$!^bG#_@Pg({r1V
z53l_0d#Y{zS7*lh7~WL<nS1@stv~nrn_GYG^><(Kzd9<`OLOJ_ne_AjOvm?D{#P%m
z3*{a9>526Oys7#lcmB<-pC9D;n9t$O=KryhpIiS|e>aU^nyx?T;pe&be{TJC^{F~{
zw(tMpkyGcYd*SqI;dtcW`i0|Z;dtcWcyM~_j+0Z@lM^SgKkI~eihTI%)8GGV{a?NR
zV;;>nzrgtwj;Dp=k%QyG`K9hSxu*P|Tfcb_&oMsX-v1M~C#3zY$bTa36R3A0AAWuM
iURrMbb@mzOd(-#d*zNCx#qoImt4#Zykq`IX|Nj8R*4M58

literal 0
HcmV?d00001

diff --git a/out/solution_000000018.h5 b/out/solution_000000018.h5
new file mode 100644
index 0000000000000000000000000000000000000000..26c964491c4f766eea11b3b5185726cbf7939fd5
GIT binary patch
literal 34816
zcmeHQdt6P~7eDuwN|LMe_<2p3Mlug(80LqzM@3VTktWf*ON|QM%9Q5eF)5GHAdg|f
zFyoccLwQs*%?P73l1c`Z)I}vyx1{c@bMHQk^KpJ9W&b|o^hc-GS^KPY&RXrY?pbSn
zi=1L^*<G`@CPjZ~Xi%z@Nbw`{O^&~)2qopo@0h_29|ib`*{FR&B|$*-py+dk_`C>T
z-$+GC@Wy=hpIBJX6z2NO5ijEpulW<0VllyveguyAIcVVHS?oXfSB0<9KP#9S;!7ME
zf8C;{ya;n$Eozv8phBrpT9gJQ@b~iZlTaeiOR3uUC#pQ7p%yg-uhRrZSSMZ-<R0j=
zz+Xb?iYV$0ILZ7e+4%U2-2-NLy8F5ZOz;X8d(vzcnEicor9Myf&`|2JaDi`#mw4WD
zHNk&#^vb8})Up1q!R`UizX3*}XV-OO)`@+^eq#T?=dWe%=BYe<-+p3=w`*X?Lb2k$
z>1W1>u@4Z7r@~oU+&}dE=eoH*zRGp#|5!&~@m=qcie<2FwBiR9steVNe$IhDeqySp
z5qe^t5#6W}a(E()SmnAJ#H?rRtD;=*75G0_yOpoB$XpbQp5A9y<@2HCN&>|TpKA~u
zr7-v|pwyL`V08p|#}o)K0h~keg^wz5o<8%uJEHlz(Dw{Sp-`x}fCj#viU2<Viu#)7
zW4?4SnP$$rkUUq5GEtnOpP#}{{oUdD(Jo5=At){hM~T8Y6^z)gl|R4ycj1$W!9)V)
zsJu4q4U85vMU8o0B0)EBWFL1AU$Mk>ctn?%SLI)Qg@E(F9lt7<A-jhU0$OTiNPbw{
zaOIE;Nyl&P-rGxtY~H1;&jVf6;m^YQ92uHAyr4w~=xUqO+`6mEP;h{I%sJq1(%R!T
z4CuLkin<w|tVf&9>~D!2TaVUmdS~_(PZ@geQ#GAODl)Wi{uV74prgz=cQp3ZA;-yS
zjgR5{t|fkF`oi`0M)V#3C7d6r;TUBM=eP73t@#kvjglT5R{3o`iqA4VJOJo=kvdwZ
zgmq}s_WRc)z+Vx&G<6ivjl;9%mu!59tOh#9^ewDKB@N9b^?U2k43Ez)#sa_b5y8yn
zf?8xb#!>Whk6PmLQPX}KJ<xme9TTIopkH}jpZmQ)@8D9`A*rC>63JrCIH2EN73@FR
z_8~I;{l@0wK%Z$fDSfcn160vpJcQ=oJl#h7G|)5BVp@71zmJ5TzM+3zt3ke^sII$#
zZhpG3-T?SRrs)P=0QxqRW_e$J50RF?bJf1UB+$?PC!KYspkJ4$6a9)nzuM3_S176m
zr5wFiZwz$78z1()B)*Sy50-Wj06k>IlIW@C)o6+1326!N@7$Z`DTn(RSW*7f>=9My
zgZWMNhnC+(y0<BdMxbwuD!Vrp_}$Yo3lo4|wo@7#ReT$fmVaBr(7G$!Pc!xV2@B9~
z(XuE5T0bN7X39#?uXoW%=XpTSS*0C!Zfzx6X`a-R);l6|$(SqCq-fj<<Fh%ypMN-E
zdmzwfEs8iG@vJ~Se;C^*;L$Dg=;K8uzXCmQ-{RAkfM4sHsdEg_&CZrfFCQ-B?N8YM
zNZX&ZeD&t&rUdZMHAAcPe*^z4^>P`u4E)ogaawI6_-FNXXO*sYQWW+biXNL^jvUah
z=T8DXK|8Oz`pTPVYnLyo4gkL>{|mEqK%cnwExVzmr6|X3{PHYW3Hox1Z9~tvGBjet
zvWq={|HksAae56UsQ>3hb5%!|Aky-u%}GzK0lhuN?>?*uzlkPixI34D-f1QGj%1H2
zL75BOb0>q|nnN|5E-Wuaqiu3cuL8YdL(Dt2dx}wzuA5v0_cKrIA&UWe`47=~dTzz&
zxO_(+4Z~uTD*1A#4czBPn{I9M1pdYG=Cj2>zb-O5J((&-y#2}gA8GqD^4+ojS$C5E
zjl}-9gynw^-~Zz;UGPu&`b<qz@Xw<MhRFtkf7<tV%}@dV9Hc*`Y$MQDC;Ze{V^)lk
zKXcO80s6-SZr(KpzQ=<9{AvrwY#pb>ALTbt{@!U-wZK;t8mTc1)<sB5w}%72gSuvd
zCGdR`|Bc-`pyyuIdq)IvNy{&g=AE(v{lwqs+6h5F?-ds2HlUxg*803G_<N3z*&FHu
zbPL}ce;v1CG*cq<E&_U^>;9-)GI(w}mkZ_tzrXm#M>~M7@3hhEPuM3`99WWi7xt5*
z@frK4!+m}hd%Ch~LkYTC(==-$>?b>tA8FJ34r6~(`w#Q{RW>N?pGe65q~w#fKWX`w
zhZw)t1^iRB{0E)0;Ln>LSJ|t8f9|n8p-0F6;U>MmH3Yxi?G(|^1me2t&5gTLf$nxJ
z<F};{_dnkedCq|L=dJ#sQ9vIvDr?a=$P=8JO=B}5pBU23!(AKVf95Y+Q-212yHiUg
z>X1)_3~^XP=N+WwpRs!TK>_F;uGM(Q5A@qLFSc(!=&j~ubs+-ub8CsJ+X3{T=&-vB
zA+Ohe%Oi0O&=2Z|-;0C1e|K6~dLi&12wH6w4D{W$CSmmFghGM(9837zuwYMmY8ucZ
zR-C(cf&Scc>DJd?Krge5^U{FN9en)f<3Gv#li&aUzvKV!9G^dK;?Ex;9G^eZ+W7gS
zHT%EOO8!R)J0}0*ivOK0pLF~sEuVD!B`u$H{O9dYHvc4Tf70@i+UorEke`OH&0KpB
z^3%-2AqCcupI)!+Wlq1(uS*#;ClTmAnRBye!F&B9_12SGfIgx6_O(KI@Bc@-vv3^n
zH)*c7rRxPf%B7aMP$wAOe0akUs2e1X>sL)*$9#oZ+)d!uH*{D{*AvSA`u^i&s5_9B
z|EcDtqDuNcvo^@+dQE6?STtSNSsEZenF9L7i>WhNKp${4o9Y8~og{~3mt#P8nRV#A
zF4TQ;ER%-Qb)6sAU%y7ze~#!?>RCda=<>tw5;LK0q%OT`x(nzh&cp`Ob)sE^<jdlL
z9(3%Jf?iN};_Xk?|47@PwETH@hl}=te;&$?7||d6(<ZjR8(r7UNgb2%7Wn6Y#yiPP
zP_LaRl$6tT-5TqDJB@&DR23f)26f+6K?nLA1b)YAOOw$+Khkr;+B;Auc5*Dx?hW<i
z1>WiQbx{9J=|9rxbKu|N5@k^b_2hoGu3B=aJCl~5@>Td{HPGAO?Vw3(K)-;Gy4%rp
zbrYv4+YO<fzM{#xW;W1$q-&z-y85@e7Y*t9`^AW`;6kXon`O^xkiq@9<obO<*WW!V
z(nWPpr_bFTEqn}h`!fdl%ZuPXLyYY_E&zY=q^Y*qKo5;ky-@^pe?I>6@t@@TC(`kk
zw0zS0FCYK;_)jwb?CAKLL^%GxsQk9ZfA3C;|D^XH-u`XN|2oV5r2Vg6KUY)_`Dy-w
zoXR_ppLW~x{h(^dPd}absJ0pMQ$fbDH<N(AZM);vmBlhNWAN<4dY~V=wAio;_+ocR
zk_zGI_MmWv0o8!Ut}D?A2fiEfeN7dxuJXI((i6bnHDU0ktH77n+duax&<(chRLejv
zY59v-{kJW+|F)c?|JIt-e;b7RZyQ<tw;OT)?P^y4tv~L+HN*Y4y#2}gA8GrOmhZso
z|NaU0e~)1Ge}9PkzvFQKw?6Lw-iZ6Z%Q^bLi*f(AJ*)qFGVcE-E#D0He{W&+e=o-U
z-(jr&@1Y$1-&;8Pzx`V2|K{U=$G(3emH)Z6^ZgIM|G(b*e_QR}S??d)n}7Ynk^hb1
z$p17uIREP?`Aqzsj_0q~`QKF5{8elD-$gus)z<uvk<ZM3y{vp@{;RF>NylH(@=3@4
z*6h!o|7*?vc>6y)f7n|7N!tFT<?kXie^rj>|47T{>c1VxG5;6IG5=S{n*Y<q^M8+7
z^MBf``9I$NWadv9|08XG((=2q=KsoGaQ<%_$Nb+#j`_b5JpX5e=l@8{XV2eQ;rXjm
z9P@ww#q)o+aR2RM+<)tX=l`U5{x1dh-=4(tf1~jHRS2HHx`F5a`1(Iz|L5!f<oXXu
zzklQF|9t(Qum5*k{~@XOzwPb+BfbCd_5T;E|F7%d{!jk>Z~NZ=|6AVww$=Wf)qm37
z_y2tVhZ2A7W&1yv_rF~A2mbx95`V6<`VUCIf8*o-zoY)&-u#c>|JxV;J1YJ#`!84h
zwX^qM-u|z*{Ym>D>H0fq`K04NZ~tfUhhKI7DQWwYmQT9=PFg<c`YUgLzW)C_{OMQK
ze?q$cPFg<c`acu@x!(VDcK<2q{EM`F()%xI`K047Y59EoCmVlXcl(p}Khp7+w0zR>
zpSS<B_`k2Z|CF@-Nz1Rt`2PZo|DJ^L-)Ca{{{bxg*P|H!-wNaZw_yDL9E|^8fbsvY
zVf_C{jQ>9i<Nr^``2UwM{y%B?g&6<e8{_}a!}$L$82>*G<NsS?{C@|G|KEV||Ao3r
z{Qm_k{Qs30|Nj`q|MzVp{{O4)KW^{ypW2H3uXM8gJ4$|Q_CLa5|0Znz7X1Go*!E}s
z|3ho`CoMl8<NvF1;QxPs@&EU*@c%tH@c+kg;Qyax;s1wX{D0E&nf}`d7XH6F#{c(V
z;r|<O;Qt@w!2e&^PW*rV{ZHHe|J2LQpOB8fq~-JbzY_nSEB<%({{Onze<sfVW_D2i
z_oDHiOFk2Sy>R|FpOyc0l>L2h{&$`A|DPE7r2LO`{3R`)bo}S-ublrQm4A}9KWX{n
z3E}_W!ubEB<#YAlF5|%ePv*e?-_64R&*#AZ{}JQ=2eR=0Nqznx{rpMGCvAVy^4arW
zQ!)Pkg;wzY`(yn7)osN8KaTPLH)H(&z8v`fuUbCSf2+sB{~wI;|CeF>|9KeycP(rF
ZYaGV^cV^+gUBdYPr&;*_B`o~^zX5sDY5V{H

literal 0
HcmV?d00001

diff --git a/out/solution_000000040.h5 b/out/solution_000000040.h5
new file mode 100644
index 0000000000000000000000000000000000000000..0f18c4428a804c407dd86c76c9e0ce0e6ab85b40
GIT binary patch
literal 2560
zcmeD5aB<`1lHy_j0S*oZ76t(@6Gr@p0tSc*1~?z8AOT7<Lg@}Dy@CnCU}WHD0Le)}
z#RZ_|D=<R@GN9^0f?QoeA~5q|G*sItHNr!{!__eW6a*2_aNvN_C!pyE$_1x25F4x*
zCP)DdORqE7Acnxq<zX-Y3otXVG4L>OFfit&Wabt#2msv$ahe0%DX^>oR?oxW0oBI^
zR0!0UT3DJ`l9``Z%pfelz#tElgsJz*%u7uy3QkVUNi1?qDN9WTDRYG?H*^6T1GF3F
z4i1Pr3i5L*Q&Q8hx&)f%Kspin*rD!^FH0;U>;Z&6ZkWE*oYdUZyb|nsfqnu5gg&Uh
zb5n~m;!7$EQo;TPg&EA>Aw{XFzCc+#`qDCUAo|#m^npy&%P<9N1G>!u%wT3<We@>{
zb4g}yDubYcJv@3$xEW0T149I6Ekxf3dzgM$xH3cZr<CB=JrfM<5&q<agb&Ex;*!(?
ztP!FCu?$EuGO$Bj0(A}}w8Ug&ga|l5Rl;aypgc$hn~y+BSwVgV(o9TDU;{Xy<}ouu
z0{|!mQU?P<P+^c_-?Q7&co-mV1BD+ro#rJ1!w-~1K^XyTB#<r!Cr*eN9uSI3Bon;G
zfYk&r8e#z`6~SBq5@CadmjEL-Fd=6qCgr3S#~XIAQpr8T)tKw=SeCr||K2{^=cDf3
zyC3Wee{cP)ZuiOFWVhPHJ^Y{Trze|iUJj&x^*)Q1z46(eTaH8Q0Ml1{2W!g?^%Y<3
zt75X3a{}q#@o8b*Hs9<+yOb{Rv3$4Ber<F^D&@QVh5So=Ka;=PN2SN^d%5<z{heYT
zzjvR%+nXzMzyAOKyZyp{zE1ambQ({4{$HSYkiKrgEucCB_LwUm{l@Z^6aEABpXE53
zBn>q8@`WP@Yk~HdI~<Y$*?*4Z{1Q<hZL1gY^f=I61~PZ}MS$+ya(bH=$o;#79<^)-
z(o>Qqn`i<3)tjaza}(%KXQ@N(|AGGA*FA3y2g3*Zs;i5)xVpTvXAZyleSPz5`>n!n
znA#t_w0AUX=-2?H?^YjSJ?j6$K07RVPW`W^_H#vymv31A#NMzy)&32To^&+7>Sfzw
z`}N<RM?Cue(0*yu+0Q4<AKHIO-G7YF;-P)$qf`x@B@gWXTJ~`$OFXb&`Ow>qf#HGu
zMyoV6XCN(V^!hOiP`tG;Tge-!j+f`omQz6ed={%j7=ijXI?da65@@bv%|}m`_Q&?G
zgI_Uj1loHha<S$WApJFYCUXJMU5_=lZn+M0=ez%=dqM6`Ig<Hm0+4=~QnY6^&|eHG
ROAh%0{rPxynH&Sq-vAYpj_?2g

literal 0
HcmV?d00001

diff --git a/src/TrixiGPU.jl b/src/TrixiGPU.jl
index 24051033..5087cfb3 100644
--- a/src/TrixiGPU.jl
+++ b/src/TrixiGPU.jl
@@ -1,29 +1,29 @@
 module TrixiGPU
 
-# Include other packages that are used in TrixiGPU.jl (? reorder)
+# Include other packages that are used in TrixiGPU.jl (# FIXME: Remember to reorder)
 # using Reexport: @reexport
 
 using CUDA: @cuda, CuArray, HostKernel,
             threadIdx, blockIdx, blockDim, similar,
             launch_configuration
 
-using Trixi: AbstractEquations, TreeMesh, VolumeIntegralWeakForm,
-             SemidiscretizationHyperbolic,
-             BoundaryConditionPeriodic, DGSEM,
+using Trixi: AbstractEquations, TreeMesh, DGSEM,
+             BoundaryConditionPeriodic, SemidiscretizationHyperbolic,
+             VolumeIntegralWeakForm,
              flux, ntuple, nvariables,
              True, False,
-             compute_coefficients, wrap_array, have_nonconservative_terms
+             wrap_array, compute_coefficients, have_nonconservative_terms
 
 import Trixi: get_node_vars, get_node_coords, get_surface_node_vars
 
+using SciMLBase: ODEProblem, FullSpecialize
+
 using StrideArrays: PtrArray
 
 using StaticArrays: SVector
 
 using SimpleUnPack: @unpack
 
-using SciMLBase: ODEProblem, FullSpecialize
-
 # Include other source files
 include("function.jl")
 include("auxiliary/auxiliary.jl")
diff --git a/src/solvers/common.jl b/src/solvers/common.jl
index bacbbca7..93c5d31e 100644
--- a/src/solvers/common.jl
+++ b/src/solvers/common.jl
@@ -1,7 +1,7 @@
 # Define common functions for solvers
 
 # Copy matrices from host to device
-function copy_to_device!(du::PtrArray, u::PtrArray) # ? PtrArray{Float64}
+function copy_to_device!(du::PtrArray, u::PtrArray)
     du = CuArray{Float32}(zero(du))
     u = CuArray{Float32}(u)
 
@@ -9,8 +9,8 @@ function copy_to_device!(du::PtrArray, u::PtrArray) # ? PtrArray{Float64}
 end
 
 # Copy matrices from device to host 
-function copy_to_host!(du::CuArray, u::CuArray) # ? CuArray{Float32}
-    # ? direct CuArray to PtrArray conversion is impossible
+function copy_to_host!(du::CuArray, u::CuArray)
+    # FIXME: Maybe direct CuArray to PtrArray conversion is possible
     du = PtrArray(Array{Float64}(du))
     u = PtrArray(Array{Float64}(u))
 
diff --git a/src/solvers/dg_1d.jl b/src/solvers/dg_1d.jl
index b93c0dd6..9e1af96d 100644
--- a/src/solvers/dg_1d.jl
+++ b/src/solvers/dg_1d.jl
@@ -1,7 +1,10 @@
 # Everything related to a DG semidiscretization in 1D
 
-# TODO: Please check whether `equations::AbstractEquations{1}` is needed for each function here!!
-# Functions end with `_kernel` are CUDA kernels that are going to be launed by the `@cuda` macro.
+# TODO: Check whether `equations::AbstractEquations{1}` is needed for each function here
+
+# Functions that end with `_kernel` are CUDA kernels that are going to be launched by 
+# the @cuda macro with parameters from the kernel configurator. They are purely run on 
+# the device (i.e., GPU).
 
 # Kernel for calculating fluxes along normal direction
 function flux_kernel!(flux_arr, u, equations::AbstractEquations{1}, flux::Function)
@@ -97,7 +100,8 @@ function interface_flux_kernel!(surface_flux_values, surface_flux_arr, neighbor_
 end
 
 # Kernel for calculating surface integrals
-function surface_integral_kernel!(du, factor_arr, surface_flux_values)
+function surface_integral_kernel!(du, factor_arr, surface_flux_values,
+                                  equations::AbstractEquations{1})
     i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
     j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
     k = (blockIdx().z - 1) * blockDim().z + threadIdx().z
@@ -114,7 +118,7 @@ function surface_integral_kernel!(du, factor_arr, surface_flux_values)
 end
 
 # Kernel for applying inverse Jacobian 
-function jacobian_kernel!(du, inverse_jacobian)
+function jacobian_kernel!(du, inverse_jacobian, equations::AbstractEquations{1})
     i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
     j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
     k = (blockIdx().z - 1) * blockDim().z + threadIdx().z
@@ -126,8 +130,9 @@ function jacobian_kernel!(du, inverse_jacobian)
     return nothing
 end
 
-# Functions begin with `cuda_` are the functions that pack CUDA kernels together, 
-# calling them from the host (i.e., CPU) and running them on the device (i.e., GPU).
+# Functions that begin with `cuda_` are the functions that pack CUDA kernels together to do 
+# partial work in semidiscretization. They are used to invoke kernels from the host (i.e., CPU) 
+# and run them on the device (i.e., GPU).
 
 # Pack kernels for calculating volume integrals
 function cuda_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms, equations,
@@ -203,17 +208,15 @@ function cuda_prolong2boundaries!(u, mesh::TreeMesh{1},
 end
 
 # Dummy function for asserting boundary fluxes
-function cuda_boundary_flux!(t,
-                             mesh::TreeMesh{1},
-                             boundary_condition::BoundaryConditionPeriodic,
-                             equations,
-                             dg::DGSEM,
-                             cache)
+function cuda_boundary_flux!(t, mesh::TreeMesh{1},
+                             boundary_condition::BoundaryConditionPeriodic, equations,
+                             dg::DGSEM, cache)
     @assert iszero(length(cache.boundaries.orientations))
 end
 
 # Pack kernels for calculating surface integrals
-function cuda_surface_integral!(du, mesh::TreeMesh{1}, dg::DGSEM, cache)
+function cuda_surface_integral!(du, mesh::TreeMesh{1}, equations, dg::DGSEM, cache)
+    # FIXME: `surface_integral`
     factor_arr = CuArray{Float32}([
                                       dg.basis.boundary_interpolation[1, 1],
                                       dg.basis.boundary_interpolation[size(du, 2), 2]
@@ -223,21 +226,21 @@ function cuda_surface_integral!(du, mesh::TreeMesh{1}, dg::DGSEM, cache)
     size_arr = CuArray{Float32}(undef, size(du, 1), size(du, 2), size(du, 3))
 
     surface_integral_kernel = @cuda launch=false surface_integral_kernel!(du, factor_arr,
-                                                                          surface_flux_values)
-    surface_integral_kernel(du,
-                            factor_arr,
-                            surface_flux_values;
+                                                                          surface_flux_values,
+                                                                          equations)
+    surface_integral_kernel(du, factor_arr, surface_flux_values, equations;
                             configurator_3d(surface_integral_kernel, size_arr)...,)
 
     return nothing
 end
 
 # Pack kernels for applying Jacobian to reference element
-function cuda_jacobian!(du, mesh::TreeMesh{1}, cache)
+function cuda_jacobian!(du, mesh::TreeMesh{1}, equations, cache)
     inverse_jacobian = CuArray{Float32}(cache.elements.inverse_jacobian)
 
-    jacobian_kernel = @cuda launch=false jacobian_kernel!(du, inverse_jacobian)
-    jacobian_kernel(du, inverse_jacobian; configurator_3d(jacobian_kernel, du)...)
+    jacobian_kernel = @cuda launch=false jacobian_kernel!(du, inverse_jacobian, equations)
+    jacobian_kernel(du, inverse_jacobian, equations;
+                    configurator_3d(jacobian_kernel, du)...)
 
     return nothing
 end
@@ -248,7 +251,7 @@ function cuda_sources!(du, u, t, source_terms::Nothing, equations::AbstractEquat
     return nothing
 end
 
-# Put everything together in single function 
+# Put everything together into a single function 
 # Ref: `rhs!` function in Trixi.jl
 
 function rhs_gpu!(du_cpu, u_cpu, t, mesh::TreeMesh{1}, equations, initial_condition,
@@ -267,9 +270,9 @@ function rhs_gpu!(du_cpu, u_cpu, t, mesh::TreeMesh{1}, equations, initial_condit
 
     cuda_boundary_flux!(t, mesh, boundary_conditions, equations, dg, cache)
 
-    cuda_surface_integral!(du, mesh, dg, cache)
+    cuda_surface_integral!(du, mesh, equations, dg, cache)
 
-    cuda_jacobian!(du, mesh, cache)
+    cuda_jacobian!(du, mesh, equations, cache)
 
     cuda_sources!(du, u, t, source_terms, equations, cache)
 
@@ -278,23 +281,3 @@ function rhs_gpu!(du_cpu, u_cpu, t, mesh::TreeMesh{1}, equations, initial_condit
 
     return nothing
 end
-
-function rhs_gpu!(du_ode, u_ode, semi::SemidiscretizationHyperbolic, t)
-    @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi
-
-    u = wrap_array(u_ode, mesh, equations, solver, cache)
-    du = wrap_array(du_ode, mesh, equations, solver, cache)
-
-    rhs_gpu!(du, u, t, mesh, equations, initial_condition, boundary_conditions,
-             source_terms, solver, cache)
-
-    return nothing
-end
-
-function semidiscretize_gpu(semi::SemidiscretizationHyperbolic, tspan)
-    u0_ode = compute_coefficients(first(tspan), semi)
-
-    iip = true
-    specialize = FullSpecialize
-    return ODEProblem{iip, specialize}(rhs_gpu!, u0_ode, tspan, semi)
-end
diff --git a/src/solvers/dg_2d.jl b/src/solvers/dg_2d.jl
new file mode 100644
index 00000000..65ec5293
--- /dev/null
+++ b/src/solvers/dg_2d.jl
@@ -0,0 +1,338 @@
+# Everything related to a DG semidiscretization in 2D
+
+# Functions that end with `_kernel` are CUDA kernels that are going to be launched by 
+# the @cuda macro with parameters from the kernel configurator. They are purely run on 
+# the device (i.e., GPU).
+
+# Kernel for calculating fluxes along normal direction
+function flux_kernel!(flux_arr1, flux_arr2, u, equations::AbstractEquations{2},
+                      flux::Function)
+    j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    k = (blockIdx().y - 1) * blockDim().y + threadIdx().y
+
+    if (j <= size(u, 2)^2 && k <= size(u, 4))
+        j1 = div(j - 1, size(u, 2)) + 1
+        j2 = rem(j - 1, size(u, 2)) + 1
+
+        u_node = get_node_vars(u, equations, j1, j2, k)
+
+        flux_node1 = flux(u_node, 1, equations)
+        flux_node2 = flux(u_node, 2, equations)
+
+        @inbounds begin
+            for ii in axes(u, 1)
+                flux_arr1[ii, j1, j2, k] = flux_node1[ii]
+                flux_arr2[ii, j1, j2, k] = flux_node2[ii]
+            end
+        end
+    end
+
+    return nothing
+end
+
+# Kernel for calculating weak form
+function weak_form_kernel!(du, derivative_dhat, flux_arr1, flux_arr2)
+    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
+    k = (blockIdx().z - 1) * blockDim().z + threadIdx().z
+
+    if (i <= size(du, 1) && j <= size(du, 2)^2 && k <= size(du, 4))
+        j1 = div(j - 1, size(du, 2)) + 1
+        j2 = rem(j - 1, size(du, 2)) + 1
+
+        @inbounds begin
+            for ii in axes(du, 2)
+                du[i, j1, j2, k] += derivative_dhat[j1, ii] * flux_arr1[i, ii, j2, k]
+                du[i, j1, j2, k] += derivative_dhat[j2, ii] * flux_arr2[i, j1, ii, k]
+            end
+        end
+    end
+
+    return nothing
+end
+
+# Kernel for prolonging two interfaces 
+function prolong_interfaces_kernel!(interfaces_u, u, neighbor_ids, orientations)
+    j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    k = (blockIdx().y - 1) * blockDim().y + threadIdx().y
+
+    if (j <= size(interfaces_u, 2) * size(interfaces_u, 3) && k <= size(interfaces_u, 4))
+        j1 = div(j - 1, size(interfaces_u, 3)) + 1
+        j2 = rem(j - 1, size(interfaces_u, 3)) + 1
+
+        orientation = orientations[k]
+        left_element = neighbor_ids[1, k]
+        right_element = neighbor_ids[2, k]
+
+        @inbounds begin
+            interfaces_u[1, j1, j2, k] = u[j1,
+                                           isequal(orientation, 1) * size(u, 2) + isequal(orientation, 2) * j2,
+                                           isequal(orientation, 1) * j2 + isequal(orientation, 2) * size(u,
+                                                                                                         2),
+                                           left_element]
+            interfaces_u[2, j1, j2, k] = u[j1,
+                                           isequal(orientation, 1) * 1 + isequal(orientation, 2) * j2,
+                                           isequal(orientation, 1) * j2 + isequal(orientation, 2) * 1,
+                                           right_element]
+        end
+    end
+
+    return nothing
+end
+
+# Kernel for calculating surface fluxes 
+function surface_flux_kernel!(surface_flux_arr, interfaces_u, orientations,
+                              equations::AbstractEquations{2}, surface_flux::Any)
+    j2 = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    k = (blockIdx().y - 1) * blockDim().y + threadIdx().y
+
+    if (j2 <= size(surface_flux_arr, 3) && k <= size(surface_flux_arr, 4))
+        u_ll, u_rr = get_surface_node_vars(interfaces_u, equations, j2, k)
+        orientation = orientations[k]
+
+        surface_flux_node = surface_flux(u_ll, u_rr, orientation, equations)
+
+        @inbounds begin
+            for j1j1 in axes(surface_flux_arr, 2)
+                surface_flux_arr[1, j1j1, j2, k] = surface_flux_node[j1j1]
+            end
+        end
+    end
+
+    return nothing
+end
+
+# Kernel for setting interface fluxes
+function interface_flux_kernel!(surface_flux_values, surface_flux_arr, neighbor_ids,
+                                orientations)
+    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
+    k = (blockIdx().z - 1) * blockDim().z + threadIdx().z
+
+    if (i <= size(surface_flux_values, 1) &&
+        j <= size(surface_flux_arr, 3) &&
+        k <= size(surface_flux_arr, 4))
+        left_id = neighbor_ids[1, k]
+        right_id = neighbor_ids[2, k]
+
+        left_direction = 2 * orientations[k]
+        right_direction = 2 * orientations[k] - 1
+
+        @inbounds begin
+            surface_flux_values[i, j, left_direction, left_id] = surface_flux_arr[1, i, j,
+                                                                                  k]
+            surface_flux_values[i, j, right_direction, right_id] = surface_flux_arr[1, i, j,
+                                                                                    k]
+        end
+    end
+
+    return nothing
+end
+
+# Kernel for calculating surface integrals
+function surface_integral_kernel!(du, factor_arr, surface_flux_values,
+                                  equations::AbstractEquations{2})
+    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
+    k = (blockIdx().z - 1) * blockDim().z + threadIdx().z
+
+    if (i <= size(du, 1) && j <= size(du, 2)^2 && k <= size(du, 4))
+        j1 = div(j - 1, size(du, 2)) + 1
+        j2 = rem(j - 1, size(du, 2)) + 1
+
+        @inbounds begin
+            du[i, j1, j2, k] -= (surface_flux_values[i, j2, 1, k] * isequal(j1, 1) +
+                                 surface_flux_values[i, j1, 3, k] * isequal(j2, 1)) *
+                                factor_arr[1]
+            du[i, j1, j2, k] += (surface_flux_values[i, j2, 2, k] *
+                                 isequal(j1, size(du, 2)) +
+                                 surface_flux_values[i, j1, 4, k] *
+                                 isequal(j2, size(du, 2))) * factor_arr[2]
+        end
+    end
+
+    return nothing
+end
+
+# CUDA kernel for applying inverse Jacobian 
+function jacobian_kernel!(du, inverse_jacobian, equations::AbstractEquations{2})
+    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
+    j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
+    k = (blockIdx().z - 1) * blockDim().z + threadIdx().z
+
+    if (i <= size(du, 1) && j <= size(du, 2)^2 && k <= size(du, 4))
+        j1 = div(j - 1, size(du, 2)) + 1
+        j2 = rem(j - 1, size(du, 2)) + 1
+
+        @inbounds du[i, j1, j2, k] *= -inverse_jacobian[k]
+    end
+
+    return nothing
+end
+
+# Functions that begin with `cuda_` are the functions that pack CUDA kernels together to do 
+# partial work in semidiscretization. They are used to invoke kernels from the host (i.e., CPU) 
+# and run them on the device (i.e., GPU).
+
+# Pack kernels for calculating volume integrals
+function cuda_volume_integral!(du, u, mesh::TreeMesh{2}, nonconservative_terms, equations,
+                               volume_integral::VolumeIntegralWeakForm, dg::DGSEM)
+    derivative_dhat = CuArray{Float32}(dg.basis.derivative_dhat)
+    flux_arr1 = similar(u)
+    flux_arr2 = similar(u)
+
+    size_arr = CuArray{Float32}(undef, size(u, 2)^2, size(u, 4))
+
+    flux_kernel = @cuda launch=false flux_kernel!(flux_arr1, flux_arr2, u, equations, flux)
+    flux_kernel(flux_arr1, flux_arr2, u, equations;
+                configurator_2d(flux_kernel, size_arr)...,)
+
+    size_arr = CuArray{Float32}(undef, size(du, 1), size(du, 2)^2, size(du, 4))
+
+    weak_form_kernel = @cuda launch=false weak_form_kernel!(du, derivative_dhat, flux_arr1,
+                                                            flux_arr2)
+    weak_form_kernel(du, derivative_dhat, flux_arr1, flux_arr2;
+                     configurator_3d(weak_form_kernel, size_arr)...,)
+
+    return nothing
+end
+
+# Pack kernels for prolonging solution to interfaces
+function cuda_prolong2interfaces!(u, mesh::TreeMesh{2}, cache)
+    neighbor_ids = CuArray{Int64}(cache.interfaces.neighbor_ids)
+    orientations = CuArray{Int64}(cache.interfaces.orientations)
+    interfaces_u = CuArray{Float32}(cache.interfaces.u)
+
+    size_arr = CuArray{Float32}(undef,
+                                size(interfaces_u, 2) * size(interfaces_u, 3),
+                                size(interfaces_u, 4))
+
+    prolong_interfaces_kernel = @cuda launch=false prolong_interfaces_kernel!(interfaces_u,
+                                                                              u,
+                                                                              neighbor_ids,
+                                                                              orientations)
+    prolong_interfaces_kernel(interfaces_u, u, neighbor_ids, orientations;
+                              configurator_2d(prolong_interfaces_kernel, size_arr)...,)
+
+    cache.interfaces.u = interfaces_u  # copy back to host automatically
+
+    return nothing
+end
+
+# Pack kernels for calculating interface fluxes
+function cuda_interface_flux!(mesh::TreeMesh{2}, nonconservative_terms::False, equations,
+                              dg::DGSEM, cache)
+    surface_flux = dg.surface_integral.surface_flux
+
+    neighbor_ids = CuArray{Int64}(cache.interfaces.neighbor_ids)
+    orientations = CuArray{Int64}(cache.interfaces.orientations)
+    interfaces_u = CuArray{Float32}(cache.interfaces.u)
+    surface_flux_arr = CuArray{Float32}(undef, 1, size(interfaces_u)[2:end]...)
+    surface_flux_values = CuArray{Float32}(cache.elements.surface_flux_values)
+
+    size_arr = CuArray{Float32}(undef, size(interfaces_u, 3), size(interfaces_u, 4))
+
+    surface_flux_kernel = @cuda launch=false surface_flux_kernel!(surface_flux_arr,
+                                                                  interfaces_u,
+                                                                  orientations, equations,
+                                                                  surface_flux)
+    surface_flux_kernel(surface_flux_arr, interfaces_u, orientations, equations,
+                        surface_flux; configurator_2d(surface_flux_kernel, size_arr)...,)
+
+    size_arr = CuArray{Float32}(undef, size(surface_flux_values, 1), size(interfaces_u, 3),
+                                size(interfaces_u, 4))
+
+    interface_flux_kernel = @cuda launch=false interface_flux_kernel!(surface_flux_values,
+                                                                      surface_flux_arr,
+                                                                      neighbor_ids,
+                                                                      orientations)
+    interface_flux_kernel(surface_flux_values, surface_flux_arr, neighbor_ids, orientations;
+                          configurator_3d(interface_flux_kernel, size_arr)...,)
+
+    cache.elements.surface_flux_values = surface_flux_values # copy back to host automatically
+
+    return nothing
+end
+
+# Dummy function for asserting boundaries
+function cuda_prolong2boundaries!(u, mesh::TreeMesh{2},
+                                  boundary_condition::BoundaryConditionPeriodic, cache)
+    @assert iszero(length(cache.boundaries.orientations))
+end
+
+# Dummy function for asserting boundary fluxes
+function cuda_boundary_flux!(t, mesh::TreeMesh{2},
+                             boundary_condition::BoundaryConditionPeriodic, equations,
+                             dg::DGSEM, cache)
+    @assert iszero(length(cache.boundaries.orientations))
+end
+
+# Pack kernels for calculating surface integrals
+function cuda_surface_integral!(du, mesh::TreeMesh{2}, equations, dg::DGSEM, cache)
+    # FIXME: `surface_integral::SurfaceIntegralWeakForm`
+    factor_arr = CuArray{Float32}([
+                                      dg.basis.boundary_interpolation[1, 1],
+                                      dg.basis.boundary_interpolation[size(du, 2), 2]
+                                  ])
+    surface_flux_values = CuArray{Float32}(cache.elements.surface_flux_values)
+
+    size_arr = CuArray{Float32}(undef, size(du, 1), size(du, 2)^2, size(du, 4))
+
+    surface_integral_kernel = @cuda launch=false surface_integral_kernel!(du, factor_arr,
+                                                                          surface_flux_values,
+                                                                          equations)
+    surface_integral_kernel(du, factor_arr, surface_flux_values, equations;
+                            configurator_3d(surface_integral_kernel, size_arr)...,)
+
+    return nothing
+end
+
+# Pack kernels for applying Jacobian to reference element
+function cuda_jacobian!(du, mesh::TreeMesh{2}, equations, cache)
+    inverse_jacobian = CuArray{Float32}(cache.elements.inverse_jacobian)
+
+    size_arr = CuArray{Float32}(undef, size(du, 1), size(du, 2)^2, size(du, 4))
+
+    jacobian_kernel = @cuda launch=false jacobian_kernel!(du, inverse_jacobian, equations)
+    jacobian_kernel(du, inverse_jacobian, equations;
+                    configurator_3d(jacobian_kernel, size_arr)...)
+
+    return nothing
+end
+
+# Dummy function returning nothing              
+function cuda_sources!(du, u, t, source_terms::Nothing, equations::AbstractEquations{2},
+                       cache)
+    return nothing
+end
+
+# Put everything together into a single function 
+# Ref: `rhs!` function in Trixi.jl
+
+function rhs_gpu!(du_cpu, u_cpu, t, mesh::TreeMesh{2}, equations, initial_condition,
+                  boundary_conditions, source_terms::Source, dg::DGSEM,
+                  cache) where {Source}
+    du, u = copy_to_device!(du_cpu, u_cpu)
+
+    cuda_volume_integral!(du, u, mesh, have_nonconservative_terms(equations), equations,
+                          dg.volume_integral, dg)
+
+    cuda_prolong2interfaces!(u, mesh, cache)
+
+    cuda_interface_flux!(mesh, have_nonconservative_terms(equations), equations, dg, cache)
+
+    cuda_prolong2boundaries!(u, mesh, boundary_conditions, cache)
+
+    cuda_boundary_flux!(t, mesh, boundary_conditions, equations, dg, cache)
+
+    cuda_surface_integral!(du, mesh, equations, dg, cache)
+
+    cuda_jacobian!(du, mesh, equations, cache)
+
+    cuda_sources!(du, u, t, source_terms, equations, cache)
+
+    du_computed, _ = copy_to_host!(du, u)
+    du_cpu .= du_computed
+
+    return nothing
+end
diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl
index da6623dc..a1373776 100644
--- a/src/solvers/solvers.jl
+++ b/src/solvers/solvers.jl
@@ -1,2 +1,25 @@
 include("common.jl")
 include("dg_1d.jl")
+include("dg_2d.jl")
+
+# Ref: `rhs!` function in Trixi.jl
+function rhs_gpu!(du_ode, u_ode, semi::SemidiscretizationHyperbolic, t)
+    @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi
+
+    u = wrap_array(u_ode, mesh, equations, solver, cache)
+    du = wrap_array(du_ode, mesh, equations, solver, cache)
+
+    rhs_gpu!(du, u, t, mesh, equations, initial_condition, boundary_conditions,
+             source_terms, solver, cache)
+
+    return nothing
+end
+
+# Ref: `semidiscretize` function in Trixi.jl
+function semidiscretize_gpu(semi::SemidiscretizationHyperbolic, tspan)
+    u0_ode = compute_coefficients(first(tspan), semi)
+
+    iip = true
+    specialize = FullSpecialize
+    return ODEProblem{iip, specialize}(rhs_gpu!, u0_ode, tspan, semi)
+end
diff --git a/test/test_solvers.jl b/test/test_solvers.jl
index 76467056..881b3454 100644
--- a/test/test_solvers.jl
+++ b/test/test_solvers.jl
@@ -4,6 +4,8 @@ using Test
 
 # @testset "Test solver functions" begin
 
+# Test step by step
+
 # The header part for testing basic kernels in 1D
 # advection_velocity = 1.0f0
 # equations = LinearScalarAdvectionEquation1D(advection_velocity)
@@ -54,32 +56,59 @@ using Test
 
 # TrixiGPU.cuda_sources!(du, u, t, source_terms, equations, cache)
 
-advection_velocity = 1.0f0
-equations = LinearScalarAdvectionEquation1D(advection_velocity)
-
-solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+# 2D Linear Advection Equation
 
-coordinates_min = -1.0f0
-coordinates_max = 1.0f0
+# The header part for testing basic kernels in 2D
+advection_velocity = (0.2f0, -0.7f0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
 
+coordinates_min = (-1.0f0, -1.0f0)
+coordinates_max = (1.0f0, 1.0f0)
 mesh = TreeMesh(coordinates_min,
                 coordinates_max,
                 initial_refinement_level = 4,
                 n_cells_max = 30_000)
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
 semi = SemidiscretizationHyperbolic(mesh,
                                     equations,
                                     initial_condition_convergence_test,
                                     solver)
 
+@unpack mesh,
+equations,
+initial_condition,
+boundary_conditions,
+source_terms,
+solver,
+cache = semi
+
+t = 0.0f0
 tspan = (0.0f0, 1.0f0)
 
-ode_gpu = TrixiGPU.semidiscretize_gpu(semi, tspan)
+ode = semidiscretize(semi, tspan)
+u_ode = copy(ode.u0)
+du_ode = similar(u_ode)
+u = Trixi.wrap_array(u_ode, mesh, equations, solver, cache)
+du = Trixi.wrap_array(du_ode, mesh, equations, solver, cache)
+
+du, u = TrixiGPU.copy_to_device!(du, u)
+
+TrixiGPU.cuda_volume_integral!(du, u, mesh,
+                               Trixi.have_nonconservative_terms(equations), equations,
+                               solver.volume_integral, solver)
+
+TrixiGPU.cuda_prolong2interfaces!(u, mesh, cache)
+
+TrixiGPU.cuda_interface_flux!(mesh, Trixi.have_nonconservative_terms(equations), equations,
+                              solver, cache)
+
+TrixiGPU.cuda_prolong2boundaries!(u, mesh, boundary_conditions, cache)
+
+TrixiGPU.cuda_boundary_flux!(t, mesh, boundary_conditions, equations, solver, cache)
+
+TrixiGPU.cuda_surface_integral!(du, mesh, equations, solver, cache)
+
+TrixiGPU.cuda_jacobian!(du, mesh, equations, cache)
 
-sol_gpu = OrdinaryDiffEq.solve(ode_gpu,
-                               BS3(),
-                               adaptive = false,
-                               dt = 0.01;
-                               abstol = 1.0e-6,
-                               reltol = 1.0e-6,
-                               ode_default_options()...,)
+TrixiGPU.cuda_sources!(du, u, t, source_terms, equations, cache)