diff --git a/config.py b/config.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Y29uZmlnLnB5..a02e313946d37e785ae750c7e61647ed861c80d0_Y29uZmlnLnB5 100644 --- a/config.py +++ b/config.py @@ -11,8 +11,6 @@ """ -from __future__ import print_function - import os import sys import subprocess @@ -28,10 +26,9 @@ except ImportError: from distutils.command.build_ext import build_ext from setuptools import Extension -try: # python 3 - from configparser import ConfigParser -except ImportError: # python 2.7 - from ConfigParser import ConfigParser + +from configparser import ConfigParser + try: import colorlog as logging @@ -222,19 +219,6 @@ def get_default_config(): """Generate default configuration.""" config = ConfigParser() - section = "exclude_pythran" - config.add_section(section) - exclude_pythran = ( - "fluidsim.solvers.plate2d", - "fluidsim.solvers.plate2d.output", - "fluidsim.solvers.sw1l", - "fluidsim.solvers.sw1l.output", - "fluidsim.solvers.ns2d.strat", - "fluidsim.solvers.ns2d.bouss", - ) - for excluded in exclude_pythran: - config.set(section, excluded, "True") - config.add_section("environ") return config diff --git a/fluidsim/__init__.py b/fluidsim/__init__.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vX19pbml0X18ucHk=..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vX19pbml0X18ucHk= 100644 --- a/fluidsim/__init__.py +++ b/fluidsim/__init__.py @@ -24,6 +24,31 @@ """ from pathlib import Path +import os +import sys + +if any( + any(test_tool in arg for arg in sys.argv) + for test_tool in ("pytest", "unittest") +): + print( + "Fluidsim guesses that it is tested so it" + " loads the Agg Matplotlib backend." + ) + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + def _show(*args, **kwargs): + pass + + plt.show = _show + + if "FLUID_COMPILE_CACHEDJIT" not in os.environ: + from fluidpythran import set_compile_cachedjit + + set_compile_cachedjit(False) from ._version import __version__, get_local_version diff --git a/fluidsim/base/output/increments.py b/fluidsim/base/output/increments.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvaW5jcmVtZW50cy5weQ==..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvaW5jcmVtZW50cy5weQ== 100644 --- a/fluidsim/base/output/increments.py +++ b/fluidsim/base/output/increments.py @@ -1,7 +1,4 @@ -from __future__ import division, print_function - -from builtins import range import h5py import os import numpy as np @@ -4,7 +1,10 @@ import h5py import os import numpy as np +# pythran import numpy as np + +from fluidpythran import pythran_def, Array from fluiddyn.util import mpi from .base import SpecificOutput @@ -8,7 +8,28 @@ from fluiddyn.util import mpi from .base import SpecificOutput -from .util_pythran import strfunc_from_pdf + +A = Array[float, "2d"] +# pythran def strfunc_from_pdf( +# int32[], float64[][], float64[][], float, bool) + + +@pythran_def +def strfunc_from_pdf( + rxs: A, pdf: A, values: A, order: float, absolute: bool = False +): + """Compute structure function of specified order from pdf for increments + module. + + """ + S_order = np.empty(rxs.shape) + if absolute: + values = abs(values) + for irx in range(rxs.size): + deltainc = abs(values[irx, 1] - values[irx, 0]) + S_order[irx] = deltainc * np.sum(pdf[irx] * values[irx] ** order) + + return S_order class Increments(SpecificOutput): diff --git a/fluidsim/base/output/util_pythran.py b/fluidsim/base/output/util_pythran.py deleted file mode 100644 index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vYmFzZS9vdXRwdXQvdXRpbF9weXRocmFuLnB5..0000000000000000000000000000000000000000 --- a/fluidsim/base/output/util_pythran.py +++ /dev/null @@ -1,23 +0,0 @@ -# pythran export strfunc_from_pdf( -# float64[][], float64[][], float64[][], float, bool) - -# pythran export strfunc_from_pdf( -# int32[], float64[][], float64[][], float, bool) - - -import numpy as np - - -def strfunc_from_pdf(rxs, pdf, values, order, absolute=False): - """Compute structure function of specified order from pdf for increments - module. - - """ - S_order = np.empty(rxs.shape) - if absolute: - values = abs(values) - for irx in range(rxs.size): - deltainc = abs(values[irx, 1] - values[irx, 0]) - S_order[irx] = deltainc * np.sum(pdf[irx] * values[irx] ** order) - - return S_order diff --git a/fluidsim/solvers/ns2d/bouss/solver.py b/fluidsim/solvers/ns2d/bouss/solver.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9uczJkL2JvdXNzL3NvbHZlci5weQ==..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vc29sdmVycy9uczJkL2JvdXNzL3NvbHZlci5weQ== 100644 --- a/fluidsim/solvers/ns2d/bouss/solver.py +++ b/fluidsim/solvers/ns2d/bouss/solver.py @@ -6,8 +6,5 @@ :private-members: """ - -from __future__ import division - import numpy as np @@ -12,3 +9,5 @@ import numpy as np +from fluidpythran import cachedjit, Array + from fluidsim.base.setofvariables import SetOfVariables @@ -14,4 +13,3 @@ from fluidsim.base.setofvariables import SetOfVariables - from fluidsim.solvers.ns2d.solver import InfoSolverNS2D, Simul as SimulNS2D @@ -16,6 +14,16 @@ from fluidsim.solvers.ns2d.solver import InfoSolverNS2D, Simul as SimulNS2D -from .util_pythran import tendencies_nonlin_ns2dbouss + +AF = Array[np.float64, "2d"] + + +@cachedjit +def tendencies_nonlin_ns2dbouss( + ux: AF, uy: AF, px_rot: AF, py_rot: AF, px_b: AF, py_b: AF +): + Frot = -ux * px_rot - uy * py_rot + px_b + Fb = -ux * px_b - uy * py_b + return Frot, Fb class InfoSolverNS2DBouss(InfoSolverNS2D): diff --git a/fluidsim/solvers/ns2d/bouss/util_pythran.py b/fluidsim/solvers/ns2d/bouss/util_pythran.py deleted file mode 100644 index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9uczJkL2JvdXNzL3V0aWxfcHl0aHJhbi5weQ==..0000000000000000000000000000000000000000 --- a/fluidsim/solvers/ns2d/bouss/util_pythran.py +++ /dev/null @@ -1,11 +0,0 @@ -# pythran export tendencies_nonlin_ns2dbouss( -# float64[][], float64[][], float64[][], float64[][], -# float64[][], float64[][]) - - -def tendencies_nonlin_ns2dbouss(ux, uy, px_rot, py_rot, px_b, py_b): - - Frot = -ux * px_rot - uy * py_rot + px_b - Fb = -ux * px_b - uy * py_b - - return Frot, Fb diff --git a/fluidsim/solvers/ns2d/forcing.py b/fluidsim/solvers/ns2d/forcing.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9uczJkL2ZvcmNpbmcucHk=..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vc29sdmVycy9uczJkL2ZvcmNpbmcucHk= 100644 --- a/fluidsim/solvers/ns2d/forcing.py +++ b/fluidsim/solvers/ns2d/forcing.py @@ -10,7 +10,7 @@ from fluidsim.base.forcing import ForcingBasePseudoSpectral from fluidsim.base.forcing.anisotropic import ( - TimeCorrelatedRandomPseudoSpectralAnisotropic, + TimeCorrelatedRandomPseudoSpectralAnisotropic ) diff --git a/fluidsim/solvers/ns2d/strat/solver.py b/fluidsim/solvers/ns2d/strat/solver.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3NvbHZlci5weQ==..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3NvbHZlci5weQ== 100644 --- a/fluidsim/solvers/ns2d/strat/solver.py +++ b/fluidsim/solvers/ns2d/strat/solver.py @@ -6,7 +6,6 @@ :private-members: """ -from __future__ import division import numpy as np @@ -10,7 +9,9 @@ import numpy as np +from fluidpythran import cachedjit, Array + from fluidsim.base.setofvariables import SetOfVariables from fluidsim.solvers.ns2d.solver import InfoSolverNS2D, Simul as SimulNS2D @@ -13,8 +14,17 @@ from fluidsim.base.setofvariables import SetOfVariables from fluidsim.solvers.ns2d.solver import InfoSolverNS2D, Simul as SimulNS2D -from .util_pythran import tendencies_nonlin_ns2dstrat +AF = Array[np.float64, "2d"] + + +@cachedjit +def tendencies_nonlin_ns2dstrat( + ux: AF, uy: AF, px_rot: AF, py_rot: AF, px_b: AF, py_b: AF, N: float +): + Frot = -ux * px_rot - uy * py_rot + Fb = -ux * px_b - uy * py_b - N ** 2 * uy + return Frot, Fb class InfoSolverNS2DStrat(InfoSolverNS2D): diff --git a/fluidsim/solvers/ns2d/strat/util_pythran.py b/fluidsim/solvers/ns2d/strat/util_pythran.py deleted file mode 100644 index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9uczJkL3N0cmF0L3V0aWxfcHl0aHJhbi5weQ==..0000000000000000000000000000000000000000 --- a/fluidsim/solvers/ns2d/strat/util_pythran.py +++ /dev/null @@ -1,11 +0,0 @@ -# pythran export tendencies_nonlin_ns2dstrat( -# float64[][], float64[][], float64[][], float64[][], -# float64[][], float64[][], float) - - -def tendencies_nonlin_ns2dstrat(ux, uy, px_rot, py_rot, px_b, py_b, N): - - Frot = -ux * px_rot - uy * py_rot - Fb = -ux * px_b - uy * py_b - N ** 2 * uy - - return Frot, Fb diff --git a/fluidsim/solvers/plate2d/operators.py b/fluidsim/solvers/plate2d/operators.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL29wZXJhdG9ycy5weQ==..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL29wZXJhdG9ycy5weQ== 100644 --- a/fluidsim/solvers/plate2d/operators.py +++ b/fluidsim/solvers/plate2d/operators.py @@ -9,4 +9,5 @@ """ +import numpy as np @@ -12,5 +13,5 @@ -import numpy as np +from fluidpythran import cachedjit, Array from ...operators.operators2d import OperatorsPseudoSpectral2D @@ -14,7 +15,26 @@ from ...operators.operators2d import OperatorsPseudoSpectral2D -from .util_oper_pythran import monge_ampere_step0, monge_ampere_step1 +AC2 = Array[np.complex128, "2d"] +AF2 = Array[np.float64, "2d"] + + +@cachedjit +def monge_ampere_step0(a_fft: AC2, b_fft: AC2, KX2: AF2, KY2: AF2, KXKZ: AF2): + pxx_a_fft = -a_fft * KX2 + pyy_a_fft = -a_fft * KY2 + pxy_a_fft = -a_fft * KXKZ + pxx_b_fft = -b_fft * KX2 + pyy_b_fft = -b_fft * KY2 + pxy_b_fft = -b_fft * KXKZ + return pxx_a_fft, pyy_a_fft, pxy_a_fft, pxx_b_fft, pyy_b_fft, pxy_b_fft + + +@cachedjit +def monge_ampere_step1( + pxx_a: AF2, pyy_a: AF2, pxy_a: AF2, pxx_b: AF2, pyy_b: AF2, pxy_b: AF2 +): + return pxx_a * pyy_b + pyy_a * pxx_b - 2 * pxy_a * pxy_b class OperatorsPseudoSpectralPlate2D(OperatorsPseudoSpectral2D): @@ -37,10 +57,6 @@ self.tmp_pyy_b = np.empty_like(self.XX) self.tmp_pxy_b = np.empty_like(self.XX) - # def monge_ampere_from_fft2(self, a_fft, b_fft): - # return monge_ampere_from_fft( - # a_fft, b_fft, self.KX, self.KY, self.ifft2) - def monge_ampere_from_fft(self, a_fft, b_fft): """Compute the Monge-Ampere operator""" @@ -65,7 +81,7 @@ ) -def monge_ampere_from_fft2(a_fft, b_fft, KX, KY, ifft2): +def monge_ampere_from_fft_numpy(a_fft, b_fft, KX, KY, ifft2): pxx_a = -ifft2(a_fft * KX ** 2) pyy_a = -ifft2(a_fft * KY ** 2) diff --git a/fluidsim/solvers/plate2d/output/correlations_freq.py b/fluidsim/solvers/plate2d/output/correlations_freq.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9jb3JyZWxhdGlvbnNfZnJlcS5weQ==..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC9jb3JyZWxhdGlvbnNfZnJlcS5weQ== 100644 --- a/fluidsim/solvers/plate2d/output/correlations_freq.py +++ b/fluidsim/solvers/plate2d/output/correlations_freq.py @@ -17,6 +17,8 @@ from builtins import range import os import numpy as np + +# pythran import numpy as np import h5py import matplotlib.pyplot as plt @@ -20,8 +22,9 @@ import h5py import matplotlib.pyplot as plt +from fluidpythran import cachedjit from fluiddyn.util import mpi from fluidsim.base.output.base import SpecificOutput from fluiddyn.calcul.easypyfft import FFTW1DReal2Complex @@ -23,7 +26,38 @@ from fluiddyn.util import mpi from fluidsim.base.output.base import SpecificOutput from fluiddyn.calcul.easypyfft import FFTW1DReal2Complex -# from fluidsim.operators.miscellaneous import compute_correl4, compute_correl2 + +@cachedjit +def compute_correl4_seq( + q_fftt: "complex128[][]", iomegas1: "int32[]", nb_omegas: int, nb_xs_seq: int +): + r"""Compute the correlations 4. + + .. math:: + C_4(\omega_1, \omega_2, \omega_3, \omega_4) = + \langle + \tilde w(\omega_1, \mathbf{x}) + \tilde w(\omega_2, \mathbf{x}) + \tilde w(\omega_3, \mathbf{x})^* + \tilde w(\omega_4, \mathbf{x})^* + \rangle_\mathbf{x}, + + where + + .. math:: + \omega_2 = \omega_3 + \omega_4 - \omega_1 + + and :math:`\omega_1 > 0`, :math:`\omega_3 > 0` and + :math:`\omega_4 > 0`. Thus, this function produces an array + :math:`C_4(\omega_1, \omega_3, \omega_4)`. + + """ + q_fftt_conj = np.conj(q_fftt) + + nx = q_fftt.shape[0] + n0 = iomegas1.shape[0] + + corr4 = np.zeros((len(iomegas1), nb_omegas, nb_omegas), dtype=np.complex128) @@ -29,5 +63,71 @@ -from .util_pythran import compute_correl2_seq, compute_correl4_seq + for i1 in range(n0): + io1 = iomegas1[i1] + # this loop could be parallelized (OMP) + for io3 in range(nb_omegas): + # we use the symmetry omega_3 <--> omega_4 + for io4 in range(io3 + 1): + io2 = io3 + io4 - io1 + if io2 < 0: + io2 = -io2 + for ix in range(nx): + corr4[i1, io3, io4] += ( + q_fftt[ix, io1] + * q_fftt_conj[ix, io3] + * q_fftt_conj[ix, io4] + * q_fftt_conj[ix, io2] + ) + elif io2 >= nb_omegas: + io2 = 2 * nb_omegas - 1 - io2 + for ix in range(nx): + corr4[i1, io3, io4] += ( + q_fftt[ix, io1] + * q_fftt_conj[ix, io3] + * q_fftt_conj[ix, io4] + * q_fftt_conj[ix, io2] + ) + else: + for ix in range(nx): + corr4[i1, io3, io4] += ( + q_fftt[ix, io1] + * q_fftt_conj[ix, io3] + * q_fftt_conj[ix, io4] + * q_fftt[ix, io2] + ) + # symmetry omega_3 <--> omega_4: + corr4[i1, io4, io3] = corr4[i1, io3, io4] + return corr4 + + +@cachedjit +def compute_correl2_seq( + q_fftt: "complex128[][]", iomegas1: "int32[]", nb_omegas: int, nb_xs_seq: int +): + r"""Compute the correlations 2. + + .. math:: + C_2(\omega_1, \omega_2) = + \langle + \tilde w(\omega_1, \mathbf{x}) + \tilde w(\omega_2, \mathbf{x})^* + \rangle_\mathbf{x}, + + where :math:`\omega_1 = \omega_2`. Thus, this function + produces an array :math:`C_2(\omega)`. + + """ + q_fftt_conj = np.conj(q_fftt) + nx = q_fftt.shape[0] + + corr2 = np.zeros((nb_omegas, nb_omegas), dtype=np.complex128) + + for io3 in range(nb_omegas): + for io4 in range(io3 + 1): + for ix in range(nx): + corr2[io3, io4] += q_fftt[ix, io3] * q_fftt_conj[ix, io4] + corr2[io4, io3] = np.conj(corr2[io3, io4]) + + return corr2 def compute_correl4(q_fftt, iomegas1, nb_omegas, nb_xs_seq): diff --git a/fluidsim/solvers/plate2d/output/util_pythran.py b/fluidsim/solvers/plate2d/output/util_pythran.py deleted file mode 100644 index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL291dHB1dC91dGlsX3B5dGhyYW4ucHk=..0000000000000000000000000000000000000000 --- a/fluidsim/solvers/plate2d/output/util_pythran.py +++ /dev/null @@ -1,102 +0,0 @@ -import numpy as np - - -# pythran export compute_correl4_seq(complex128[][], int32[], int, int) - - -def compute_correl4_seq(q_fftt, iomegas1, nb_omegas, nb_xs_seq): - r"""Compute the correlations 4. - - .. math:: - C_4(\omega_1, \omega_2, \omega_3, \omega_4) = - \langle - \tilde w(\omega_1, \mathbf{x}) - \tilde w(\omega_2, \mathbf{x}) - \tilde w(\omega_3, \mathbf{x})^* - \tilde w(\omega_4, \mathbf{x})^* - \rangle_\mathbf{x}, - - where - - .. math:: - \omega_2 = \omega_3 + \omega_4 - \omega_1 - - and :math:`\omega_1 > 0`, :math:`\omega_3 > 0` and - :math:`\omega_4 > 0`. Thus, this function produces an array - :math:`C_4(\omega_1, \omega_3, \omega_4)`. - - """ - q_fftt_conj = np.conj(q_fftt) - - nx = q_fftt.shape[0] - n0 = iomegas1.shape[0] - - corr4 = np.zeros((len(iomegas1), nb_omegas, nb_omegas), dtype=np.complex128) - - for i1 in range(n0): - io1 = iomegas1[i1] - # this loop could be parallelized (OMP) - for io3 in range(nb_omegas): - # we use the symmetry omega_3 <--> omega_4 - for io4 in range(io3 + 1): - io2 = io3 + io4 - io1 - if io2 < 0: - io2 = -io2 - for ix in range(nx): - corr4[i1, io3, io4] += ( - q_fftt[ix, io1] - * q_fftt_conj[ix, io3] - * q_fftt_conj[ix, io4] - * q_fftt_conj[ix, io2] - ) - elif io2 >= nb_omegas: - io2 = 2 * nb_omegas - 1 - io2 - for ix in range(nx): - corr4[i1, io3, io4] += ( - q_fftt[ix, io1] - * q_fftt_conj[ix, io3] - * q_fftt_conj[ix, io4] - * q_fftt_conj[ix, io2] - ) - else: - for ix in range(nx): - corr4[i1, io3, io4] += ( - q_fftt[ix, io1] - * q_fftt_conj[ix, io3] - * q_fftt_conj[ix, io4] - * q_fftt[ix, io2] - ) - # symmetry omega_3 <--> omega_4: - corr4[i1, io4, io3] = corr4[i1, io3, io4] - return corr4 - - -# pythran export compute_correl2_seq(complex128[][], int32[], int, int) - - -def compute_correl2_seq(q_fftt, iomegas1, nb_omegas, nb_xs_seq): - r"""Compute the correlations 2. - - .. math:: - C_2(\omega_1, \omega_2) = - \langle - \tilde w(\omega_1, \mathbf{x}) - \tilde w(\omega_2, \mathbf{x})^* - \rangle_\mathbf{x}, - - where :math:`\omega_1 = \omega_2`. Thus, this function - produces an array :math:`C_2(\omega)`. - - """ - q_fftt_conj = np.conj(q_fftt) - nx = q_fftt.shape[0] - - corr2 = np.zeros((nb_omegas, nb_omegas), dtype=np.complex128) - - for io3 in range(nb_omegas): - for io4 in range(io3 + 1): - for ix in range(nx): - corr2[io3, io4] += q_fftt[ix, io3] * q_fftt_conj[ix, io4] - corr2[io4, io3] = np.conj(corr2[io3, io4]) - - return corr2 diff --git a/fluidsim/solvers/plate2d/util_oper_pythran.py b/fluidsim/solvers/plate2d/util_oper_pythran.py deleted file mode 100644 index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9wbGF0ZTJkL3V0aWxfb3Blcl9weXRocmFuLnB5..0000000000000000000000000000000000000000 --- a/fluidsim/solvers/plate2d/util_oper_pythran.py +++ /dev/null @@ -1,22 +0,0 @@ -# pythran export monge_ampere_step0( -# complex128[][], complex128[][], -# float64[][], float64[][], float64[][]) - - -def monge_ampere_step0(a_fft, b_fft, KX2, KY2, KXKZ): - pxx_a_fft = -a_fft * KX2 - pyy_a_fft = -a_fft * KY2 - pxy_a_fft = -a_fft * KXKZ - pxx_b_fft = -b_fft * KX2 - pyy_b_fft = -b_fft * KY2 - pxy_b_fft = -b_fft * KXKZ - return pxx_a_fft, pyy_a_fft, pxy_a_fft, pxx_b_fft, pyy_b_fft, pxy_b_fft - - -# pythran export monge_ampere_step1( -# float64[][], float64[][], float64[][], -# float64[][], float64[][], float64[][]) - - -def monge_ampere_step1(pxx_a, pyy_a, pxy_a, pxx_b, pyy_b, pxy_b): - return pxx_a * pyy_b + pyy_a * pxx_b - 2 * pxy_a * pxy_b diff --git a/fluidsim/solvers/sw1l/init_fields.py b/fluidsim/solvers/sw1l/init_fields.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2luaXRfZmllbGRzLnB5..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vc29sdmVycy9zdzFsL2luaXRfZmllbGRzLnB5 100644 --- a/fluidsim/solvers/sw1l/init_fields.py +++ b/fluidsim/solvers/sw1l/init_fields.py @@ -33,7 +33,7 @@ from fluidsim.base.init_fields import InitFieldsBase, SpecificInitFields from fluidsim.solvers.ns2d.init_fields import ( - InitFieldsNoise as InitFieldsNoiseNS2D, + InitFieldsNoise as InitFieldsNoiseNS2D ) from fluidsim.solvers.ns2d.init_fields import InitFieldsJet, InitFieldsDipole diff --git a/fluidsim/solvers/sw1l/operators.py b/fluidsim/solvers/sw1l/operators.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9zdzFsL29wZXJhdG9ycy5weQ==..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vc29sdmVycy9zdzFsL29wZXJhdG9ycy5weQ== 100644 --- a/fluidsim/solvers/sw1l/operators.py +++ b/fluidsim/solvers/sw1l/operators.py @@ -9,5 +9,11 @@ """ +import numpy as np + +# pythran import numpy as np + +from fluidpythran import cachedjit, Array + from fluidsim.operators.operators2d import OperatorsPseudoSpectral2D, rank @@ -12,6 +18,90 @@ from fluidsim.operators.operators2d import OperatorsPseudoSpectral2D, rank -from . import util_oper_pythran +AC = Array[np.complex128, "2d"] +AF = Array[np.float64, "2d"] + + +@cachedjit +def _qapamfft_from_uxuyetafft( + ux_fft: AC, + uy_fft: AC, + eta_fft: AC, + n0: int, + n1: int, + KX: AF, + KY: AF, + K2: AF, + Kappa_over_ic: AC, + f: float, + c2: float, + rank: int, +): + """Calculate normal modes from primitive variables.""" + freq_Corio = f + f_over_c2 = freq_Corio / c2 + + q_fft = np.empty([n0, n1], dtype=np.complex128) + ap_fft = np.empty([n0, n1], dtype=np.complex128) + am_fft = np.empty([n0, n1], dtype=np.complex128) + + if freq_Corio != 0: + for i0 in range(n0): + for i1 in range(n1): + if i0 == 0 and i1 == 0 and rank == 0: + q_fft[i0, i1] = 0 + ap_fft[i0, i1] = ux_fft[0, 0] + 1.0j * uy_fft[0, 0] + am_fft[i0, i1] = ux_fft[0, 0] - 1.0j * uy_fft[0, 0] + else: + + rot_fft = 1j * ( + KX[i0, i1] * uy_fft[i0, i1] - KY[i0, i1] * ux_fft[i0, i1] + ) + + q_fft[i0, i1] = rot_fft - freq_Corio * eta_fft[i0, i1] + + a_over2_fft = 0.5 * ( + K2[i0, i1] * eta_fft[i0, i1] + f_over_c2 * rot_fft + ) + + Deltaa_over2_fft = ( + 0.5j + * Kappa_over_ic[i0, i1] + * ( + KX[i0, i1] * ux_fft[i0, i1] + + KY[i0, i1] * uy_fft[i0, i1] + ) + ) + + ap_fft[i0, i1] = a_over2_fft + Deltaa_over2_fft + am_fft[i0, i1] = a_over2_fft - Deltaa_over2_fft + + else: # (freq_Corio == 0.) + for i0 in range(n0): + for i1 in range(n1): + if i0 == 0 and i1 == 0 and rank == 0: + q_fft[i0, i1] = 0 + ap_fft[i0, i1] = ux_fft[0, 0] + 1.0j * uy_fft[0, 0] + am_fft[i0, i1] = ux_fft[0, 0] - 1.0j * uy_fft[0, 0] + else: + q_fft[i0, i1] = 1j * ( + KX[i0, i1] * uy_fft[i0, i1] - KY[i0, i1] * ux_fft[i0, i1] + ) + + a_over2_fft = 0.5 * K2[i0, i1] * eta_fft[i0, i1] + + Deltaa_over2_fft = ( + 0.5j + * Kappa_over_ic[i0, i1] + * ( + KX[i0, i1] * ux_fft[i0, i1] + + KY[i0, i1] * uy_fft[i0, i1] + ) + ) + + ap_fft[i0, i1] = a_over2_fft + Deltaa_over2_fft + am_fft[i0, i1] = a_over2_fft - Deltaa_over2_fft + + return q_fft, ap_fft, am_fft class OperatorsPseudoSpectralSW1L(OperatorsPseudoSpectral2D): @@ -31,7 +121,7 @@ f = float(params.f) c2 = float(params.c2) - return util_oper_pythran.qapamfft_from_uxuyetafft( + return _qapamfft_from_uxuyetafft( ux_fft, uy_fft, eta_fft, diff --git a/fluidsim/solvers/sw1l/output/__init__.py b/fluidsim/solvers/sw1l/output/__init__.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9fX2luaXRfXy5weQ==..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9fX2luaXRfXy5weQ== 100644 --- a/fluidsim/solvers/sw1l/output/__init__.py +++ b/fluidsim/solvers/sw1l/output/__init__.py @@ -15,7 +15,6 @@ normal_mode """ -from __future__ import division import numpy as np @@ -19,4 +18,6 @@ import numpy as np +from fluidpythran import cachedjit + from fluidsim.base.output import OutputBasePseudoSpectral @@ -22,5 +23,23 @@ from fluidsim.base.output import OutputBasePseudoSpectral -from .util_pythran import linear_eigenmode_from_values_1k + + +@cachedjit +def linear_eigenmode_from_values_1k( + ux_fft: np.complex128, + uy_fft: np.complex128, + eta_fft: np.complex128, + kx: float, + ky: float, + f: "float or int", + c2: "float or int", +): + """Compute q, d, a (fft) for a single wavenumber.""" + div_fft = 1j * (kx * ux_fft + ky * uy_fft) + rot_fft = 1j * (kx * uy_fft - ky * ux_fft) + q_fft = rot_fft - f * eta_fft + k2 = kx ** 2 + ky ** 2 + ageo_fft = f * rot_fft / c2 + k2 * eta_fft + return q_fft, div_fft, ageo_fft class OutputBaseSW1L(OutputBasePseudoSpectral): diff --git a/fluidsim/solvers/sw1l/output/normal_mode.py b/fluidsim/solvers/sw1l/output/normal_mode.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9ub3JtYWxfbW9kZS5weQ==..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC9ub3JtYWxfbW9kZS5weQ== 100644 --- a/fluidsim/solvers/sw1l/output/normal_mode.py +++ b/fluidsim/solvers/sw1l/output/normal_mode.py @@ -20,4 +20,8 @@ import numpy as np +# pythran import numpy as np + +from fluidpythran import cachedjit, Array + from fluiddyn.util import mpi @@ -23,5 +27,41 @@ from fluiddyn.util import mpi -from .util_pythran import get_qmat + + +AF = Array[float, "2d"] + + +@cachedjit +def get_qmat( + f: "float or int", + c: "float or int", + sigma: AF, + KX: AF, + KY: AF, + K: AF, + K2: AF, + K_not0: AF, +): + """Compute Q matrix to transform q, ap, am (fft) -> b0, b+, b- (fft) with + dimensions of velocity. + + """ + ck = c * K_not0 + qmat = np.array( + [ + [ + -1j * 2.0 ** 0.5 * ck * KY, + +1j * f * KY + KX * sigma, + +1j * f * KY - KX * sigma, + ], + [ + +1j * 2.0 ** 0.5 * ck * KX, + -1j * f * KX + KY * sigma, + -1j * f * KX - KY * sigma, + ], + [2.0 ** 0.5 * f * K, c * K2, c * K2], + ] + ) / (2.0 ** 0.5 * sigma * K_not0) + return qmat class NormalModeBase(object): @@ -79,7 +119,7 @@ return bvec_fft def bvecfft_from_uxuyetafft(self, ux_fft, uy_fft, eta_fft): - r""" Compute normal mode vector, :math:`\mathbf{B}` + r""" Compute normal mode vector, :math:`\mathbf{B}` with dimensions of velocity from primitive variables. """ q_fft, ap_fft, am_fft = self.oper.qapamfft_from_uxuyetafft( ux_fft, uy_fft, eta_fft diff --git a/fluidsim/solvers/sw1l/output/util_pythran.py b/fluidsim/solvers/sw1l/output/util_pythran.py deleted file mode 100644 index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9zdzFsL291dHB1dC91dGlsX3B5dGhyYW4ucHk=..0000000000000000000000000000000000000000 --- a/fluidsim/solvers/sw1l/output/util_pythran.py +++ /dev/null @@ -1,45 +0,0 @@ -import numpy as np - - -# pythran export get_qmat( -# float or int, float or int, float64[][], float64[][], float64[][], -# float64[][], float64[][], float64[][]) - - -def get_qmat(f, c, sigma, KX, KY, K, K2, K_not0): - """Compute Q matrix to transform q, ap, am (fft) -> b0, b+, b- (fft) with - dimensions of velocity. - - """ - ck = c * K_not0 - qmat = np.array( - [ - [ - -1j * 2.0 ** 0.5 * ck * KY, - +1j * f * KY + KX * sigma, - +1j * f * KY - KX * sigma, - ], - [ - +1j * 2.0 ** 0.5 * ck * KX, - -1j * f * KX + KY * sigma, - -1j * f * KX - KY * sigma, - ], - [2.0 ** 0.5 * f * K, c * K2, c * K2], - ] - ) / (2.0 ** 0.5 * sigma * K_not0) - return qmat - - -# pythran export linear_eigenmode_from_values_1k( -# complex128, complex128, complex128, float, float, -# float or int, float or int) - - -def linear_eigenmode_from_values_1k(ux_fft, uy_fft, eta_fft, kx, ky, f, c2): - """Compute q, d, a (fft) for a single wavenumber.""" - div_fft = 1j * (kx * ux_fft + ky * uy_fft) - rot_fft = 1j * (kx * uy_fft - ky * ux_fft) - q_fft = rot_fft - f * eta_fft - k2 = kx ** 2 + ky ** 2 - ageo_fft = f * rot_fft / c2 + k2 * eta_fft - return q_fft, div_fft, ageo_fft diff --git a/fluidsim/solvers/sw1l/solver.py b/fluidsim/solvers/sw1l/solver.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3NvbHZlci5weQ==..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3NvbHZlci5weQ== 100644 --- a/fluidsim/solvers/sw1l/solver.py +++ b/fluidsim/solvers/sw1l/solver.py @@ -16,7 +16,9 @@ """ -from __future__ import division + +from fluidpythran import cachedjit, Array +from fluiddyn.util import mpi from fluidsim.base.setofvariables import SetOfVariables from fluidsim.base.solvers.pseudo_spect import ( @@ -24,5 +26,6 @@ InfoSolverPseudoSpectral, ) -from fluiddyn.util import mpi +A = Array[float, "2d"] + @@ -28,5 +31,16 @@ -from .util_pythran import compute_Frot +@cachedjit +def compute_Frot(rot: A, ux: A, uy: A, f: "float or int"): + """Compute cross-product of absolute potential vorticity with velocity.""" + if f != 0: + rot_abs = rot + f + else: + rot_abs = rot + + F1x = rot_abs * uy + F1y = -rot_abs * ux + + return F1x, F1y class InfoSolverSW1L(InfoSolverPseudoSpectral): diff --git a/fluidsim/solvers/sw1l/util_oper_pythran.py b/fluidsim/solvers/sw1l/util_oper_pythran.py deleted file mode 100644 index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3V0aWxfb3Blcl9weXRocmFuLnB5..0000000000000000000000000000000000000000 --- a/fluidsim/solvers/sw1l/util_oper_pythran.py +++ /dev/null @@ -1,78 +0,0 @@ -import numpy as np - - -# pythran export qapamfft_from_uxuyetafft( -# complex128[][], complex128[][], complex128[][], -# int, int, float64[][], float64[][], float64[][], -# complex128[][], float, float, int) - - -def qapamfft_from_uxuyetafft( - ux_fft, uy_fft, eta_fft, n0, n1, KX, KY, K2, Kappa_over_ic, f, c2, rank -): - """Calculate normal modes from primitive variables.""" - freq_Corio = f - f_over_c2 = freq_Corio / c2 - - q_fft = np.empty([n0, n1], dtype=np.complex128) - ap_fft = np.empty([n0, n1], dtype=np.complex128) - am_fft = np.empty([n0, n1], dtype=np.complex128) - - if freq_Corio != 0: - for i0 in range(n0): - for i1 in range(n1): - if i0 == 0 and i1 == 0 and rank == 0: - q_fft[i0, i1] = 0 - ap_fft[i0, i1] = ux_fft[0, 0] + 1.0j * uy_fft[0, 0] - am_fft[i0, i1] = ux_fft[0, 0] - 1.0j * uy_fft[0, 0] - else: - - rot_fft = 1j * ( - KX[i0, i1] * uy_fft[i0, i1] - KY[i0, i1] * ux_fft[i0, i1] - ) - - q_fft[i0, i1] = rot_fft - freq_Corio * eta_fft[i0, i1] - - a_over2_fft = 0.5 * ( - K2[i0, i1] * eta_fft[i0, i1] + f_over_c2 * rot_fft - ) - - Deltaa_over2_fft = ( - 0.5j - * Kappa_over_ic[i0, i1] - * ( - KX[i0, i1] * ux_fft[i0, i1] - + KY[i0, i1] * uy_fft[i0, i1] - ) - ) - - ap_fft[i0, i1] = a_over2_fft + Deltaa_over2_fft - am_fft[i0, i1] = a_over2_fft - Deltaa_over2_fft - - else: # (freq_Corio == 0.) - for i0 in range(n0): - for i1 in range(n1): - if i0 == 0 and i1 == 0 and rank == 0: - q_fft[i0, i1] = 0 - ap_fft[i0, i1] = ux_fft[0, 0] + 1.0j * uy_fft[0, 0] - am_fft[i0, i1] = ux_fft[0, 0] - 1.0j * uy_fft[0, 0] - else: - q_fft[i0, i1] = 1j * ( - KX[i0, i1] * uy_fft[i0, i1] - KY[i0, i1] * ux_fft[i0, i1] - ) - - a_over2_fft = 0.5 * K2[i0, i1] * eta_fft[i0, i1] - - Deltaa_over2_fft = ( - 0.5j - * Kappa_over_ic[i0, i1] - * ( - KX[i0, i1] * ux_fft[i0, i1] - + KY[i0, i1] * uy_fft[i0, i1] - ) - ) - - ap_fft[i0, i1] = a_over2_fft + Deltaa_over2_fft - am_fft[i0, i1] = a_over2_fft - Deltaa_over2_fft - - return q_fft, ap_fft, am_fft diff --git a/fluidsim/solvers/sw1l/util_pythran.py b/fluidsim/solvers/sw1l/util_pythran.py deleted file mode 100644 index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vc29sdmVycy9zdzFsL3V0aWxfcHl0aHJhbi5weQ==..0000000000000000000000000000000000000000 --- a/fluidsim/solvers/sw1l/util_pythran.py +++ /dev/null @@ -1,16 +0,0 @@ -# pythran export compute_Frot( -# float64[][], float64[][], float64[][], -# float or int) - - -def compute_Frot(rot, ux, uy, f): - """Compute cross-product of absolute potential vorticity with velocity.""" - if f != 0: - rot_abs = rot + f - else: - rot_abs = rot - - F1x = rot_abs * uy - F1y = -rot_abs * ux - - return F1x, F1y diff --git a/fluidsim/util/util.py b/fluidsim/util/util.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_Zmx1aWRzaW0vdXRpbC91dGlsLnB5..a02e313946d37e785ae750c7e61647ed861c80d0_Zmx1aWRzaW0vdXRpbC91dGlsLnB5 100644 --- a/fluidsim/util/util.py +++ b/fluidsim/util/util.py @@ -196,7 +196,7 @@ for ii, name in enumerate(name_files): tmp = ".".join(name[ind_start_time:].split(".")[:2]) if "_" in tmp: - tmp = tmp[:tmp.index("_")] + tmp = tmp[: tmp.index("_")] times[ii] = float(tmp) if t_approx is None: t_approx = times.max() diff --git a/pyproject.toml b/pyproject.toml index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_cHlwcm9qZWN0LnRvbWw=..a02e313946d37e785ae750c7e61647ed861c80d0_cHlwcm9qZWN0LnRvbWw= 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,2 +1,2 @@ [build-system] -requires = ["setuptools", "wheel", "numpy", "fluidpythran>=0.1.0.post0"] \ No newline at end of file +requires = ["setuptools", "wheel", "numpy", "fluidpythran>=0.1.2"] \ No newline at end of file diff --git a/setup.cfg b/setup.cfg index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_c2V0dXAuY2Zn..a02e313946d37e785ae750c7e61647ed861c80d0_c2V0dXAuY2Zn 100644 --- a/setup.cfg +++ b/setup.cfg @@ -21,5 +21,5 @@ [options] setup_requires = - fluidpythran>=0.1.1 + fluidpythran>=0.1.2 numpy diff --git a/setup.py b/setup.py index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_c2V0dXAucHk=..a02e313946d37e785ae750c7e61647ed861c80d0_c2V0dXAucHk= 100644 --- a/setup.py +++ b/setup.py @@ -53,8 +53,11 @@ from fluidpythran.dist import make_pythran_files -paths = ["fluidsim/base/time_stepping/pseudo_spect.py"] +paths = [ + "fluidsim/base/time_stepping/pseudo_spect.py", + "fluidsim/base/output/increments.py", +] make_pythran_files( [here / path for path in paths], mocked_modules=( "psutil", @@ -57,8 +60,9 @@ make_pythran_files( [here / path for path in paths], mocked_modules=( "psutil", + "h5py", "fluiddyn", "fluiddyn.io", "fluiddyn.util", "fluiddyn.util.paramcontainer", @@ -61,8 +65,8 @@ "fluiddyn", "fluiddyn.io", "fluiddyn.util", "fluiddyn.util.paramcontainer", - "h5py", + "fluiddyn.output", ), ) @@ -133,7 +137,7 @@ "future >= 0.16", "h5py", "h5netcdf", - "fluidpythran >= 0.1.1", + "fluidpythran>=0.1.2", "setuptools_scm", "xarray", ] @@ -148,9 +152,7 @@ def make_pythran_extensions(modules): - exclude_pythran = tuple( - key for key, value in config["exclude_pythran"].items() if value - ) + exclude_pythran = tuple() if len(exclude_pythran) > 0: logger.info( "Pythran files in the packages " diff --git a/site.cfg.default b/site.cfg.default index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_c2l0ZS5jZmcuZGVmYXVsdA==..a02e313946d37e785ae750c7e61647ed861c80d0_c2l0ZS5jZmcuZGVmYXVsdA== 100644 --- a/site.cfg.default +++ b/site.cfg.default @@ -1,11 +1,3 @@ -[exclude_pythran] -fluidsim.solvers.plate2d = True -fluidsim.solvers.plate2d.output = True -fluidsim.solvers.sw1l = True -fluidsim.solvers.sw1l.output = True -fluidsim.solvers.ns2d.strat = True -fluidsim.solvers.ns2d.bouss = True - [environ] diff --git a/tox.ini b/tox.ini index 208ccc560904f1b8d4e321003ea2fa36a1721f7e_dG94LmluaQ==..a02e313946d37e785ae750c7e61647ed861c80d0_dG94LmluaQ== 100644 --- a/tox.ini +++ b/tox.ini @@ -22,7 +22,7 @@ deps = pip==18.0 coverage - fluidpythran>=0.1.1 + fluidpythran>=0.1.2 mako numpy matplotlib