Skip to content
Snippets Groups Projects
Commit 1af1fb70 authored by Pierre Augier's avatar Pierre Augier
Browse files

Towards ns3d solver.

parent a8ec70e9
No related branches found
No related tags found
No related merge requests found
......@@ -52,6 +52,17 @@
'class_name': 'OperatorsPseudoSpectral2D'})
class InfoSolverPseudoSpectral3D(InfoSolverPseudoSpectral):
"""Contain the information on a solver."""
def _init_root(self):
super(InfoSolverPseudoSpectral3D, self)._init_root()
self.Operators.module_name = 'fluidsim.operators.operators3d'
self.Operators.class_name = 'OperatorsPseudoSpectral3D'
class SimulBasePseudoSpectral(SimulBase):
InfoSolver = InfoSolverPseudoSpectral
......
import numpy as np
from math import pi
from fluidsim.operators.fft.easypyfft import FFTW3DReal2Complex
class OperatorsPseudoSpectral3D(object):
"""Provides fast Fourier transform functions and 3D operators.
"""
@staticmethod
def _complete_params_with_default(params):
"""This static method is used to complete the *params* container.
"""
attribs = {'type_fft': type_fft,
'TRANSPOSED_OK': True,
'coef_dealiasing': 2./3,
'nx': 48,
'ny': 48,
'nz': 48,
'Lx': 2*pi,
'Ly': 2*pi,
'Lz': 2*pi}
params._set_child('oper', attribs=attribs)
def __init__(self, params=None, SEQUENTIAL=None):
self.params = params
nx = self.nx = params.oper.nx
ny = self.ny = params.oper.ny
nz = self.nz = params.oper.nz
self.shape_phys = (nz, ny, nx)
Lx = self.Lx = params.oper.Lx
Ly = self.Ly = params.oper.Ly
Lz = self.Lz = params.oper.Lz
if nx % 2 != 0 or ny % 2 != 0 or nz % 2 != 0:
raise ValueError('nx, ny and nz have to be even.')
self._op_fft = FFTW3DReal2Complex(nx, ny, nz)
self.ifft3d = self._op_fft.ifft3d
self.fft3d = self._op_fft.fft3d
kx_adim_max = nx/2
ky_adim_max = ny/2
kz_adim_max = nz/2
self.nkx = kx_adim_max + 1
self.nky = ny
self.nkz = nz
self.nk0 = self.nkz
self.nk1 = self.nky
self.nk2 = self.nkx
self.shape_fft = (self.nk0, self.nk1, self.nk2)
self.deltakx = 2*pi/Lx
self.deltaky = 2*pi/Ly
self.deltakz = 2*pi/Lz
self.k0 = self.deltakz * np.r_[0:kz_adim_max+1, -kz_adim_max+1:0]
self.k1 = self.deltakz * np.r_[0:ky_adim_max+1, -ky_adim_max+1:0]
self.k2 = self.deltakx * np.arange(self.nk2)
K1, K0, K2 = np.meshgrid(self.k1, self.k0, self.k2, copy=False)
self.Kz = K0
self.Ky = K1
self.Kx = K2
self.K_square_nozero = K0**2 + K1**2 + K2**2
self.K_square_nozero[0, 0, 0] = 1e-14
def project_perpk3d(self, vx_fft, vy_fft, vz_fft):
Kx = self.Kx
Ky = self.Ky
Kz = self.Kz
tmp = (Kx * vx_fft + Ky * vy_fft + Kz * vz_fft) / self.K_square_nozero
return (vx_fft - Kx * tmp, vy_fft - Ky * tmp, vz_fft - Kz * tmp)
def vgradv_from_v(self, vx, vy, vz, vx_fft=None, vy_fft=None, vz_fft=None):
ifft3d = self.ifft3d
if vx_fft is None:
vx_fft = self.fft3d(vx)
vy_fft = self.fft3d(vy)
vz_fft = self.fft3d(vz)
Kx = self.Kx
Ky = self.Ky
Kz = self.Kz
px_vx_fft = 1j * Kx * vx_fft
py_vx_fft = 1j * Ky * vx_fft
pz_vx_fft = 1j * Kz * vx_fft
px_vy_fft = 1j * Kx * vy_fft
py_vy_fft = 1j * Ky * vy_fft
pz_vy_fft = 1j * Kz * vy_fft
px_vz_fft = 1j * Kx * vz_fft
py_vz_fft = 1j * Ky * vz_fft
pz_vz_fft = 1j * Kz * vz_fft
vgradvx = (vx * ifft3d(px_vx_fft) +
vy * ifft3d(py_vx_fft) +
vz * ifft3d(pz_vx_fft))
vgradvy = (vx * ifft3d(px_vy_fft) +
vy * ifft3d(py_vy_fft) +
vz * ifft3d(pz_vy_fft))
vgradvz = (vx * ifft3d(px_vz_fft) +
vy * ifft3d(py_vz_fft) +
vz * ifft3d(pz_vz_fft))
return vgradvx, vgradvy, vgradvz
if __name__ == '__main__':
n = 4
from fluidsim.base.params import Parameters
p = Parameters(tag='params')
p._set_child(
'oper', {'nx': n, 'ny': n, 'nz': 2*n,
'Lx': 2*pi, 'Ly': 2*pi, 'Lz': 2*pi})
oper = OperatorsPseudoSpectral3D(params=p)
field = np.ones(oper.shape_phys)
field_fft = oper.fft3d(field)
assert field_fft.shape == oper.shape_fft
oper.vgradv_from_v(field, field, field)
oper.project_perpk3d(field_fft, field_fft, field_fft)
......@@ -17,6 +17,6 @@
from fluidsim.base.setofvariables import SetOfVariables
from fluidsim.base.solvers.pseudo_spect import (
SimulBasePseudoSpectral, InfoSolverPseudoSpectral)
SimulBasePseudoSpectral, InfoSolverPseudoSpectral3d)
......@@ -21,6 +21,6 @@
class InfoSolverNS3D(InfoSolverPseudoSpectral):
class InfoSolverNS3D(InfoSolverPseudoSpectral3d):
def _init_root(self):
super(InfoSolverNS3D, self)._init_root()
......@@ -100,8 +100,6 @@
"""This static method is used to complete the *params* container.
"""
SimulBasePseudoSpectral._complete_params_with_default(params)
# attribs = {'beta': 0.}
# params._set_attribs(attribs)
def tendencies_nonlin(self, state_fft=None):
oper = self.oper
......@@ -155,5 +153,5 @@
params.short_name_type_run = 'test'
nh = 32
n = 16
Lh = 2*np.pi
......@@ -159,5 +157,6 @@
Lh = 2*np.pi
params.oper.nx = nh
params.oper.ny = nh
params.oper.nx = n
params.oper.ny = n
params.oper.nz = n
params.oper.Lx = Lh
params.oper.Ly = Lh
......@@ -162,7 +161,6 @@
params.oper.Lx = Lh
params.oper.Ly = Lh
# params.oper.type_fft = 'FFTWPY'
params.oper.Lz = Lh
delta_x = params.oper.Lx/params.oper.nx
params.nu_8 = 2.*10e-1*params.forcing.forcing_rate**(1./3)*delta_x**8
......@@ -179,5 +177,5 @@
# params.output.periods_print.print_stdout = 0.25
params.output.periods_save.phys_fields = 1.
params.output.periods_save.spectra = 0.5
# params.output.periods_save.spectra = 0.5
params.output.periods_save.spatial_means = 0.05
......@@ -183,6 +181,5 @@
params.output.periods_save.spatial_means = 0.05
params.output.periods_save.spect_energy_budg = 0.5
params.output.periods_save.increments = 0.5
# params.output.periods_save.spect_energy_budg = 0.5
params.output.periods_plot.phys_fields = 0.0
......@@ -191,7 +188,6 @@
# params.output.spectra.HAS_TO_PLOT_SAVED = True
# params.output.spatial_means.HAS_TO_PLOT_SAVED = True
# params.output.spect_energy_budg.HAS_TO_PLOT_SAVED = True
# params.output.increments.HAS_TO_PLOT_SAVED = True
params.output.phys_fields.field_to_plot = 'rot'
......
"""State for the NS3D solver (:mod:`fluidsim.solvers.ns3d.state`)
=======================================================================
"""
from fluidsim.base.state import StatePseudoSpectral
from fluiddyn.util import mpi
class StateNS2D(StatePseudoSpectral):
"""Contains the variables corresponding to the state and handles the
access to other fields for the solver NS3D.
"""
@staticmethod
def _complete_info_solver(info_solver):
"""Complete the ParamContainer info_solver.
This is a static method!
"""
keys_state_phys = ['vx', 'vy', 'vz']
info_solver.classes.State._set_attribs({
'keys_state_fft': [k + '_fft' for k in keys_state_phys],
'keys_state_phys': keys_state_phys,
'keys_phys_needed': keys_state_phys,
'keys_linear_eigenmodes': ['rot_fft']})
def compute(self, key, SAVE_IN_DICT=True, RAISE_ERROR=True):
it = self.sim.time_stepping.it
if (key in self.vars_computed and it == self.it_computed[key]):
return self.vars_computed[key]
if key == 'vx_fft':
result = self.oper.fft3d(self.state_phys.get_var('vx'))
elif key == 'vy_fft':
result = self.oper.fft3d(self.state_phys.get_var('vy'))
elif key == 'vz_fft':
result = self.oper.fft3d(self.state_phys.get_var('vz'))
else:
to_print = 'Do not know how to compute "' + key + '".'
if RAISE_ERROR:
raise ValueError(to_print)
else:
if mpi.rank == 0:
print(to_print +
'\nreturn an array of zeros.')
result = self.oper.constant_arrayX(value=0.)
if SAVE_IN_DICT:
self.vars_computed[key] = result
self.it_computed[key] = it
return result
def statephys_from_statefft(self):
vx_fft = self.state_fft.get_var('vx_fft')
vy_fft = self.state_fft.get_var('vy_fft')
vz_fft = self.state_fft.get_var('vz_fft')
self.state_phys.set_var('vx', self.oper.ifft3d(vx_fft))
self.state_phys.set_var('vy', self.oper.ifft3d(vy_fft))
self.state_phys.set_var('vz', self.oper.ifft3d(vz_fft))
def statefft_from_statephys(self):
vx = self.state_phys.get_var('vx')
vy = self.state_phys.get_var('vy')
vz = self.state_phys.get_var('vz')
self.state_fft.get_var('vx_fft', self.oper.fft3d(vx))
self.state_fft.get_var('vy_fft', self.oper.fft3d(vy))
self.state_fft.get_var('vz_fft', self.oper.fft3d(vz))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment