18561934460

Committed 16 Oct 2025 12:54PM UTC coverage: 65.961% (-0.03%) from 65.988%

Build # 18561934460

Build Type

push

github

Committed by

web-flow

Commit Message

Removed `test_case` and made several dispatches more precise (#104)

* Removed `test_case`.

* Made several dispatches more precise.

* Removed `examples.md` (was comitted by accident).

* Renamed `AbstractVariable` to `AbstractPredictand` and removed some of its subtypes (these are direct subtypes of `Any` now).

* Some simplifications and additional comments in initialize_rays!.

* Corrected a few typos.

* Moved the namelist parameter "model" to AtmosphereNamelist and removed SettingNamelist. Reformatted the code.

Run Details

104 of 119 new or added lines in 44 files covered. (87.39%)

4 existing lines in 4 files now uncovered.

4951 of 7506 relevant lines covered (65.96%)

1423740.77 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

93.6

/src/MSGWaM/RayUpdate/initialize_rays!.jl

"""
```julia
initialize_rays!(state::State)
```

Complete the initialization of MSGWaM by dispatching to a WKB-mode-specific method.

```julia
initialize_rays!(state::State, wkb_mode::NoWKB)
```

Return for non-WKB configurations.

```julia
initialize_rays!(
    state::State,
    wkb_mode::Union{SteadyState, SingleColumn, MultiColumn},
)
```

Complete the initialization of MSGWaM.

In each grid cell, `wave_modes` wave modes are computed, using `state.namelists.wkb.initial_wave_field`, as well as `activate_orographic_source!` for mountain waves. For each of these modes, `nrx * nry * nrz * nrk * nrl * nrm` ray volumes are then defined such that they evenly divide the volume one would get for `nrx = nry = nrz = nrk = nrl = nrm = 1` (the parameters are taken from `state.namelists.wkb`). Finally, the maximum group velocities are determined for the corresponding CFL condition that is used in the computation of the time step.

# Arguments

  - `state`: Model state.

  - `wkb_mode`: Approximations used by MSGWaM.

# See also

  - [`PinCFlow.MSGWaM.RaySources.activate_orographic_source!`](@ref)

  - [`PinCFlow.MSGWaM.Interpolation.interpolate_stratification`](@ref)

  - [`PinCFlow.MSGWaM.Interpolation.interpolate_mean_flow`](@ref)
"""
function initialize_rays! end

function initialize_rays!(state::State)
    (; wkb_mode) = state.namelists.wkb
    initialize_rays!(state, wkb_mode)
    return
end

function initialize_rays!(state::State, wkb_mode::NoWKB)
    return
end

function initialize_rays!(
    state::State,
    wkb_mode::Union{SteadyState, SingleColumn, MultiColumn},
)
    (; x_size, y_size) = state.namelists.domain
    (; coriolis_frequency) = state.namelists.atmosphere
    (;
        nrx,
        nry,
        nrz,
        nrk,
        nrl,
        nrm,
        wave_modes,
        dkr_factor,
        dlr_factor,
        dmr_factor,
        wkb_mode,
        wave_modes,
        initial_wave_field,
    ) = state.namelists.wkb
    (; lref, tref, rhoref, uref) = state.constants
    (; comm, master, nxx, nyy, nzz, io, jo, ko, i0, i1, j0, j1, k0, k1) =
        state.domain
    (; dx, dy, dz, x, y, zc, jac) = state.grid
    (;
        nray_max,
        nray_wrk,
        n_sfc,
        nray,
        rays,
        surface_indices,
        cgx_max,
        cgy_max,
        cgz_max,
    ) = state.wkb

    # Set Coriolis parameter.
    fc = coriolis_frequency * tref

    # Initialize local arrays.
    omi_ini = zeros(wave_modes, nxx, nyy, nzz)
    wnk_ini = zeros(wave_modes, nxx, nyy, nzz)
    wnl_ini = zeros(wave_modes, nxx, nyy, nzz)
    wnm_ini = zeros(wave_modes, nxx, nyy, nzz)
    wad_ini = zeros(wave_modes, nxx, nyy, nzz)

    # Compute initial wavenumbers, intrinsic frequencies and wave-action
    # densities with initial_wave_field.
    if wkb_mode != SteadyState()
        for k in k0:k1, j in j0:j1, i in i0:i1, alpha in 1:wave_modes
            (kdim, ldim, mdim, omegadim, adim) = initial_wave_field(
                alpha,
                x[io + i] * lref,
                y[jo + j] * lref,
                zc[i, j, k] * lref,
            )
            wnk_ini[alpha, i, j, k] = kdim * lref
            wnl_ini[alpha, i, j, k] = ldim * lref
            wnm_ini[alpha, i, j, k] = mdim * lref
            omi_ini[alpha, i, j, k] = omegadim * tref
            wad_ini[alpha, i, j, k] = adim / rhoref / uref^2 / tref
        end
    else
        println(
            "Warning: MSGWaM's steady-state mode currently ignores non-orographic initializations!",
        )
        println("")
    end

    # Add orographic wave modes.
    activate_orographic_source!(
        state,
        omi_ini,
        wnk_ini,
        wnl_ini,
        wnm_ini,
        wad_ini,
    )

    # Set initial spectral extents (these will be overwritten in the loop).
    dk_ini_nd = 0.0
    dl_ini_nd = 0.0
    dm_ini_nd = 0.0

    # Set vertical index bounds.
    kmin = ko == 0 ? k0 - 1 : k0
    kmax = k1

    # Loop over all grid cells with ray volumes.
    @ivy for k in kmin:kmax, j in j0:j1, i in i0:i1
        r = 0
        s = 0

        # Loop over all ray volumes within a spatial cell.
        for ix in 1:nrx,
            ik in 1:nrk,
            jy in 1:nry,
            jl in 1:nrl,
            kz in 1:nrz,
            km in 1:nrm,
            alpha in 1:wave_modes

            # Set ray-volume indices.
            if ko == 0 && k == k0 - 1
                s += 1

                # Set surface indices.
                surface_indices.ixs[s] = ix
                surface_indices.jys[s] = jy
                surface_indices.kzs[s] = kz
                surface_indices.iks[s] = ik
                surface_indices.jls[s] = jl
                surface_indices.kms[s] = km
                surface_indices.alphas[s] = alpha

                # Set surface ray-volume index.
                if wad_ini[alpha, i, j, k] == 0.0
                    surface_indices.rs[s, i, j] = -1
                    continue
                else
                    r += 1
                    surface_indices.rs[s, i, j] = r
                end
            else
                if wad_ini[alpha, i, j, k] == 0.0
                    continue
                end
                r += 1
            end

            # Set ray-volume positions.
            rays.x[r, i, j, k] = (x[io + i] - 0.5 * dx + (ix - 0.5) * dx / nrx)
            rays.y[r, i, j, k] = (y[jo + j] - 0.5 * dy + (jy - 0.5) * dy / nry)
            rays.z[r, i, j, k] = (
                zc[i, j, k] - 0.5 * jac[i, j, k] * dz +
                (kz - 0.5) * jac[i, j, k] * dz / nrz
            )

            xr = rays.x[r, i, j, k]
            yr = rays.y[r, i, j, k]
            zr = rays.z[r, i, j, k]

            # Check if ray volume is too low.
            if zr < -dz
                error(
                    "Error in initialize_rays!: Ray volume",
                    r,
                    "at",
                    i,
                    j,
                    k,
                    "is too low!",
                )
            end

            # Compute local stratification.
            n2r = interpolate_stratification(zr, state, N2())

            # Set spatial extents.
            rays.dxray[r, i, j, k] = dx / nrx
            rays.dyray[r, i, j, k] = dy / nry
            rays.dzray[r, i, j, k] = jac[i, j, k] * dz / nrz

            wnk0 = wnk_ini[alpha, i, j, k]
            wnl0 = wnl_ini[alpha, i, j, k]
            wnm0 = wnm_ini[alpha, i, j, k]

            # Ensure correct wavenumber extents.
            if x_size > 1
                dk_ini_nd = dkr_factor * sqrt(wnk0^2 + wnl0^2)
            end
            if y_size > 1
                dl_ini_nd = dlr_factor * sqrt(wnk0^2 + wnl0^2)
            end
            if wnm0 == 0.0
                error("Error in WKB: wnm0 = 0!")
            else
                dm_ini_nd = dmr_factor * abs(wnm0)
            end

            # Set ray-volume wavenumbers.
            rays.k[r, i, j, k] =
                (wnk0 - 0.5 * dk_ini_nd + (ik - 0.5) * dk_ini_nd / nrk)
            rays.l[r, i, j, k] =
                (wnl0 - 0.5 * dl_ini_nd + (jl - 0.5) * dl_ini_nd / nrl)
            rays.m[r, i, j, k] =
                (wnm0 - 0.5 * dm_ini_nd + (km - 0.5) * dm_ini_nd / nrm)

            # Set spectral extents.
            rays.dkray[r, i, j, k] = dk_ini_nd / nrk
            rays.dlray[r, i, j, k] = dl_ini_nd / nrl
            rays.dmray[r, i, j, k] = dm_ini_nd / nrm

            # Set spectral volume.
            pspvol = dm_ini_nd
            if x_size > 1
                pspvol = pspvol * dk_ini_nd
            end
            if y_size > 1
                pspvol = pspvol * dl_ini_nd
            end

            # Set phase-space wave-action density.
            rays.dens[r, i, j, k] = wad_ini[alpha, i, j, k] / pspvol

            # Interpolate winds to ray-volume position.
            uxr = interpolate_mean_flow(xr, yr, zr, state, U())
            vyr = interpolate_mean_flow(xr, yr, zr, state, V())
            wzr = interpolate_mean_flow(xr, yr, zr, state, W())

            wnrk = rays.k[r, i, j, k]
            wnrl = rays.l[r, i, j, k]
            wnrm = rays.m[r, i, j, k]
            wnrh = sqrt(wnrk^2 + wnrl^2)
            omir = omi_ini[alpha, i, j, k]

            # Compute maximum group velocities.
            cgirx = wnrk * (n2r - omir^2) / (omir * (wnrh^2 + wnrm^2))
            if abs(uxr + cgirx) > abs(cgx_max[])
                cgx_max[] = abs(uxr + cgirx)
            end
            cgiry = wnrl * (n2r - omir^2) / (omir * (wnrh^2 + wnrm^2))
            if abs(vyr + cgiry) > abs(cgy_max[])
                cgy_max[] = abs(vyr + cgiry)
            end
            cgirz = -wnrm * (omir^2 - fc^2) / (omir * (wnrh^2 + wnrm^2))
            if abs(wzr + cgirz) > abs(cgz_max[i, j, k])
                cgz_max[i, j, k] = max(cgz_max[i, j, k], abs(wzr + cgirz))
            end
        end

        # Set ray-volume count.
        nray[i, j, k] = r
        if r > nray_wrk
            error(
                "Error in initialize_rays!: nray",
                [i, j, k],
                " > nray_wrk =",
                nray_wrk,
            )
        end

        # Check if surface ray-volume count is correct.
        if ko == 0 && k == k0 - 1
            if s != n_sfc
                error(
                    "Error in initialize_rays!: s =",
                    s,
                    "/= n_sfc =",
                    n_sfc,
                    "at (i, j, k) = ",
                    (i, j, k),
                )
            end
        end
    end

    # Compute global ray-volume count.
    @ivy local_sum = sum(nray[i0:i1, j0:j1, kmin:kmax])
    global_sum = MPI.Allreduce(local_sum, +, comm)

    # Print information.
    if master
        println("MSGWaM:")
        println("Global ray-volume count: ", global_sum)
        println("Maximum number of ray volumes per cell: ", nray_max)
        println("")
    end

    return
end

1	"""
2	```julia
3	initialize_rays!(state::State)
4	```
5
6	Complete the initialization of MSGWaM by dispatching to a WKB-mode-specific method.
7
8	```julia
9	initialize_rays!(state::State, wkb_mode::NoWKB)
10	```
11
12	Return for non-WKB configurations.
13
14	```julia
15	initialize_rays!(
16	state::State,
17	wkb_mode::Union{SteadyState, SingleColumn, MultiColumn},
18	)
19	```
20
21	Complete the initialization of MSGWaM.
22
23	In each grid cell, `wave_modes` wave modes are computed, using `state.namelists.wkb.initial_wave_field`, as well as `activate_orographic_source!` for mountain waves. For each of these modes, `nrx * nry * nrz * nrk * nrl * nrm` ray volumes are then defined such that they evenly divide the volume one would get for `nrx = nry = nrz = nrk = nrl = nrm = 1` (the parameters are taken from `state.namelists.wkb`). Finally, the maximum group velocities are determined for the corresponding CFL condition that is used in the computation of the time step.
24
25	# Arguments
26
27	- `state`: Model state.
28
29	- `wkb_mode`: Approximations used by MSGWaM.
30
31	# See also
32
33	- [`PinCFlow.MSGWaM.RaySources.activate_orographic_source!`](@ref)
34
35	- [`PinCFlow.MSGWaM.Interpolation.interpolate_stratification`](@ref)
36
37	- [`PinCFlow.MSGWaM.Interpolation.interpolate_mean_flow`](@ref)
38	"""
39	function initialize_rays! end
40
41	function initialize_rays!(state::State)	6✔
42	(; wkb_mode) = state.namelists.wkb	6✔
43	initialize_rays!(state, wkb_mode)	6✔
44	return	6✔
45	end
46
47	function initialize_rays!(state::State, wkb_mode::NoWKB)	4✔
48	return	4✔
49	end
50
51	function initialize_rays!(	2✔
52	state::State,
53	wkb_mode::Union{SteadyState, SingleColumn, MultiColumn},
54	)
55	(; x_size, y_size) = state.namelists.domain	2✔
56	(; coriolis_frequency) = state.namelists.atmosphere	2✔
57	(;	2✔
58	nrx,
59	nry,
60	nrz,
61	nrk,
62	nrl,
63	nrm,
64	wave_modes,
65	dkr_factor,
66	dlr_factor,
67	dmr_factor,
68	wkb_mode,
69	wave_modes,
70	initial_wave_field,
71	) = state.namelists.wkb
72	(; lref, tref, rhoref, uref) = state.constants	2✔
73	(; comm, master, nxx, nyy, nzz, io, jo, ko, i0, i1, j0, j1, k0, k1) =	2✔
74	state.domain
75	(; dx, dy, dz, x, y, zc, jac) = state.grid	2✔
76	(;	2✔
77	nray_max,
78	nray_wrk,
79	n_sfc,
80	nray,
81	rays,
82	surface_indices,
83	cgx_max,
84	cgy_max,
85	cgz_max,
86	) = state.wkb
87
88	# Set Coriolis parameter.
89	fc = coriolis_frequency * tref	2✔
90
91	# Initialize local arrays.
92	omi_ini = zeros(wave_modes, nxx, nyy, nzz)	8,192✔
93	wnk_ini = zeros(wave_modes, nxx, nyy, nzz)	8,192✔
94	wnl_ini = zeros(wave_modes, nxx, nyy, nzz)	8,192✔
95	wnm_ini = zeros(wave_modes, nxx, nyy, nzz)	8,192✔
96	wad_ini = zeros(wave_modes, nxx, nyy, nzz)	8,192✔
97
98	# Compute initial wavenumbers, intrinsic frequencies and wave-action
99	# densities with initial_wave_field.
100	if wkb_mode != SteadyState()	2✔
101	for k in k0:k1, j in j0:j1, i in i0:i1, alpha in 1:wave_modes	2✔
102	(kdim, ldim, mdim, omegadim, adim) = initial_wave_field(	2,000✔
103	alpha,
104	x[io + i] * lref,
105	y[jo + j] * lref,
106	zc[i, j, k] * lref,
107	)
108	wnk_ini[alpha, i, j, k] = kdim * lref	2,000✔
109	wnl_ini[alpha, i, j, k] = ldim * lref	2,000✔
110	wnm_ini[alpha, i, j, k] = mdim * lref	2,000✔
111	omi_ini[alpha, i, j, k] = omegadim * tref	2,000✔
112	wad_ini[alpha, i, j, k] = adim / rhoref / uref^2 / tref	2,000✔
113	end	2,018✔
114	else
NEW 115	println(	×
116	"Warning: MSGWaM's steady-state mode currently ignores non-orographic initializations!",
117	)
NEW 118	println("")	×
119	end
120
121	# Add orographic wave modes.
122	activate_orographic_source!(	2✔
123	state,
124	omi_ini,
125	wnk_ini,
126	wnl_ini,
127	wnm_ini,
128	wad_ini,
129	)
130
131	# Set initial spectral extents (these will be overwritten in the loop).
132	dk_ini_nd = 0.0	2✔
133	dl_ini_nd = 0.0	2✔
134	dm_ini_nd = 0.0	2✔
135
136	# Set vertical index bounds.
137	kmin = ko == 0 ? k0 - 1 : k0	2✔
138	kmax = k1	2✔
139
140	# Loop over all grid cells with ray volumes.
141	@ivy for k in kmin:kmax, j in j0:j1, i in i0:i1	22✔
142	r = 0	2,200✔
143	s = 0	2,200✔
144
145	# Loop over all ray volumes within a spatial cell.
146	for ix in 1:nrx,	2,200✔
147	ik in 1:nrk,
148	jy in 1:nry,
149	jl in 1:nrl,
150	kz in 1:nrz,
151	km in 1:nrm,
152	alpha in 1:wave_modes
153
154	# Set ray-volume indices.
155	if ko == 0 && k == k0 - 1	2,200✔
156	s += 1	200✔
157
158	# Set surface indices.
159	surface_indices.ixs[s] = ix	200✔
160	surface_indices.jys[s] = jy	200✔
161	surface_indices.kzs[s] = kz	200✔
162	surface_indices.iks[s] = ik	200✔
163	surface_indices.jls[s] = jl	200✔
164	surface_indices.kms[s] = km	200✔
165	surface_indices.alphas[s] = alpha	200✔
166
167	# Set surface ray-volume index.
168	if wad_ini[alpha, i, j, k] == 0.0	200✔
169	surface_indices.rs[s, i, j] = -1	192✔
170	continue	192✔
171	else
172	r += 1	8✔
173	surface_indices.rs[s, i, j] = r	8✔
174	end
175	else
176	if wad_ini[alpha, i, j, k] == 0.0	2,000✔
177	continue	2,000✔
178	end
UNCOV 179	r += 1	×
180	end
181
182	# Set ray-volume positions.
183	rays.x[r, i, j, k] = (x[io + i] - 0.5 * dx + (ix - 0.5) * dx / nrx)	8✔
184	rays.y[r, i, j, k] = (y[jo + j] - 0.5 * dy + (jy - 0.5) * dy / nry)	8✔
185	rays.z[r, i, j, k] = (	8✔
186	zc[i, j, k] - 0.5 * jac[i, j, k] * dz +
187	(kz - 0.5) * jac[i, j, k] * dz / nrz
188	)
189
190	xr = rays.x[r, i, j, k]	8✔
191	yr = rays.y[r, i, j, k]	8✔
192	zr = rays.z[r, i, j, k]	8✔
193
194	# Check if ray volume is too low.
195	if zr < -dz	8✔
196	error(	×
197	"Error in initialize_rays!: Ray volume",
198	r,
199	"at",
200	i,
201	j,
202	k,
203	"is too low!",
204	)
205	end
206
207	# Compute local stratification.
208	n2r = interpolate_stratification(zr, state, N2())	16✔
209
210	# Set spatial extents.
211	rays.dxray[r, i, j, k] = dx / nrx	8✔
212	rays.dyray[r, i, j, k] = dy / nry	8✔
213	rays.dzray[r, i, j, k] = jac[i, j, k] * dz / nrz	8✔
214
215	wnk0 = wnk_ini[alpha, i, j, k]	8✔
216	wnl0 = wnl_ini[alpha, i, j, k]	8✔
217	wnm0 = wnm_ini[alpha, i, j, k]	8✔
218
219	# Ensure correct wavenumber extents.
220	if x_size > 1	8✔
221	dk_ini_nd = dkr_factor * sqrt(wnk0^2 + wnl0^2)	8✔
222	end
223	if y_size > 1	8✔
224	dl_ini_nd = dlr_factor * sqrt(wnk0^2 + wnl0^2)	8✔
225	end
226	if wnm0 == 0.0	8✔
227	error("Error in WKB: wnm0 = 0!")	×
228	else
229	dm_ini_nd = dmr_factor * abs(wnm0)	8✔
230	end
231
232	# Set ray-volume wavenumbers.
233	rays.k[r, i, j, k] =	8✔
234	(wnk0 - 0.5 * dk_ini_nd + (ik - 0.5) * dk_ini_nd / nrk)
235	rays.l[r, i, j, k] =	8✔
236	(wnl0 - 0.5 * dl_ini_nd + (jl - 0.5) * dl_ini_nd / nrl)
237	rays.m[r, i, j, k] =	8✔
238	(wnm0 - 0.5 * dm_ini_nd + (km - 0.5) * dm_ini_nd / nrm)
239
240	# Set spectral extents.
241	rays.dkray[r, i, j, k] = dk_ini_nd / nrk	8✔
242	rays.dlray[r, i, j, k] = dl_ini_nd / nrl	8✔
243	rays.dmray[r, i, j, k] = dm_ini_nd / nrm	8✔
244
245	# Set spectral volume.
246	pspvol = dm_ini_nd	8✔
247	if x_size > 1	8✔
248	pspvol = pspvol * dk_ini_nd	8✔
249	end
250	if y_size > 1	8✔
251	pspvol = pspvol * dl_ini_nd	8✔
252	end
253
254	# Set phase-space wave-action density.
255	rays.dens[r, i, j, k] = wad_ini[alpha, i, j, k] / pspvol	8✔
256
257	# Interpolate winds to ray-volume position.
258	uxr = interpolate_mean_flow(xr, yr, zr, state, U())	8✔
259	vyr = interpolate_mean_flow(xr, yr, zr, state, V())	8✔
260	wzr = interpolate_mean_flow(xr, yr, zr, state, W())	8✔
261
262	wnrk = rays.k[r, i, j, k]	8✔
263	wnrl = rays.l[r, i, j, k]	8✔
264	wnrm = rays.m[r, i, j, k]	8✔
265	wnrh = sqrt(wnrk^2 + wnrl^2)	8✔
266	omir = omi_ini[alpha, i, j, k]	8✔
267
268	# Compute maximum group velocities.
269	cgirx = wnrk * (n2r - omir^2) / (omir * (wnrh^2 + wnrm^2))	8✔
270	if abs(uxr + cgirx) > abs(cgx_max[])	8✔
271	cgx_max[] = abs(uxr + cgirx)	2✔
272	end
273	cgiry = wnrl * (n2r - omir^2) / (omir * (wnrh^2 + wnrm^2))	8✔
274	if abs(vyr + cgiry) > abs(cgy_max[])	8✔
275	cgy_max[] = abs(vyr + cgiry)	3✔
276	end
277	cgirz = -wnrm * (omir^2 - fc^2) / (omir * (wnrh^2 + wnrm^2))	8✔
278	if abs(wzr + cgirz) > abs(cgz_max[i, j, k])	8✔
279	cgz_max[i, j, k] = max(cgz_max[i, j, k], abs(wzr + cgirz))	8✔
280	end
281	end	2,200✔
282
283	# Set ray-volume count.
284	nray[i, j, k] = r	2,200✔
285	if r > nray_wrk	2,200✔
286	error(	×
287	"Error in initialize_rays!: nray",
288	[i, j, k],
289	" > nray_wrk =",
290	nray_wrk,
291	)
292	end
293
294	# Check if surface ray-volume count is correct.
295	if ko == 0 && k == k0 - 1	2,200✔
296	if s != n_sfc	200✔
297	error(	×
298	"Error in initialize_rays!: s =",
299	s,
300	"/= n_sfc =",
301	n_sfc,
302	"at (i, j, k) = ",
303	(i, j, k),
304	)
305	end
306	end
307	end	×
308
309	# Compute global ray-volume count.
310	@ivy local_sum = sum(nray[i0:i1, j0:j1, kmin:kmax])	4✔
311	global_sum = MPI.Allreduce(local_sum, +, comm)	2✔
312
313	# Print information.
314	if master	2✔
315	println("MSGWaM:")	2✔
316	println("Global ray-volume count: ", global_sum)	2✔
317	println("Maximum number of ray volumes per cell: ", nray_max)	2✔
318	println("")	2✔
319	end
320
321	return	2✔
322	end

Atmospheric-Dynamics-GUF / PinCFlow.jl / 18561934460

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous