From 883953390021bb34c3d8b8d98da581e1f3460a6b Mon Sep 17 00:00:00 2001 From: ignis Date: Mon, 19 Feb 2018 17:21:43 +0900 Subject: [PATCH] initial scalar generator commit --- gen_scalar_field.py | 137 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 137 insertions(+) create mode 100644 gen_scalar_field.py diff --git a/gen_scalar_field.py b/gen_scalar_field.py new file mode 100644 index 0000000..04c4d5e --- /dev/null +++ b/gen_scalar_field.py @@ -0,0 +1,137 @@ +# coding: utf-8 +class TopHat : + def __init__ (self, ks ,k0) : + self.ks = ks + self.k0 = k0 + + def __call__ (self, k) : + import numpy as np + f = np.ones (k.shape) + f[k < self.ks - self.k0/2] = 0 + f[k > self.ks + self.k0/2] = 0 + return f + +class CutOff : + def __init__ (self, kc) : + self.kc = kc + + def __call__ (self, k) : + import numpy as np + f = np.ones (k.shape) + f[k > self.kc] = self.kc**2 / k[k > self.kc]**2 + return f + +class RandomScalarField2D: + def __init__ (self, nx, lx, ny, ly, scalar_max, scalar_min, ksk0, kcks) : + import numpy as np + dx = lx / nx + dy = ly / ny + self.X, self.Y = np.meshgrid(np.arange(0,nx*dx,dx), np.arange(0,ny*dy,dy)) + self.KX, self.KY = np.meshgrid(2*np.pi*np.fft.fftfreq(nx, dx), 2*np.pi*np.fft.fftfreq(ny, dy)) + print (np.arange(0,nx*dx,dx), np.arange(0,ny*dy,dy)) + print (2*np.pi*np.fft.fftfreq(nx, dx), 2*np.pi*np.fft.fftfreq(ny, dy)) + self.k = np.sqrt(self.KX**2 + self.KY**2) + 1e-15 + self.smax = scalar_max + self.smin = scalar_min + + self.Nmax = max(nx, ny) + # self.k0 = self.k.max() / self.Nmax + self.k0 = 2 * np.pi / lx + print self.k0 + # self.k0 = min (np.fft.fftfreq(nx, dx)[1], np.fft.fftfreq(ny, dy)[1]) + # print self.k0 + self.ks = ksk0 * self.k0 + self.kc = kcks * self.ks + + print self.ks, self.kc + print self.k + + self.f = TopHat(self.ks, self.k0) + self.F = CutOff(self.kc) + + self.Phi = np.sqrt(self.f(self.k)/(4*np.pi*self.k**2)) * np.exp(2*np.pi*(1j)*np.random.uniform(0, 1, self.k.shape)) + + print np.count_nonzero(self.Phi) + print np.nonzero(self.Phi) + print self.Phi[np.nonzero(self.Phi)] + + self.iPhi = np.fft.ifft2 (self.Phi) + print self.iPhi + print np.angle(self.iPhi) + print np.angle(self.iPhi).min() + print np.angle(self.iPhi).max() + self.phi = np.zeros(self.iPhi.shape) + self.phi[np.angle(self.iPhi) >= 0 ] = self.smax + self.phi[np.angle(self.iPhi) < 0 ] = self.smin + + self.Phi_ = self.F(self.k) * np.fft.fft2(self.phi) + self.phi_ = np.fft.ifft2(self.Phi_) + + def energy_spectrum (self) : + import numpy as np + E = np.zeros (self.Nmax + 1) + for i in np.arange(self.Nmax + 1): + k = self.k[np.logical_and( + self.k >= (i-0.5)*self.k0 , self.k < (i+0.5)*self.k0 )] + Phi_n = self.f(k)/(4*np.pi*k**2) + E[i] = Phi_n.sum() + E /= E.sum() + E += 10e-8 + + E_ = np.zeros (self.Nmax + 1) + Phi_ = np.fft.fft2(self.phi) + for i in np.arange(self.Nmax + 1): + Phi_n = Phi_[np.logical_and( + self.k >= (i-0.5)*self.k0 , self.k < (i+0.5)*self.k0 )] + E_[i] = (np.abs(Phi_n)**2).sum() + E_ /= E_.sum() + E_ += 10e-8 + + E__ = np.zeros (self.Nmax + 1) + for i in np.arange(self.Nmax + 1): + Phi_n = self.Phi_[np.logical_and( + self.k >= (i-0.5)*self.k0 , self.k < (i+0.5)*self.k0 )] + E__[i] = (np.abs(Phi_n)**2).sum() + E__ /= E__.sum() + E__ += 10e-8 + + return E, E_, E__ + + +class RandomScalarField3D: + def __init__ (self, nx, dx, scalar_max, scalar_min, ksk0, kcks) : + import numpy as np + self.X, self.Y, self.Z = np.meshgrid(np.arange(0,nx*dx,dx), np.arange(0,nx*dx,dx), np.arange(0,nx*dx,dx)) + self.KX, self.KY, self.KZ = np.meshgrid(np.fft.fftfreq(nx, dx), np.fft.fftfreq(nx, dx), np.fft.fftfreq(nx, dx)) + + print (np.arange(0,nx*dx,dx)) + + print (np.fft.fftfreq(nx, dx)) + + self.k = np.sqrt(self.KX**2 + self.KY**2 + self.KZ**2) + 1e-15 + + self.smax = scalar_max + + self.smin = scalar_min + + self.Nmax = nx + + self.k0 = self.k.max() / self.Nmax + self.ks = ksk0 * self.k0 + self.kc = kcks * self.ks + + self.f = TopHat(self.ks, self.k0) + self.F = CutOff(self.kc) + + self.Phi = np.sqrt(self.f(self.k)/(4*np.pi*self.k**2)) * np.exp(2*np.pi*(1j)*np.random.uniform(0, 1, self.k.shape)) + + self.iPhi = np.fft.ifftn (self.Phi) + self.phi = np.zeros(self.iPhi.shape) + self.phi[self.iPhi >= 0 ] = self.smax + self.phi[self.iPhi < 0 ] = self.smin + + self.Phi_ = self.F(self.k) * np.fft.fftn(self.phi) + self.phi_ = np.fft.ifftn(self.Phi_) + ''' + ''' +