Numerical solution of partial differential equations using the example of an acoustic wave

In this example, we will consider the numerical solution of a partial differential equation using the example of an acoustic wave oscillation.
The topic of embedded functions will be touched upon to show how you can improve the readability of the code without reducing the speed of its execution.

Let's consider the solution of the string oscillation equation without delving into the code, but only paying attention to its structure.
(The code will be explained in the next chapter).

@views function acoustic_wave_1D_orig()
# physical parameters
    lx =
20.0 r, β = 1.0, 1.0
# counts
    nx = 200
    nvis = 2
    # differentials
    dx = lx / nx
    dt = dx*√(ρ / β)
    nt = nx^2 ÷ 100
    xc = LinRange(dx / 2, lx - dx / 2, nx)
    # initializing arrays
    Pr = @. exp(-(xc - lx/4)^2)
Sample = copy(Pr)
Vx = zeros(Float64, nx - 1)
# algorithm cycle
    @gif for it = 1:nt
            Vx .-= dt ./ ρ .* diff(Pr) ./ dx
            Pr[2:end-1] .-= dt ./ β .* diff(Vx) ./ dx
            plot(xc, [Sample, Pr]; xlims=(0, lx), ylims=(-1.1, 1.1),
                 xlabel="lx", ylabel="Concentration",
                 title="Time = $(round(it*dt,digits=1))")
    end every nvis
end

Несмотря на то, что код сопровождается комментариями, можно явно выделить этапы работы функции.
Может возникнуть предположение, что внедрение функций в код под каждый этап замедлит выполнение главной функции.

Действительно, вызов функции занимает значительное время, в сравнении с арифметическими операциями. Однако компиляторы умеют оптимизировать код и делать функции встраиваемыми.

Чтобы проверить, происходит ли вызов функции - воспользуемся макросом @code_llvm, который покажет код, близкий к ассемблеру. (Для большей читабельности можно включить тёмную тему в Engee)

В строчке %3 = add i64 %1, %0 происходит сложение чисел a и b.
Для того чтобы явно указать компилятору, что не нужно делать оптимизацию встраивания, воспользуемся макросом @noinline

Теперь же происходит вызов функции fa, которая складывает a и b:

%3 = call i64 @j_fa_7789(i64 signext %0, i64 signext %1)

Для того чтобы убедиться в том, что встроенная функция работает быстрее, чем вызываемая - напишем следующий код:

И подключим библиотеку BenchmarkTools для оценки времени выполнения кода:

Как видим, встроенная функция работает, быстрее.

Уравнение акустической волны

Уравнение акустической волны имеет вид:

где - давление, - время, а - положение частиц в некоторой среде. - скорость звука в этой среде.

Введём некоторые обозначения:

- объёмный модуль сжимаемости (βulk modulus)

- плотность среды.
Тогда исходное уравнение второго порядка:

можно переписать как систему двух уравнений первого порядка:

Numerical solution of the acoustic wave equation

Let's introduce a function defining the physical parameters of our equation:

We have set it in such a way as to demonstrate the following scenarios:

Note that a tuple is a specific data type that should be understood when to apply.

However, for a small set of data that is worth unpacking, a tuple is one of the best options.
In the function below foo we can see that the compiler has seen that "unpacking" a tuple is just an assignment. a,b,c values x,y and z. And you don't need to make any calls to the elements of the tuple.

If only physical_parameters I didn't return the motorcade (lx, ρ, β), and the vector [lx, ρ, β] then the assembly code would be much larger and take longer to execute.

The finite element method

To solve the system of equations, we will use the [finite difference] method (https://en.wikipedia.org/wiki/Finite_difference_method ) in conjunction with the Euler method.
(Often, to solve such problems, [apply] (https://en.wikipedia.org/wiki/FTCS_scheme ) more accurate [central difference](https://ru.wikipedia.org/wiki/Finite differences))

The essence of the method is to represent a partial differential equation in terms of small but finite differences.:

Selecting the sampling step

Having the length of the "string" lx we can obtain a spatial discretization by dividing the length by the number of partitions, which we determine ourselves.

• nx - set by the user
• dx = lx / nx - spatial sampling step
• - the time step of the sampling. His choice affects whether the solution will be sustainable. The values of this step are selected based on Courant conditions. In our case, the speed of sound in the medium acts as a constant.
• nt It can be set by the user, or it can depend on the length, so that at least 1 period always passes.
• xc - a set of grid centers along which spatial sampling is performed

Initializing arrays

Let's set arrays corresponding to pressure and velocity.

We will also save the initial pressure array for subsequent visualization.

Before proceeding to the solution, we'll connect a library for visualizing graphs of solutions and animating them.:

Finite difference solution

a = a - 1 ~ a -= 1

Let's rewrite it

how

- a scalar showing the value at a given time at the point

For more information, see the picture below.

Similarly, for pressure P.

Then the final code will look like this:

Let's compare the execution time with our source code.

From which we can conclude that wrapping small code in a function does not significantly affect the speed of its execution. However, the following advantages of dividing the code into functions can be noted:

1. Documentation of functions. Documenting functions is a much more effective and expressive tool than commenting on code.
2. Some functions may be used by you in different scenarios.
3. Creating functions allows you to decompose complex tasks into simple ones, as well as gradually increase the complexity and accuracy of your program. (First, we will insert a simple function, but we will check its operability, and only then we will select more accurate and effective ways to implement it)