from typing import Any
import os
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
else:
import numpy as np
from g3d.defs import G3DWorld
from g3d.defs import G3DSizes
from g3d.flow.placeholders import ConsField
from g3d.flow.domain import Domain
from g3d.flow.solver.step import G3DStep
[docs]
class DualTimeStepper(G3DStep):
def __init__(self, *args: Any, **kwargs: Any) -> None:
super().__init__(*args, **kwargs)
self.__r_max = 1.0
self.__r_max_old = 1.0
self.__cons0 = ConsField()
self.__cons_ref = ConsField()
self._cons_incr = ConsField()
self._cons_curr = ConsField()
g3d_sizes = G3DSizes()
self.__q_res = np.zeros(g3d_sizes.n_var)
self.__q_res_ref = np.zeros(g3d_sizes.n_var)
self._cons_incr_pointers = self._cons_incr.get_pointers()
self._cons_keys = self._cons.get_keys()
def _allocate(self) -> None:
super()._allocate()
g3d_sizes = G3DSizes()
n_cells = g3d_sizes.n_cells
self.__cons0.setup(n_cells=n_cells)
self.__cons_ref.setup(n_cells=n_cells)
self._cons_incr.setup(n_cells=n_cells)
self._cons_curr.setup(n_cells=n_cells)
[docs]
def associate(self, dom: Domain) -> None:
"""
:param dom:
"""
super().associate(dom)
g3d_sizes = G3DSizes()
self.__q_res = np.zeros(g3d_sizes.n_var)
self.__q_res_ref = np.zeros(g3d_sizes.n_var)
self._cons_incr_pointers = self._cons_incr.get_pointers()
self._cons_keys = self._cons.get_keys()
def __nullify(self) -> None:
g3d_sizes = G3DSizes()
self.__cons_ref.nullify()
self._cons_incr.nullify()
self.__r_max = 1.0
self.__r_max_old = self.__r_max
self.__q_res[0:g3d_sizes.n_var] = 0.
self.__q_res_ref[0:g3d_sizes.n_var] = 0.
[docs]
def fwd(self) -> None:
"""
"""
g3d_world = G3DWorld()
self.__nullify()
self._cons_curr.set_from_cons_field(self._cons)
self._calculate_time_step()
for sub_iter in range(g3d_world.num.time.n_sub_iter):
self.__r_max_old = self.__r_max
self.__cons0.set_from_cons_field(self._cons_incr)
# RK 1
self._flux_inc(1)
self._cons_incr.WAXPY(1.0, self._dcons, self.__cons0)
self.__cons0.AXPBY(0.5, 0.5, self._cons_incr)
self._domain.auxvr()
# RK 2
self._flux_inc(2)
self._cons_incr.WAXPY(0.5, self._dcons, self.__cons0)
self._domain.auxvr()
# RK 3
self._flux_inc(3)
self._cons_incr.WAXPY(0.5, self._dcons, self.__cons0)
self._calculate_residual(sub_iter)
self._cons.WAXPY(1.0, self._cons_incr, self._cons_curr)
g3d_world.status.time = g3d_world.status.time + g3d_world.num.time.dtime
def _calculate_time_step(self) -> None:
g3d_world = G3DWorld()
self._calculate_s_max()
self._dt[0:g3d_world.sizes.n_cells] = (
g3d_world.num.time.cfl /
np.maximum(1.0, 0.5 * g3d_world.num.time.dtime * self._s_max[0:g3d_world.sizes.n_cells]))
def _calculate_residual(self, sub_iter: int) -> None:
g3d_world = G3DWorld()
self.__cons_ref.AXPY(-1.0, self._cons_incr)
self.__cons_ref.get_max_norm()
self._cons.get_max_norm()
q_max = self._cons.get_extreme_value_array('max_norm')
self.__cons_ref.normalize_extreme_values(np.maximum(q_max, 1.e-20), 'max_norm')
self.__q_res[0:g3d_world.sizes.n_var] = self.__cons_ref.get_extreme_value_array('max_norm')
if sub_iter == 0:
self.__q_res_ref[0:g3d_world.sizes.n_var] = np.where(self.__q_res == 0., 1.0, self.__q_res)
self.__q_res[0:g3d_world.sizes.n_var] /= self.__q_res_ref[0:g3d_world.sizes.n_var]
self.__r_max = max(self.__q_res)
trend = '='
if self.__r_max != self.__r_max_old and sub_iter > 0:
trend = '+'
if self.__r_max < self.__r_max_old:
trend = '-'
print('sub_iter %d/%d: residual = % 10.5E (%s %s)' % (sub_iter + 1, g3d_world.num.time.n_sub_iter,
self.__r_max, self._cons_keys[np.argmax(self.__q_res)],
trend))
# self.__cons_ref.print_extreme_values_by_type('max_norm')
self.__cons_ref.set_from_cons_field(self._cons_incr)
def _flux_inc(self, RK_stage: int = 1) -> None:
# base class DualTimeStepper does not have an associated flux_inc method
pass
def _flux_inc_common(self, RK_stage: int = 1) -> None:
g3d_world = G3DWorld()
idx1 = 0
idx2 = g3d_world.sizes.n_cells
self._domain.auxvr()
self._dcons.evaluate(RK_stage)
for i in range(g3d_world.sizes.n_var):
self._dcons_pointers[i][idx1:idx2] = (self._dcons_pointers[i][idx1:idx2] * g3d_world.num.time.dtime
- self._cons_incr_pointers[i][idx1:idx2] * self._dt[idx1:idx2])