-
Notifications
You must be signed in to change notification settings - Fork 0
/
run_intro_transport.py
executable file
·232 lines (193 loc) · 12.8 KB
/
run_intro_transport.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
"""Generate samples with configuration specified in 'config.py'."""
from firedrake import *
import numpy as np
import logging
from time import process_time_ns
from functools import partial
from src.discretisation.space import get_space_discretisation_from_CONFIG, SpaceDiscretisation
from src.discretisation.time import TimeDiscretisation
from src.data_dump.setup import update_logfile
from src.algorithms.select import Algorithm, select_algorithm
from src.noise import SamplingStrategy, select_sampling
from src.predefined_data import get_function
from src.string_formatting import format_runtime, format_header
from src.utils import logstring_to_logger
from src.math.distances.space import l2_distance, h1_distance, V_distance, V_sym_distance
from src.math.distances.Bochner_time import linf_X_distance, l2_X_distance, end_time_X_distance, h_minus1_X_distance, w_minus1_inf_X_distance
from src.math.norms.space import l2_space, h1_space, hdiv_space
from src.math.norms.Bochner_time import linf_X_norm, l2_X_norm, end_time_X_norm, h_minus1_X_norm
from src.math.energy import kinetic_energy, potential_energy
from src.postprocess.time_convergence import TimeComparison
from src.postprocess.stability_check import StabilityCheck
from src.postprocess.energy_check import Energy
from src.postprocess.statistics import StatisticsObject
from src.postprocess.processmanager import ProcessManager
#load global and lokal configs
from configs import intro_transport as cf
from configs import intro_global as gcf
def generate_one(time_disc: TimeDiscretisation,
space_disc: SpaceDiscretisation,
initial_condition: Function,
noise_coefficient: Function,
ref_to_time_to_det_forcing: dict[int,dict[float,Function]],
algorithm: Algorithm,
sampling_strategy: SamplingStrategy) -> tuple[dict[int,list[float]], dict[int,dict[float,Function]], dict[int,dict[float,Function]]]:
"""Run the numerical experiment once.
Return noise and solution."""
### Generate noise on all refinement levels
ref_to_noise_increments = sampling_strategy(time_disc.refinement_levels,time_disc.initial_time,time_disc.end_time)
ref_to_time_to_velocity = dict()
ref_to_time_to_velocity_midpoints = dict()
ref_to_time_to_pressure = dict()
ref_to_time_to_pressure_midpoints = dict()
for level in ref_to_noise_increments:
### Solve algebraic system
ref_to_time_to_velocity[level], ref_to_time_to_pressure[level], ref_to_time_to_velocity_midpoints[level], ref_to_time_to_pressure_midpoints[level] = algorithm(
space_disc=space_disc,
time_grid=time_disc.ref_to_time_grid[level],
noise_steps= ref_to_noise_increments[level],
noise_coefficient=noise_coefficient,
initial_condition=initial_condition,
time_to_det_forcing = ref_to_time_to_det_forcing[level],
Reynolds_number=1
)
return ref_to_noise_increments, ref_to_time_to_velocity, ref_to_time_to_pressure, ref_to_time_to_velocity_midpoints, ref_to_time_to_pressure_midpoints
def generate() -> None:
"""Runs the experiment multiple times."""
update_logfile(gcf.DUMP_LOCATION,cf.NAME_LOGFILE_GENERATE)
logging.basicConfig(filename=cf.NAME_LOGFILE_GENERATE,format='%(asctime)s| \t %(message)s', datefmt='%d/%m/%Y %I:%M:%S %p',
level=logstring_to_logger(gcf.LOG_LEVEL),force=True)
# define discretisation
space_disc = get_space_discretisation_from_CONFIG(name_mesh=gcf.MESH_NAME,
space_points=gcf.NUMBER_SPACE_POINTS,
velocity_element=gcf.VELOCITY_ELEMENT,
velocity_degree=gcf.VELOCITY_DEGREE,
pressure_element=gcf.PRESSURE_ELEMENT,
pressure_degree=gcf.PRESSURE_DEGREE,
name_bc=gcf.NAME_BOUNDARY_CONDITION
)
logging.info(space_disc)
time_disc = TimeDiscretisation(initial_time=gcf.INITIAL_TIME, end_time=gcf.END_TIME,refinement_levels=gcf.REFINEMENT_LEVELS)
logging.info(time_disc)
# define data
logging.info(f"\nINITIAL INTENSITY:\t{gcf.INITIAL_INTENSITY}\nNOISE INTENSITY:\t{gcf.NOISE_INTENSITY}")
initial_condition = gcf.INITIAL_INTENSITY*get_function(gcf.INITIAL_CONDITION_NAME,space_disc,gcf.INITIAL_FREQUENZY_X,gcf.INITIAL_FREQUENZY_Y)
noise_coefficient = gcf.NOISE_INTENSITY*get_function(gcf.NOISE_COEFFICIENT_NAME,space_disc,gcf.NOISE_FREQUENZY_X,gcf.NOISE_FREQUENZY_Y)
ref_to_time_to_det_forcing = {level: {time: gcf.FORCING_INTENSITY*get_function(gcf.FORCING,space_disc,gcf.FORCING_FREQUENZY_X,gcf.FORCING_FREQUENZY_Y) for time in time_disc.ref_to_time_grid[level]} for level in time_disc.refinement_levels}
# select algorithm
algorithm = select_algorithm(gcf.MODEL_NAME,cf.ALGORITHM_NAME)
#select sampling
sampling_strategy = select_sampling(gcf.NOISE_INCREMENTS)
#Initialise process managers to handle data processing
if gcf.TIME_CONVERGENCE:
time_convergence_velocity = ProcessManager([
TimeComparison(time_disc.ref_to_time_stepsize,"Linf_L2_velocity",linf_X_distance,l2_distance,gcf.TIME_COMPARISON_TYPE),
TimeComparison(time_disc.ref_to_time_stepsize,"End_time_L2_velocity",end_time_X_distance,l2_distance,gcf.TIME_COMPARISON_TYPE),
TimeComparison(time_disc.ref_to_time_stepsize,"L2_H1_velocity",l2_X_distance,h1_distance,gcf.TIME_COMPARISON_TYPE),
])
time_convergence_pressure = ProcessManager([
TimeComparison(time_disc.ref_to_time_stepsize,"L2_L2_pressure",l2_X_distance,l2_distance,gcf.TIME_COMPARISON_TYPE),
TimeComparison(time_disc.ref_to_time_stepsize,"H-1_L2_pressure",h_minus1_X_distance,l2_distance,gcf.TIME_COMPARISON_TYPE),
TimeComparison(time_disc.ref_to_time_stepsize,"W-1_inf_L2_pressure",w_minus1_inf_X_distance,l2_distance,gcf.TIME_COMPARISON_TYPE)
])
if gcf.STABILITY_CHECK:
stability_check_velocity = ProcessManager([
StabilityCheck(time_disc.ref_to_time_stepsize,"Linf_L2_velocity",linf_X_norm,l2_space),
StabilityCheck(time_disc.ref_to_time_stepsize,"End_time_L2_velocity",end_time_X_norm,l2_space),
StabilityCheck(time_disc.ref_to_time_stepsize,"L2_H1_velocity",l2_X_norm,h1_space),
StabilityCheck(time_disc.ref_to_time_stepsize,"L2_Hdiv_velocity",l2_X_norm,hdiv_space)
])
stability_check_pressure = ProcessManager([
StabilityCheck(time_disc.ref_to_time_stepsize,"L2_L2_pressure",l2_X_norm,l2_space),
StabilityCheck(time_disc.ref_to_time_stepsize,"H-1_L2_pressure",h_minus1_X_norm,l2_space)
])
if gcf.ENERGY_CHECK:
energy_check_velocity = ProcessManager([
Energy(time_disc,"kinetic_energy",kinetic_energy),
#Energy(time_disc,"potential_energy",potential_energy)
])
if gcf.IND_ENERGY_CHECK:
sample_to_energy_check_velocity = dict()
if gcf.STATISTICS_CHECK:
statistics_velocity = StatisticsObject("velocity",time_disc.ref_to_time_grid,space_disc.velocity_space)
statistics_velocity_midpoints = StatisticsObject("velocity_midpoints",time_disc.ref_to_time_grid,space_disc.velocity_space)
statistics_pressure = StatisticsObject("pressure",time_disc.ref_to_time_grid,space_disc.pressure_space)
statistics_pressure_midpoints = StatisticsObject("pressure_midpoints",time_disc.ref_to_time_grid,space_disc.pressure_space)
runtimes = {"solving": 0,"comparison": 0, "stability": 0, "energy": 0, "statistics": 0}
### start MC iteration
print(format_header("START MONTE CARLO ITERATION") + f"\nRequested samples:\t{gcf.MC_SAMPLES}")
new_seeds = range(gcf.MC_SAMPLES)
for k in new_seeds:
### get solution
print(f"{k*100/len(new_seeds):4.2f}% completed")
time_mark = process_time_ns()
ref_to_noise_increments, ref_to_time_to_velocity, ref_to_time_to_pressure, ref_to_time_to_velocity_midpoints, ref_to_time_to_pressure_midpoints = generate_one(time_disc,space_disc,initial_condition,noise_coefficient,ref_to_time_to_det_forcing,algorithm,sampling_strategy)
runtimes["solving"] += process_time_ns()-time_mark
#update data using solution
if gcf.TIME_CONVERGENCE:
time_mark = process_time_ns()
time_to_fine_velocity = ref_to_time_to_velocity[time_disc.refinement_levels[-1]]
time_convergence_velocity.update(ref_to_time_to_velocity,time_to_fine_velocity)
time_to_fine_pressure = ref_to_time_to_pressure[time_disc.refinement_levels[-1]]
time_convergence_pressure.update(ref_to_time_to_pressure,time_to_fine_pressure)
runtimes["comparison"] += process_time_ns()-time_mark
if gcf.STABILITY_CHECK:
time_mark = process_time_ns()
stability_check_velocity.update(ref_to_time_to_velocity)
stability_check_pressure.update(ref_to_time_to_pressure)
runtimes["stability"] += process_time_ns()-time_mark
if gcf.ENERGY_CHECK:
time_mark = process_time_ns()
energy_check_velocity.update(ref_to_time_to_velocity,ref_to_noise_increments)
runtimes["energy"] += process_time_ns()-time_mark
if gcf.IND_ENERGY_CHECK and k <= gcf.IND_ENERGY_NUMBER:
time_mark = process_time_ns()
ind_energy_check_velocity = ProcessManager([
Energy(time_disc,f"ind_kinetic_energy_{k}",kinetic_energy),
#Energy(time_disc,f"ind_potential_energy_{k}",potential_energy)
])
ind_energy_check_velocity.update(ref_to_time_to_velocity,ref_to_noise_increments)
sample_to_energy_check_velocity[k] = ind_energy_check_velocity
runtimes["energy"] += process_time_ns()-time_mark
if gcf.STATISTICS_CHECK:
time_mark = process_time_ns()
statistics_velocity.update(ref_to_time_to_velocity)
statistics_velocity_midpoints.update(ref_to_time_to_velocity_midpoints)
statistics_pressure.update(ref_to_time_to_pressure)
statistics_pressure_midpoints.update(ref_to_time_to_pressure_midpoints)
runtimes["statistics"] += process_time_ns()-time_mark
### storing processed data
if gcf.TIME_CONVERGENCE:
logging.info(format_header("TIME CONVERGENCE") + f"\nComparisons are stored in:\t {cf.TIME_DIRECTORYNAME}/")
logging.info(time_convergence_velocity)
logging.info(time_convergence_pressure)
time_convergence_velocity.save(cf.TIME_DIRECTORYNAME)
time_convergence_pressure.save(cf.TIME_DIRECTORYNAME)
if gcf.STABILITY_CHECK:
logging.info(format_header("STABILITY CHECK") + f"\nStability checks are stored in:\t {cf.STABILITY_DIRECTORYNAME}/")
logging.info(stability_check_velocity)
logging.info(stability_check_pressure)
stability_check_velocity.save(cf.STABILITY_DIRECTORYNAME)
stability_check_pressure.save(cf.STABILITY_DIRECTORYNAME)
if gcf.ENERGY_CHECK:
logging.info(format_header("ENERGY CHECK") + f"\nEnergy checks are stored in:\t {cf.ENERGY_DIRECTORYNAME}/")
energy_check_velocity.save(cf.ENERGY_DIRECTORYNAME)
energy_check_velocity.plot(cf.ENERGY_DIRECTORYNAME)
#energy_check_velocity.plot_individual(cf.ENERGY_DIRECTORYNAME)
if gcf.IND_ENERGY_CHECK:
logging.info(format_header("ENERGY CHECK") + f"\nEnergy checks are stored in:\t {cf.ENERGY_DIRECTORYNAME}/")
for sample in sample_to_energy_check_velocity.keys():
sample_to_energy_check_velocity[sample].save(cf.ENERGY_DIRECTORYNAME + "/individual")
#sample_to_energy_check_velocity[sample].plot(cf.ENERGY_DIRECTORYNAME + "/individual")
energy_check_velocity.save(cf.ENERGY_DIRECTORYNAME)
if gcf.STATISTICS_CHECK:
logging.info(format_header("STATISTICS") + f"\nStatistics are stored in:\t {cf.VTK_DIRECTORY + '/' + cf.STATISTICS_DIRECTORYNAME}/")
statistics_velocity.save(cf.VTK_DIRECTORY + "/" + cf.STATISTICS_DIRECTORYNAME)
statistics_velocity_midpoints.save(cf.VTK_DIRECTORY + "/" + cf.STATISTICS_DIRECTORYNAME)
statistics_pressure.save(cf.VTK_DIRECTORY + "/" + cf.STATISTICS_DIRECTORYNAME)
statistics_pressure_midpoints.save(cf.VTK_DIRECTORY + "/" + cf.STATISTICS_DIRECTORYNAME)
#show runtimes
logging.info(format_runtime(runtimes))
print(f"Logs saved in:\t {cf.NAME_LOGFILE_GENERATE}")
if __name__ == "__main__":
generate()