test_turtle_NSVEparticles.py

#! /usr/bin/env python
#######################################################################
#                                                                     #
#  Copyright 2019 Max Planck Institute                                #
#                 for Dynamics and Self-Organization                  #
#                                                                     #
#  This file is part of TurTLE.                                         #
#                                                                     #
#  TurTLE is free software: you can redistribute it and/or modify       #
#  it under the terms of the GNU General Public License as published  #
#  by the Free Software Foundation, either version 3 of the License,  #
#  or (at your option) any later version.                             #
#                                                                     #
#  TurTLE is distributed in the hope that it will be useful,            #
#  but WITHOUT ANY WARRANTY; without even the implied warranty of     #
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      #
#  GNU General Public License for more details.                       #
#                                                                     #
#  You should have received a copy of the GNU General Public License  #
#  along with TurTLE.  If not, see <http://www.gnu.org/licenses/>       #
#                                                                     #
# Contact: Cristian.Lalescu@ds.mpg.de                                 #
#                                                                     #
#######################################################################



import os
import numpy as np
import h5py
import sys

import TurTLE
from TurTLE import DNS


def main():
    """
        Run short DNS, test against archived data.

        This test runs two MPI jobs, generating a short DNS from initial
        conditions stored as python package data.
        This is a full DNS of Navier-Stokes with particles, and we check that
        the new fields, trajectories and velocity sampled along trajectories
        agree to within numerical precision (single precision).

        Implicitly, this checks that the entire compilation mechanism is set
        up properly, as well as the correctness of the code.
        4 MPI processes with one OpenMP thread per process are used by default,
        but this may be changed with command line parameters as desired.
    """
    niterations = 32
    nparticles = 10000
    njobs = 2
    c = DNS()
    c.launch(
            ['NSVEparticles',
             '-n', '32',
             '--src-simname', 'B32p1e4',
             '--forcing_type', 'linear',
             '--src-wd', TurTLE.lib_dir + '/test',
             '--src-iteration', '0',
             '--simname', 'dns_nsveparticles',
             '--np', '4',
             '--ntpp', '1',
             '--fftw_plan_rigor', 'FFTW_PATIENT',
             '--niter_todo', '{0}'.format(niterations),
             '--niter_out', '{0}'.format(niterations),
             '--niter_stat', '1',
             '--checkpoints_per_file', '{0}'.format(3),
             '--nparticles', '{0}'.format(nparticles),
             '--particle-rand-seed', '2',
             '--njobs', '{0}'.format(njobs),
             '--wd', './'] +
             sys.argv[1:])
    f0 = h5py.File(
            os.path.join(
                os.path.join(TurTLE.lib_dir, 'test'),
                'B32p1e4_checkpoint_0.h5'),
            'r')
    f1 = h5py.File(c.get_checkpoint_0_fname(), 'r')
    for iteration in [0, 32, 64]:
        field0 = f0['vorticity/complex/{0}'.format(iteration)][...]
        field1 = f1['vorticity/complex/{0}'.format(iteration)][...]
        field_error = np.max(np.abs(field0 - field1))
        x0 = f0['tracers0/state/{0}'.format(iteration)][...]
        x1 = f1['tracers0/state/{0}'.format(iteration)][...]
        traj_error = np.max(np.abs(x0 - x1))
        y0 = f0['tracers0/rhs/{0}'.format(iteration)][...]
        y1 = f1['tracers0/rhs/{0}'.format(iteration)][...]
        rhs_error = np.max(np.abs(y0 - y1))
        assert(field_error < 1e-5)
        assert(traj_error < 1e-5)
        assert(rhs_error < 1e-5)
    print('SUCCESS! Basic test passed.')
    return None

if __name__ == '__main__':
    main()