particle-simulator/main.py at master · ddkx21/particle-simulator · GitHub

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
import time
import numpy as np
from datetime import datetime
from collision_detector import SpatialHashCollisionDetector
from force_calculator import *
from particle_generator import *
from particle_state import *
from post_processor import *
from solution import *
from solver import *

import taichi as ti


n_of_threads = 16
ti.init(arch=ti.cpu, cpu_max_num_threads=n_of_threads, default_fp=ti.f64)  # Параллельный расчёт на CPU

# import cProfile
# import pstats
# import tracemalloc

def main():
    """
    Главная функция, запускающая моделирование.

    Загружает начальные условия из файла Excel или создаёт новые,
    создает необходимые объекты,
    запускает решатель задачи,
    сохраняет результаты в файл.
    """
    print('\n'*3)

    # tracemalloc.start()

    # Параметры симуляции ___________________________________________________________________________________
    radii_range = np.array([2.5e-6, 7.5e-6])
    num_particles = 1000
    water_volume_content = 0.02

    # Вычисляем размер области моделирования (см заметки, там формула)
    box_size = np.cbrt(  (np.pi*num_particles*np.sum(radii_range)*np.sum(np.square(radii_range)))  /  (3*water_volume_content)  )
    coord_range = (0, box_size)

    t_stop = 100  # Время окончания симуляции
    dt = 0.04  # Временной шаг для интеграции

    real_time_visualization = False

    should_use_saved_initial_data = True
    filenumber = 0
    filename_to_save = f"results/droplets_N{num_particles}_vol{water_volume_content}_0.xlsx"
    filename_to_load = f"results/droplets_N{num_particles}_vol{water_volume_content}_{filenumber}.xlsx"

    # Физические параметры
    eps_oil = 2.85
    eta_oil = 0.065
    eta_water = 0.001
    rho_water = 1000
    rho_oil = 900
    E = 3e5

    # Граничные условия:
    #   "periodic" - minimum image convention + оборачивание позиций (ОБЯЗАТЕЛЬНО при use_periodic_correction=True)
    #   "open"     - открытая коробка, капли могут покидать границы
    boundary_mode = "periodic"

    # Периодическая поправка из COMSOL (псевдо-периодичность)
    # При True: конвекция = прямой стоклет + поправка от периодических образов = полная периодическая функция Грина
    # Требует boundary_mode = "periodic". Параметры в periodic_correction/data/comsol_params.json
    use_periodic_correction = True

    # ___________________________________________________________________________________

    initial_state = None
    if should_use_saved_initial_data:
        initial_state = DropletState(filename=filename_to_load)
    else:
        particle_generator = UniformDropletGenerator(coord_range=coord_range, radii_range=radii_range, num_particles=num_particles, minimum_distance=1e-6)
        initial_state = particle_generator.generate()
        initial_state.export_to_xlsx(filename_to_save)
        #particle_generator.plot()
        particle_generator.print()
        time.sleep(2)

    num_particles = len(initial_state.radii)


    # Загрузка периодической поправки COMSOL (если включена)
    lattice_correction = None
    if use_periodic_correction:
        if boundary_mode != "periodic":
            raise ValueError("use_periodic_correction=True требует boundary_mode='periodic'")
        from periodic_correction import COMSOLLatticeCorrection
        lattice_correction = COMSOLLatticeCorrection.load_default()

    # Создаем необходимые объекты
    corr_res = lattice_correction.grid_resolution if lattice_correction else 0
    force_calculator = DirectDropletForceCalculator(num_particles=num_particles, eps_oil=eps_oil, eta_oil=eta_oil, eta_water=eta_water, rho_water=rho_water, rho_oil=rho_oil, E=E, L=box_size, boundary_mode=boundary_mode, correction_grid_resolution=corr_res)

    if lattice_correction:
        force_calculator.load_periodic_correction(lattice_correction, L_sim=box_size)

    collision_detector = SpatialHashCollisionDetector(num_particles=num_particles, L=box_size, boundary_mode=boundary_mode)

    solution = DropletSolution(initial_droplet_state=initial_state, real_time_visualization=real_time_visualization)  # Объект для хранения траекторий

    post_processor = DropletPostProcessor(solution, box_size=box_size)  # Объект для визуализации

    solver = EulerDropletSolver(force_calculator=force_calculator, solution=solution, post_processor=post_processor, collision_detector=collision_detector)  # Основной решатель задачи

    # Запуск решения
    # profiler = cProfile.Profile()
    # profiler.enable()
    solver.solve(dt, t_stop)
    # profiler.disable()
    # stats = pstats.Stats(profiler).sort_stats('tottime')
    # stats.dump_stats('results/profile.prof')
    # print("\n Function call statistics:")
    # stats.print_stats(20)  # Вывести топ-20 затратных функций


    # Снимок текущего состояния памяти
    # snapshot = tracemalloc.take_snapshot()
    # Анализируем топ-10 самых больших источников потребления памяти
    # top_stats = snapshot.statistics('lineno')
    # print("[ Top 10 memory consumers ]")
    # for stat in top_stats[:10]:
    #    print(stat)


    # Сохраняем решение в файл
    #stamp = datetime.now().strftime("%Y%m%d_%H%M%S")
    stamp = f'N{num_particles}_vol{water_volume_content}_dt{dt}_{filenumber}_d1'
    results_filename = f"results/results_{stamp}.npz"
    solution.save_chain_to_file(results_filename, precision='float32')

if __name__ == "__main__":
    main()