"""flow solver (Solver class definition)"""
from typing import Any
import os
import shutil
import re
import numpy
import cProfile
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
from cupy import ndarray
else:
import numpy as np
from numpy import ndarray
import matplotlib.pyplot as plt
from g3d.defs import G3DWorld
from g3d.defs import TimeStepperMode
from g3d.defs import clean_path
from g3d.flow.io import G3DIO
from g3d.flow.solver import MOC
from g3d.flow.solver import G3DStep
from g3d.flow.solver import RungeKutta
from g3d.flow.solver import BackwardEuler
from g3d.flow.solver import CrankNicolson
from g3d.flow.solver import PMin
from g3d.flow.solver import PMax
from g3d.gas import Gas
from g3d.flow.placeholders import ConsField, AuxvField, PostField, VGradField
from g3d.flow.domain import Domain
from g3d.flow.geo import Mesh
from g3d.flow.domain import TurbulenceModel
from g3d.defs import G3DEntity
[docs]
class Solver(G3DEntity):
def __init__(self, *args: Any, **kwargs: Any) -> None:
super().__init__(*args, **kwargs)
self.__first_run = True
self.__gas = Gas()
self.__turb = TurbulenceModel()
self.__residual_log = []
self.__moc = MOC()
self.__ro_ref = np.zeros(0)
self.__step = G3DStep()
self.__io = G3DIO()
self.__limiters = []
self.__nstep_next_save = 0
self.__nlog = 0
self.__base_name = ''
self.__curiter_name = ''
self.__res_log_name = ''
self.__domain = Domain()
self.__turb = TurbulenceModel()
self.__mesh = Mesh()
self.__cons = ConsField()
self.__auxv = AuxvField()
self.__post = PostField()
self.__vgrad = VGradField()
self.__ro = np.zeros(0)
self.__allocate()
def __allocate(self) -> None:
g3d_world = G3DWorld()
if g3d_world.num.limiters.p_min is not None:
limiter = PMin()
self.__limiters.append(limiter)
if g3d_world.num.limiters.p_max is not None:
limiter = PMax()
self.__limiters.append(limiter)
self.__ro_ref = np.zeros(g3d_world.sizes.n_cells)
[docs]
def setup_and_allocate(self, dom: Domain) -> None:
"""
:param dom:
"""
g3d_world = G3DWorld()
g3d_world.sizes.n_faces = dom.get_max_FMP_size()
self.__allocate()
if g3d_world.num.time.time_stepper_mode is TimeStepperMode.EXPLICIT_RUNGE_KUTTA:
self.__step = RungeKutta()
elif g3d_world.num.time.time_stepper_mode is TimeStepperMode.BACKWARD_EULER_DUAL_TIME_STEPPER:
self.__step = BackwardEuler()
elif g3d_world.num.time.time_stepper_mode is TimeStepperMode.CRANK_NICOLSON_DUAL_TIME_STEPPER:
self.__step = CrankNicolson()
self.__step.setup_and_allocate()
[docs]
def flux_test(self) -> None: # test function
"""
"""
self.__step.flux_test()
[docs]
def flux_check(self) -> None:
"""
"""
self.__step.flux_check()
[docs]
def diff_flux_check(self) -> None:
"""
"""
self.__step.diff_flux_check()
[docs]
def ltsp_test(self) -> None: # test function
"""
"""
self.__step.ltsp_test()
[docs]
def step_test(self) -> None: # test function
"""
"""
self.__domain.status_check()
self.__step.fwd()
[docs]
def associate(self, dom: Domain) -> None:
"""
:param dom:
"""
self.__domain = dom
self.__gas = dom.gas
self.__turb = dom.turb
self.__mesh = dom.mesh
self.__cons = dom.cons
self.__auxv = dom.auxv
self.__post = dom.post
self.__vgrad = dom.vgrad
self.__step.associate(dom)
self.__ro = dom.get_cons_array('ro')
[docs]
def change_eqn_mode(self) -> None:
"""
"""
self.__step.change_eqn_mode()
[docs]
def get_res_log(self) -> list:
"""
:return:
"""
return self.__residual_log
[docs]
def reset_solver_status(self) -> None:
"""
"""
g3d_world = G3DWorld()
g3d_world.io.prev_solution_folder = ''
g3d_world.status.curiter = 0
g3d_world.status.istep_save = 0
g3d_world.status.time = 0.
g3d_world.status.residual = 1.e20
self.__first_run = True
def __get_filenames(self) -> None:
g3d_world = G3DWorld()
self.__base_name = clean_path(g3d_world.io.wd + '/' + g3d_world.io.jobname)
self.__curiter_name = self.__base_name + '.curiter'
self.__res_log_name = self.__base_name + '.log'
[docs]
def import_flow_data(self, *args: Any, **kwargs: Any) -> None:
"""
:param args:
:param kwargs:
"""
sol_dir = 'curiter'
if len(args) > 0:
if type(args[0]) is str:
sol_dir = args[0]
for key, value in kwargs.items():
if key == 'sol_dir':
sol_dir = value
if sol_dir == 'curiter':
self.__get_filenames()
sol_dir = self.__read_curiter_file()
if len(sol_dir) != 0:
self.__cons.import_data(sol_dir)
[docs]
def import_g3d_flow_data(self, *args: Any, **kwargs: Any) -> None:
"""
:param args:
:param kwargs:
"""
file_name = ''
if len(args) > 0:
if type(args[0]) is str:
file_name = args[0]
for key, value in kwargs.items():
if key == 'file_name':
file_name = value
if len(file_name) > 0:
with open(file_name, 'r+b') as f:
header = f.readline().decode('utf8').strip()
if header == 'flow data g3dflow binary format':
print(header)
header = f.readline().decode('utf8').strip()
if header == '#specified variables':
print(header)
n_var = int(f.readline().decode('utf8').strip())
print('%d variables:' % n_var)
variables = []
for ivar in range(n_var):
var_name = f.readline().decode('utf8').strip()
print('variable %d - %s' % (ivar + 1, var_name))
variables.append(var_name)
self.__cons.import_var_data_from_binary(f, variables)
def __status_report(self) -> None:
g3d_world = G3DWorld()
g3d_world.display.hline()
print(' solver status report:')
g3d_world.display.hline()
print('')
if g3d_world.num.time.time_mode.value > 1 or g3d_world.num.time.time_stepper_mode.value > 1:
print(' simulated time: %10.5E [s]' % g3d_world.status.time)
print(' iteration number: %d' % g3d_world.status.curiter)
print(' continuity residual: %10.5E' % g3d_world.status.residual)
print('')
print(' flow field data:')
print('')
self.__domain.auxvr()
self.__domain.print_flow_field_extreme_values()
if g3d_world.options.use_velocity_gradients:
print()
self.__domain.print_velocity_gradients_extreme_values()
print('')
g3d_world.display.hline()
def __read_curiter_file(self) -> str:
g3d_world = G3DWorld()
regex_curiter = re.compile(
r'\s*curiter:\s+(?P<curiter>\d+)'
)
regex_istep_save = re.compile(
r'\s*istep_save:\s+(?P<istep>\d+)'
)
regex_sol_dir = re.compile(
r'\s*sol_*dir:\s+(?P<sol_dir>[\w/\-.]+)'
)
regex_time = re.compile(
r'\s*time:\s+(?P<time>[+-]*\d+\.*\d+[eE]*[+-]*\d*)'
)
regex_residual = re.compile(
r'\s*residual:\s+(?P<res>[+-]*\d+\.*\d+[eE]*[+-]*\d*)'
)
sol_dir = ''
with open(self.__curiter_name, 'r') as f:
line = f.readline()
while line:
match = regex_curiter.search(line)
if match:
g3d_world.status.curiter = int(match.group('curiter'))
match = regex_istep_save.search(line)
if match:
g3d_world.status.istep_save = int(match.group('istep'))
match = regex_sol_dir.search(line)
if match:
sol_dir = match.group('sol_dir')
g3d_world.io.prev_solution_folder = os.path.dirname(sol_dir)
match = regex_time.search(line)
if match:
g3d_world.status.time = float(match.group('time'))
match = regex_residual.search(line)
if match:
g3d_world.status.residual = float(match.group('res'))
line = f.readline()
return sol_dir
def __solver_restart(self) -> None:
g3d_world = G3DWorld()
regex_sol_dir = re.compile(
r'\s*sol_*dir:\s+(?P<sol_dir>[\w/\-.]+)'
)
if os.path.isfile(self.__curiter_name):
if g3d_world.options.restart:
# restart from the most recent solution
sol_dir = self.__read_curiter_file()
if len(sol_dir) != 0:
self.__cons.import_data(sol_dir)
self.__status_report()
restart = g3d_world.status.curiter < g3d_world.options.nstep
if g3d_world.options.use_residual_break and restart:
restart = g3d_world.options.break_residual < g3d_world.status.residual
if restart:
self.__update_reference_density_field()
else:
self.__post.update(self.__cons, self.__auxv, self.__gas)
else:
# clean up working directory and start from scratch
with open(self.__curiter_name, 'r') as f:
line = f.readline()
while line:
match = regex_sol_dir.search(line)
if match:
shutil.rmtree(os.path.dirname(match.group('sol_dir')))
line = f.readline()
if os.path.isfile(self.__curiter_name):
os.remove(self.__curiter_name)
if os.path.isfile(self.__res_log_name):
os.remove(self.__res_log_name)
[docs]
def run(self) -> None:
"""
:return:
"""
g3d_world = G3DWorld()
if g3d_world.options.use_irsm and g3d_world.num.time.cfl > g3d_world.num.irsm.cfl_crit:
self.__domain.generate_irsm_matrix()
self.__get_filenames()
self.__update_reference_density_field()
if self.__first_run:
self.__first_run = False
self.__solver_restart()
if g3d_world.options.nstep_save == -1:
g3d_world.options.nstep_save = g3d_world.options.nstep
if g3d_world.status.curiter >= g3d_world.options.nstep:
print(' max number of iterations reached')
return
if g3d_world.options.use_residual_break and g3d_world.options.break_residual > g3d_world.status.residual:
print(' break residual reached')
return
self.__residual_log = np.zeros(min(g3d_world.options.nstep - g3d_world.status.curiter,
g3d_world.options.nstep_save - g3d_world.status.istep_save))
self.__step.initialize()
g3d_world.display.hline()
print(' start iterating ...')
g3d_world.display.hline()
self.__nlog = 0
self.__nstep_next_save = min(g3d_world.options.nstep_save, g3d_world.options.nstep - g3d_world.status.curiter)
if g3d_world.options.attach_profiler:
p = cProfile.Profile()
p.runcall(self.__iterate)
if g3d_world.display.verbose:
p.print_stats()
p.dump_stats(self.__base_name + '.pstat')
else:
self.__iterate()
def __iterate(self) -> None:
g3d_world = G3DWorld()
for i in range(g3d_world.status.curiter, g3d_world.options.nstep):
g3d_world.status.istep = i + 1
self.__domain.BC_mod()
self.__step.fwd()
for limiter in self.__limiters:
limiter.evaluate(self.__domain)
self.__calculate_continuity_residual()
self.__residual_log[self.__nlog] = g3d_world.status.residual
self.__nlog += 1
if g3d_world.display.verbose or g3d_world.display.print_iterations:
print(' iteration: %7d, residual: %10.5E' % (i + 1, g3d_world.status.residual))
g3d_world.status.istep_save = g3d_world.status.istep_save + 1
if g3d_world.display.progress_bar:
progress = np.int64(np.floor(g3d_world.status.istep_save / self.__nstep_next_save * 100.))
print('|' + '*' * progress + ' ' * (100 - progress) + '|', end='\r')
if (g3d_world.status.istep_save == g3d_world.options.nstep_save or
i == (g3d_world.options.nstep - 1) or
(g3d_world.options.use_residual_break and
g3d_world.options.break_residual > g3d_world.status.residual)):
if g3d_world.status.istep_save == g3d_world.options.nstep_save:
g3d_world.status.istep_save = 0
prev_curiter = g3d_world.status.curiter
g3d_world.status.curiter = i + 1
self.__nstep_next_save = min(g3d_world.options.nstep_save,
g3d_world.options.nstep - g3d_world.status.curiter)
self.__status_report()
self.__post.update(self.__cons, self.__auxv, self.__gas)
sol_file = clean_path(self.__base_name + '.' + str(i + 1) + '/' + g3d_world.io.jobname)
print('exporting flow data to binary files: %s' % (sol_file + '_*.flow'))
self.__io.export_flow(self.__domain, self.__cons, self.__auxv, self.__post,
name=sol_file, xmlname=self.__base_name)
if g3d_world.options.use_velocity_gradients:
self.__io.export_velocity_gradients(self.__domain, self.__vgrad, name=sol_file,
xmlname=self.__base_name)
if not g3d_world.io.save_all_files and len(g3d_world.io.prev_solution_folder) != 0:
shutil.rmtree(g3d_world.io.prev_solution_folder)
g3d_world.io.prev_solution_folder = os.path.dirname(sol_file)
with open(self.__curiter_name, 'w') as f:
f.write('curiter: %d\n' % g3d_world.status.curiter)
f.write('istep_save: %d\n' % g3d_world.status.istep_save)
f.write('time: %e\n' % g3d_world.status.time)
f.write('residual: %e\n' % g3d_world.status.residual)
f.write('sol_dir: %s\n' % clean_path(sol_file))
if prev_curiter == 0:
with open(self.__res_log_name, 'w') as f:
for tmp_res in self.__residual_log:
f.write('%10.5e\n' % tmp_res)
else:
with open(self.__res_log_name, 'a+') as f:
for tmp_res in self.__residual_log:
f.write('%10.5e\n' % tmp_res)
if g3d_world.options.auto_plot_residuals:
self.__plot_residual()
if (g3d_world.options.use_residual_break and
g3d_world.options.break_residual > g3d_world.status.residual):
break
if i < g3d_world.options.nstep:
self.__residual_log = np.zeros(min(g3d_world.options.nstep - g3d_world.status.curiter,
g3d_world.options.nstep_save - g3d_world.status.istep_save))
self.__nlog = 0
def __update_reference_density_field(self) -> None:
g3d_world = G3DWorld()
self.__ro_ref[0:g3d_world.sizes.n_cells] = self.__ro[0:g3d_world.sizes.n_cells]
def __calculate_continuity_residual(self) -> None:
g3d_world = G3DWorld()
g3d_world.status.residual = np.max(np.abs(self.__ro[0:g3d_world.sizes.n_cells] -
self.__ro_ref[0:g3d_world.sizes.n_cells]))
self.__update_reference_density_field()
def __plot_residual(self, save_plot: bool = True) -> None:
g3d_world = G3DWorld()
res = numpy.genfromtxt(self.__res_log_name, delimiter=",")
res_fig = plt.figure()
plt.semilogy(res)
plt.title('g3dflow.py: ' + g3d_world.io.jobname + ' residual log')
plt.xlabel('iteration')
plt.ylabel('continuity residual')
plt.grid(which='both')
if save_plot:
plt.savefig(g3d_world.io.wd + '/' + g3d_world.io.jobname + '_residual_log.png')
plt.close(res_fig)
else:
plt.show()
[docs]
def get_shock_capturing_field(self) -> ndarray:
"""
:return:
"""
return self.__step.get_shock_capturing_field()