PyMPDATA_examples.Arabas_and_Farhat_2020.analysis_figures_2_and_3 

 1import numpy as np
 2from joblib import Parallel, delayed, parallel_backend
 3from PyMPDATA_examples.Arabas_and_Farhat_2020.setup1_european_corridor import Settings
 4from PyMPDATA_examples.Arabas_and_Farhat_2020.simulation import Simulation
 5from PyMPDATA_examples.utils.error_norms import L2
 6
 7
 8def compute(simulation):
 9    output = []
10    for n_iters in (1, 2):
11        simulation.run(n_iters)
12        output.append(
13            {
14                "n_iters": n_iters,
15                "log2_C": np.log2(simulation.C),
16                "log2_C_opt": np.log2(simulation.settings.C_opt),
17                "log2_l2": np.log2(simulation.l2),
18                "log2_l2_opt": np.log2(simulation.settings.l2_opt),
19                "err2": error_L2_norm(
20                    simulation.solvers,
21                    simulation.settings,
22                    simulation.S,
23                    simulation.nt,
24                    n_iters,
25                ),
26            }
27        )
28    return output
29
30
31def convergence_in_space(num=8):
32    with parallel_backend("threading", n_jobs=-2):
33        data = Parallel(verbose=10)(
34            delayed(compute)(
35                Simulation(Settings(l2_opt=2**log2_l2_opt, C_opt=2**log2_C_opt))
36            )
37            for log2_C_opt in np.linspace(-9.5, -6, num=num)
38            for log2_l2_opt in range(1, 4)
39        )
40        result = {}
41        for pair in data:
42            for datum in pair:
43                label = f" $\\lambda^2\\approx2^{{{datum['log2_l2_opt']}}}$"
44                key = ("upwind" + label, "MPDATA" + label)[datum["n_iters"] - 1]
45                if key not in result:
46                    result[key] = ([], [])
47                result[key][0].append(datum["log2_C"])
48                result[key][1].append(datum["err2"])
49        return result
50
51
52def convergence_in_time(num=13):
53    with parallel_backend("threading", n_jobs=-2):
54        data = Parallel(verbose=10)(
55            delayed(compute)(
56                Simulation(Settings(l2_opt=2**log2_l2_opt, C_opt=2**log2_C_opt))
57            )
58            for log2_C_opt in np.log2((0.01, 0.005, 0.0025))
59            for log2_l2_opt in np.linspace(1.1, 3.5, num=num)
60        )
61        result = {}
62        for pair in data:
63            for datum in pair:
64                label = f" $C\\approx{2**(datum['log2_C_opt']):.4f}$"
65                key = ("upwind" + label, "MPDATA" + label)[datum["n_iters"] - 1]
66                if key not in result:
67                    result[key] = ([], [])
68                result[key][0].append(datum["log2_l2"])
69                result[key][1].append(datum["err2"])
70        return result
71
72
73def error_L2_norm(solvers, settings, S, nt, n_iters: int):
74    numerical = solvers[n_iters].advectee.get()
75    analytical = settings.analytical_solution(S)
76    return L2(numerical, analytical, nt)
def compute(simulation): 
 9def compute(simulation):
10    output = []
11    for n_iters in (1, 2):
12        simulation.run(n_iters)
13        output.append(
14            {
15                "n_iters": n_iters,
16                "log2_C": np.log2(simulation.C),
17                "log2_C_opt": np.log2(simulation.settings.C_opt),
18                "log2_l2": np.log2(simulation.l2),
19                "log2_l2_opt": np.log2(simulation.settings.l2_opt),
20                "err2": error_L2_norm(
21                    simulation.solvers,
22                    simulation.settings,
23                    simulation.S,
24                    simulation.nt,
25                    n_iters,
26                ),
27            }
28        )
29    return output
def convergence_in_space(num=8): 
32def convergence_in_space(num=8):
33    with parallel_backend("threading", n_jobs=-2):
34        data = Parallel(verbose=10)(
35            delayed(compute)(
36                Simulation(Settings(l2_opt=2**log2_l2_opt, C_opt=2**log2_C_opt))
37            )
38            for log2_C_opt in np.linspace(-9.5, -6, num=num)
39            for log2_l2_opt in range(1, 4)
40        )
41        result = {}
42        for pair in data:
43            for datum in pair:
44                label = f" $\\lambda^2\\approx2^{{{datum['log2_l2_opt']}}}$"
45                key = ("upwind" + label, "MPDATA" + label)[datum["n_iters"] - 1]
46                if key not in result:
47                    result[key] = ([], [])
48                result[key][0].append(datum["log2_C"])
49                result[key][1].append(datum["err2"])
50        return result
def convergence_in_time(num=13): 
53def convergence_in_time(num=13):
54    with parallel_backend("threading", n_jobs=-2):
55        data = Parallel(verbose=10)(
56            delayed(compute)(
57                Simulation(Settings(l2_opt=2**log2_l2_opt, C_opt=2**log2_C_opt))
58            )
59            for log2_C_opt in np.log2((0.01, 0.005, 0.0025))
60            for log2_l2_opt in np.linspace(1.1, 3.5, num=num)
61        )
62        result = {}
63        for pair in data:
64            for datum in pair:
65                label = f" $C\\approx{2**(datum['log2_C_opt']):.4f}$"
66                key = ("upwind" + label, "MPDATA" + label)[datum["n_iters"] - 1]
67                if key not in result:
68                    result[key] = ([], [])
69                result[key][0].append(datum["log2_l2"])
70                result[key][1].append(datum["err2"])
71        return result
def error_L2_norm(solvers, settings, S, nt, n_iters: int): 
74def error_L2_norm(solvers, settings, S, nt, n_iters: int):
75    numerical = solvers[n_iters].advectee.get()
76    analytical = settings.analytical_solution(S)
77    return L2(numerical, analytical, nt)