cpu_plane_jet

examples/cpu_plane_jet/cpu_planejet.py

+#!/usr/bin/env python
+
+import hysop
+
+from hysop.methods_keys import TimeIntegrator, Interpolation, Remesh,\
+    Support, Splitting, MultiScale, MultiScale, SpaceDiscretisation, \
+    GhostUpdate, Scales, dtCrit, ExtraArgs
+
+from hysop.problem.simulation        import Simulation
+from hysop.fields.variable_parameter import VariableParameter
+
+from hysop.numerics.odesolvers import RK2 as RK2
+from hysop.numerics.odesolvers import RK3 as RK3
+from hysop.numerics.odesolvers import RK4 as RK4
+
+from hysop.numerics.finite_differences import FD_C_4
+from hysop.numerics.remeshing import L6_4 as remesh_formula
+from hysop.numerics.remeshing import Linear
+
+from hysop.operator.advection import Advection
+from hysop.operator.diffusion import Diffusion
+from hysop.operator.stretching import Stretching
+from hysop.operator.poisson import Poisson
+from hysop.operator.adapt_timestep import AdaptTimeStep
+
+from hysop.operator.redistribute_inter import RedistributeInter
+from hysop.operator.redistribute_intra import RedistributeIntra
+
+from hysop.gpu.gpu_transfer import DataTransfer
+from hysop.operator.hdf_io  import HDF_Writer
+
+import utilities as utils
+
+
+## Simulation parameters
+dim = 3
+N   = 1<<8 + 1
+nb_elems = [N]*dim
+
+## Physical parameters:
+# Flow viscosity
+viscosity = 1e-4
+
+## I/O parameters
+dt_output    = IOParams(frequency=out_freq, filename='cpu_dt.dat',    fileformat=ASCII)
+velo_output  = IOParams(frequency=out_freq, filename='cpu_velo.dat',  fileformat=HDF5)
+vorti_output = IOParams(frequency=out_freq, filename='cpu_vorti.dat', fileformat=HDF5)
+
+
+# Domain
+box = hysop.Box(length=[1.0]*dim, 
+                origin=[0.0]*dim,
+                bc=PERIODIC)
+
+# Fields
+velo  = hysop.Field(domain=box, formula=utils.initialize_velocity,
+                   name='Velocity',  is_vector=True)
+vorti = hysop.Field(domain=box,
+                   name='Vorticity', is_vector=True)
+
+# Discretizations
+d_uw = Discretization(nb_elems)
+
+# Topologies
+topo_uw = box.create_topology(d_uw, dim=1)
+
+# Adaptative timestep simulation
+dt   = VariableParameter({'dt': 0.001})
+simu = Simulation(start=0.0, end=5.0, time_step=0.001, max_iter=10000)
+
+
+## Operators
+ops = {}
+
+# CPU operators
+ops['advec_w'] = Advection(velo,
+                     discretization=topo_uw,
+                     variables={vorti: topo_uw},
+                     method={Scales: 'p_64', MultiScale: 'L4_4'})
+ops['stretching'] = Stretching(velo, vorti, discretization=d_uw)
+ops['diffusion']  = Diffusion(variables={vorti: d_uw}, viscosity=viscosity)
+ops['poisson']    = Poisson(velo, vorti, discretization=d_uw)
+
+
+# other operators
+ops['dt_adapt'] = AdaptTimeStep(velo, vorti,
+                         simulation=simu,
+                         time_range=[0, np.infty],
+                         discretization=d_uw,
+                         method={TimeIntegrator: RK3,
+                                 SpaceDiscretisation: FD_C_4,
+                                 dtCrit: ['gradU', 'cfl']},
+                         lcfl=0.15,
+                         cfl=1.5,
+                         io_params=dt_output)
+
+ops['velo_io']  = HDF_Writer(variables={velo : topo_uw}, io_params=velo_output)
+ops['vorti_io'] = HDF_Writer(variables={vorti: topo_uw}, io_params=vorti_output)
+
+for op in ops:
+    op.discretize()
+
+

examples/cpu_plane_jet/create_random_arrays.py

+import os
+import numpy as np
+from hysop.constants import HYSOP_REAL, ORDER
+
+
+def random_init(shape, mpi_comm):
+    # Create a folder to store all random arrays
+    d = 'rand_init'
+    if mpi_comm.Get_rank() == 0:
+        if not os.path.exists(d):
+            os.makedirs(d)
+    mpi_comm.Barrier()
+    file_name = "{0}_{1}_{2}".format(*shape)
+    file_name += "_{0}p_{1}.dat".format(mpi_comm.Get_size(),
+                                        mpi_comm.Get_rank())
+    try:
+        randX = np.asarray(
+            np.reshape(np.fromfile(os.path.join(d, 'randX_' + file_name)),
+                       shape),
+            dtype=HYSOP_REAL, order=ORDER)
+        randY = np.asarray(
+            np.reshape(np.fromfile(os.path.join(d, 'randY_' + file_name)),
+                       shape),
+            dtype=HYSOP_REAL, order=ORDER)
+        randZ = np.asarray(
+            np.reshape(np.fromfile(os.path.join(d, 'randZ_' + file_name)),
+                       shape),
+            dtype=HYSOP_REAL, order=ORDER)
+    except IOError:
+        randX = np.asarray(np.random.random(shape),
+                           dtype=HYSOP_REAL, order=ORDER) - 0.5
+        randY = np.asarray(np.random.random(shape),
+                           dtype=HYSOP_REAL, order=ORDER) - 0.5
+        randZ = np.asarray(np.random.random(shape),
+                           dtype=HYSOP_REAL, order=ORDER) - 0.5
+        randX.tofile(os.path.join(d, 'randX_' + file_name))
+        randY.tofile(os.path.join(d, 'randY_' + file_name))
+        randZ.tofile(os.path.join(d, 'randZ_' + file_name))
+    return randX, randY, randZ

examples/cpu_plane_jet/utilities.py

+
+import numpy as np
+
+width = 0.01
+ampl1 = 0.05
+ampl3 = 0.3
+
+
+def initialize_velocity(res, x, y, z, t, rand_init=False):
+    yy  = np.abs(y - 0.5)
+    aux = (0.1 - 2.0 * yy) / (4.0 * width)
+    strg = np.exp(np.abs(aux**2))
+    if rand_init:
+        from create_random_arrays import random_init
+        randX, randY, randZ = random_init(res[0].shape, box.comm_task)
+        strg1 = strg * randX
+        strg2 = strg * randY
+        strg3 = strg * randZ
+        res[0][...] = 0.5 * (1.0 + np.tanh(aux)) * (1.0 + ampl3 * np.sin(8.0 * np.pi * x)) * (1.0 + ampl1 * strg1)
+        res[1][...] = ampl1 * strg2
+        res[2][...] = ampl1 * strg3
+    else:
+        res[0][...] = 0.5 * (1.0 + np.tanh(aux)) * (1.0 + ampl3 * np.sin(8.0 * np.pi * x))
+        res[1][...] = 0.0
+        res[2][...] = 0.0
+    return res
+