diff --git a/.github/workflows/dependencies/documentation.sh b/.github/workflows/dependencies/documentation.sh index 99696fbc6..dc66ccaee 100755 --- a/.github/workflows/dependencies/documentation.sh +++ b/.github/workflows/dependencies/documentation.sh @@ -9,13 +9,8 @@ set -eu -o pipefail sudo apt-get update +# The docs job only builds HTML with Sphinx; doxygen and LaTeX are not used. sudo apt-get install -y --no-install-recommends\ build-essential \ - pandoc \ - doxygen \ - texlive \ - texlive-latex-extra \ - texlive-lang-cjk \ - tex-gyre \ - latexmk + pandoc diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 7725dd1a3..23f578231 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -6,7 +6,7 @@ on: pull_request: concurrency: - group: ${{ github.head_ref }}-docs + group: ${{ github.ref }}-${{ github.head_ref }}-docs cancel-in-progress: true jobs: @@ -21,7 +21,7 @@ jobs: .github/workflows/dependencies/documentation.sh echo "Installing python packages for docs..." python3 -m pip install --upgrade pip - python3 -m pip install sphinx sphinx_rtd_theme breathe sphinxcontrib.bibtex docutils + python3 -m pip install sphinx sphinx_rtd_theme sphinxcontrib.bibtex docutils - name: Install and Build run: | diff --git a/.gitignore b/.gitignore index f52ead781..d704d4a98 100644 --- a/.gitignore +++ b/.gitignore @@ -7,5 +7,7 @@ plt????? d/ f/ o/ +tmp_build_dir/ +tmp_install_dir/ chk?????.old* plt?????.old* diff --git a/Docs/sphinx_documentation/Makefile b/Docs/sphinx_documentation/Makefile index 078ec6c4f..568ace415 100644 --- a/Docs/sphinx_documentation/Makefile +++ b/Docs/sphinx_documentation/Makefile @@ -16,9 +16,7 @@ help: .PHONY: help Makefile clean: - rm -rf build ../Doxygen/xml source/class source/file - rm -f source/classlist.rst source/filelist.rst - rm -rf source/*_files.rst + rm -rf build # Catch-all target: route all unknown targets to Sphinx using the new # "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). diff --git a/Docs/sphinx_documentation/make_api.py b/Docs/sphinx_documentation/make_api.py deleted file mode 100644 index 273626961..000000000 --- a/Docs/sphinx_documentation/make_api.py +++ /dev/null @@ -1,79 +0,0 @@ -""" -This script rewrites the filelist.rst file autogenerated by breathe to organize the -files in the documentation's API into the same directory structure as their -counterparts in the ../../Src directory. - -We are not going to include cpp files - doxygen assumes that functions are documented -in the header files, so that's what we'll assume here as well. - -TODO: check if doxygen has generated the corresponding files and only print their -names if so -TODO: don't make a _files.rst file if the subdir contains no doxygen-erated files -""" - -import os -import re - - -def generate_filelist(rootdir, outfile, output_data, subdir_prefix=""): - for subdir in sorted(os.listdir(rootdir)): - if not os.path.isdir(os.path.join(rootdir, subdir)): - # found a file - if (subdir[-4:] != ".F90" and subdir[-4:] != ".f90" and subdir[-2:] != ".H" or subdir[-3:] == "F.H"): - continue - - rst_name = re.sub("_", "__", subdir) - rst_name = re.sub("\.", "_8", rst_name) - - output_data += """file/{} - """.format(rst_name) - - else: - # found a subdirectory - create a new _files.rst file and call - # generate_filelist on the subdir - - # ignore F_Interfaces - if subdir.lower() in ['f_interfaces']: - continue - - output_data += """{}_files - """.format(subdir_prefix + subdir) - - subdir_file_name = "source/{}_files.rst".format( - subdir_prefix + subdir) - - with open(subdir_file_name, 'w') as subdir_file: - - subdir_output_data = "{}\n".format(subdir.capitalize()) - subdir_output_data += "=" * len(subdir) - subdir_output_data += """ - -.. toctree:: - :maxdepth: 3 - - """ - generate_filelist(os.path.join(rootdir, subdir), subdir_file, - subdir_output_data, - subdir_prefix=subdir_prefix + subdir + '_') - - outfile.write(output_data) - - -if __name__ == "__main__": - - # directory of the source files - rootdir = "../../Src" - - outfile_path = "source/filelist.rst" - - with open(outfile_path, 'w') as outfile: - - output_data = """File list -========= - -.. toctree:: - :maxdepth: 3 - - """ - - generate_filelist(rootdir, outfile, output_data) diff --git a/Docs/sphinx_documentation/source/AlgorithmOptions.rst b/Docs/sphinx_documentation/source/AlgorithmOptions.rst index 76bbdfb40..d48536a6e 100644 --- a/Docs/sphinx_documentation/source/AlgorithmOptions.rst +++ b/Docs/sphinx_documentation/source/AlgorithmOptions.rst @@ -27,23 +27,27 @@ Note that Temperature is only non-conservative. For more details, see :ref:`sec: Advection --------- -IAMR has the option to use a Method of Lines (MOL) or Godunov scheme to compute the advective terms. +IAMR computes the advective terms with an unsplit Godunov scheme (piecewise linear or piecewise +parabolic reconstruction) or with the Bell-Dawson-Shubin (BDS) scheme. The following must be +preceded by "ns." +-------------------------+-------------------------------------------------------------------------+-------------+--------------+ | | Description | Type | Default | +=========================+=========================================================================+=============+==============+ -| ns.use_godunov | If true, use Godunov, else use MOL. | bool | true | +| advection_scheme | Godunov_PLM, Godunov_PPM or BDS. Godunov_PPM and BDS are not | String | Godunov_PLM | +| | available with embedded boundaries. | | | +-------------------------+-------------------------------------------------------------------------+-------------+--------------+ +Note that the old ``ns.use_godunov`` key and the MOL scheme have been removed; setting either +aborts the run. -For problems without embedded boundaries, there are additional options when using the Godunov method. The following must + +For problems without embedded boundaries, there is an additional option for the Godunov method. The following must be preceded by "godunov." +-------------------------+-------------------------------------------------------------------------+-------------+--------------+ | | Description | Type | Default | +=========================+=========================================================================+=============+==============+ -| use_ppm | Use the Piecewise Parabolic Method to construct edge states | bool | false | -+-------------------------+-------------------------------------------------------------------------+-------------+--------------+ | use_forces_in_trans | Use external forcing terms in constructing transverse derivatives | bool | false | +-------------------------+-------------------------------------------------------------------------+-------------+--------------+ @@ -60,3 +64,56 @@ The following must be preceded by "ns." +-------------------------+-----------------------------------------------------------------------+-------------+--------------+ Note the default value of ``ns.be_cn_theta = 0.5`` corresponds to the Crank-Nicolson method. + + +.. _sec:LES: + +Large Eddy Simulation +--------------------- + +IAMR can add a subgrid-scale eddy viscosity to the viscous terms. The following must be +preceded by "ns." + ++-------------------------+-----------------------------------------------------------------------+-------------+--------------+ +| | Description | Type | Default | ++=========================+=======================================================================+=============+==============+ +| do_LES | Add a subgrid-scale eddy viscosity to the molecular viscosity | Int | 0 | ++-------------------------+-----------------------------------------------------------------------+-------------+--------------+ +| LES_model | Which model to use: Smagorinsky or Sigma. Any other value aborts. | String | Smagorinsky | +| | Sigma is 3D only. | | | ++-------------------------+-----------------------------------------------------------------------+-------------+--------------+ +| smago_Cs_cst | Model constant, used only when LES_model = Smagorinsky | Real | 0.18 | ++-------------------------+-----------------------------------------------------------------------+-------------+--------------+ +| sigma_Cs_cst | Model constant, used only when LES_model = Sigma | Real | 1.5 | ++-------------------------+-----------------------------------------------------------------------+-------------+--------------+ +| getLESVerbose | Print the model and constant in use from the LES routine | Int | 0 | ++-------------------------+-----------------------------------------------------------------------+-------------+--------------+ + +Note that each model reads its own constant, so changing ``ns.smago_Cs_cst`` has no effect +when ``ns.LES_model = Sigma``, and vice versa. The Sigma model is described in +Nicoud et al., *Using singular values to build a subgrid-scale model for large eddy +simulations*, Phys. Fluids 23, 085106 (2011). + + +.. _sec:EBOptions: + +Embedded Boundaries +------------------- + +These apply to builds with ``USE_EB=TRUE``; see :ref:`sec:EB-basics` for how the geometry +itself is constructed. The following must be preceded by "ns." + ++-------------------------+-----------------------------------------------------------------------+-------------+--------------+ +| | Description | Type | Default | ++=========================+=======================================================================+=============+==============+ +| redistribution_type | How the advective update of a small cut cell is redistributed to its | String | StateRedist | +| | neighbours: StateRedist, FluxRedist or NoRedist. Any other value | | | +| | aborts. | | | ++-------------------------+-----------------------------------------------------------------------+-------------+--------------+ +| refine_cutcells | Tag every cut cell for refinement, so that the embedded boundary | Int | 1 | +| | never crosses a coarse/fine boundary. Setting 0 allows a partially | | | +| | refined EB, which is still under development and issues a warning. | | | ++-------------------------+-----------------------------------------------------------------------+-------------+--------------+ + +Note that ``ns.advection_scheme = Godunov_PPM`` and ``BDS`` are not available with embedded +boundaries. diff --git a/Docs/sphinx_documentation/source/Contributing.rst b/Docs/sphinx_documentation/source/Contributing.rst index 7daef838c..190e1a4d1 100644 --- a/Docs/sphinx_documentation/source/Contributing.rst +++ b/Docs/sphinx_documentation/source/Contributing.rst @@ -65,7 +65,7 @@ Make your own fork ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ First, setup your local git repo. To make your own fork of the main -(`upstream`) repository, press the fork button on the `IAMR Github page `_. +(`upstream`) repository, press the fork button on the `IAMR Github page `_. Then, clone your fork on your local computer. If you plan on doing a lot of IAMR development, we recommend configuring your clone to use ssh access so you won't have to enter your Github @@ -77,7 +77,7 @@ password every time, which you can do using these commands: # Then, navigate into your repo, add a new remote for the main IAMR repo, and fetch it: cd IAMR - git remote add upstream https://github.com/IAMR-Codes/IAMR + git remote add upstream https://github.com/AMReX-Fluids/IAMR git remote set-url --push upstream git@github.com:/IAMR.git git fetch upstream @@ -95,7 +95,7 @@ If you instead prefer to use HTTPS authentication, configure your local clone as # Navigate into your repo, add a new remote for the main IAMR repo, and fetch it cd IAMR - git remote add upstream https://github.com/IAMR-Codes/IAMR + git remote add upstream https://github.com/AMReX-Fluids/IAMR git remote set-url --push upstream https://github.com//IAMR.git git fetch upstream diff --git a/Docs/sphinx_documentation/source/Introduction_Chapter.rst b/Docs/sphinx_documentation/source/Introduction_Chapter.rst index 43628616d..aac213807 100644 --- a/Docs/sphinx_documentation/source/Introduction_Chapter.rst +++ b/Docs/sphinx_documentation/source/Introduction_Chapter.rst @@ -16,7 +16,7 @@ Key software and algorithmic features of IAMR include: * Fluid velocity, density and tracers are defined at cell centroids; pressure is defined at nodes. -* Possible advection algorithms: a Method-Of-Lines (MOL) approach and a Godunov-method algorithm. Both use an intermediate MAC projection for face-centered advection velocities. +* Possible advection algorithms: an unsplit Godunov method (piecewise linear or piecewise parabolic reconstruction) and the Bell-Dawson-Shubin (BDS) scheme. All use an intermediate MAC projection for face-centered advection velocities. * Incompressibility of the fluid is imposed through the use of an approximate projection at the end of the time step. diff --git a/Docs/sphinx_documentation/source/ProblemSetup.rst b/Docs/sphinx_documentation/source/ProblemSetup.rst index 8cae7a0b1..b52a4c011 100644 --- a/Docs/sphinx_documentation/source/ProblemSetup.rst +++ b/Docs/sphinx_documentation/source/ProblemSetup.rst @@ -36,9 +36,85 @@ and problem parameters from the inputs file, and initializes the state data (velocity, density, etc.). The easiest way to get started is to create a new ``probtype`` by copying the code for an existing problem and modifying it to apply the desired initial state. +IAMR's make system puts the build directory ahead of ``IAMR/Source`` in its search path, so +a problem directory may override any IAMR source file simply by placing its own copy of that +file alongside its ``GNUmakefile``. This is the mechanism to use for a problem that needs its +own initial conditions (``prob_init.cpp``) or its own forcing function (``NS_getForce.cpp``). + +The problem is selected with ``prob.probtype``. The values recognized by the default +``Source/prob/prob_init.cpp`` are + ++------------+------------------------------------------------------------------------+-----------------------------+ +| probtype | Initial state | Example deck | ++============+========================================================================+=============================+ +| 1 | Fluid at rest with constant density | LidDrivenCavity | ++------------+------------------------------------------------------------------------+-----------------------------+ +| 2 | Density blob (bubble or drop) in a constant background | Bubble | ++------------+------------------------------------------------------------------------+-----------------------------+ +| 3 | Smoothed density/tracer jump normal to x, plus a tracer blob | -- | ++------------+------------------------------------------------------------------------+-----------------------------+ +| 4 | Constant velocity and density with a tracer blob | TracerAdvection, Poiseuille,| +| | | FlowPastCylinder | ++------------+------------------------------------------------------------------------+-----------------------------+ +| 5 | Double shear layer | DoubleShearLayer, Particles | ++------------+------------------------------------------------------------------------+-----------------------------+ +| 6 | As probtype 2, but with temperature as a state variable | HotSpot | ++------------+------------------------------------------------------------------------+-----------------------------+ +| 7 | Euler vortex tube in a triply periodic domain | Exec/run3d/regtest.3d.euler | ++------------+------------------------------------------------------------------------+-----------------------------+ +| 8 | Convected vortex | ConvectedVortex | ++------------+------------------------------------------------------------------------+-----------------------------+ +| 10 | Rayleigh-Taylor instability | RayleighTaylor | ++------------+------------------------------------------------------------------------+-----------------------------+ +| 11 | Taylor-Green vortex | TaylorGreen | ++------------+------------------------------------------------------------------------+-----------------------------+ + +Any other value aborts. Note that 9 is not used, and that a problem directory that supplies +its own ``prob_init.cpp``, as described above, defines its own set of values. + +The initial state is parameterized by the following, which must be preceded by "prob." +Which of them a given ``probtype`` actually reads varies; see ``prob_init.cpp``. + ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| | Description | Type | Default | ++========================+====================================================================+=============+===========+ +| density_ic | Background density | Real | 1.0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| velocity_ic | Background velocity, one value per coordinate direction | Reals | 0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| direction | Coordinate direction used by direction-dependent setups | Int | 0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| interface_width | Width over which density and tracer jumps are smoothed | Real | 1.0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| blob_center | Centre of the density/tracer blob, one value per direction | Reals | 0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| blob_radius | Radius of the density/tracer blob | Real | 0.1 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| rho_1, rho_2 | Densities on either side of the interface (probtypes 3 and 10) | Real | 1.0, 2.0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| tra_1, tra_2 | Tracer values on either side of the interface (probtypes 3 and 10) | Real | 0.0, 1.0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| perturbation_amplitude | Amplitude of the interface perturbation (probtype 10) | Real | 1.0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| velocity_factor | Another name for the first component of velocity_ic; it is the | Real | 0.0 | +| | velocity scale V_0 of the Taylor-Green vortex (probtype 11) | | | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| a, b, c | Problem-specific coefficients (probtype 11) | Real | 1.0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| xvort, yvort, rvort | Centre and radius of the convected vortex (probtype 8). These | Real | 0.5, 0.5, | +| | override a, b and c, which probtype 8 resets before reading them. | | 0.07 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| forcevort | Strength of the convected vortex (probtype 8) | Real | 6.0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| meanFlowDir | Direction of the mean flow the vortex is convected by (probtype 8) | Int | 0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ +| meanFlowMag | Magnitude of that mean flow (probtype 8) | Real | 0.0 | ++------------------------+--------------------------------------------------------------------+-------------+-----------+ + It is also possible to initialize the velocity using data from a previously generated plotfile. -To do this you must set ``BL_USE_VELOCITY=TRUE`` in the makefile, and provide the plotfile in -the inputs file via ``ns.velocity_plotfile = my_plotfile_name``. +To do this you must set ``USE_VELOCITY=TRUE`` in the makefile (which in turn defines the +``BL_USE_VELOCITY`` compiler macro), and provide the plotfile in the inputs file via +``ns.velocity_plotfile = my_plotfile_name``. Resolution @@ -223,7 +299,8 @@ preceded by “xlo”, “xhi”, “ylo”, “yhi”, “zlo”, and “zhi” +--------------------+---------------------------------------------------------------------------+-------------+-----------+ | temp | Sets temperature for mass inflows | Real | None | +--------------------+---------------------------------------------------------------------------+-------------+-----------+ -| pressure | Sets boundary pressure for pressure inflows, outflows and mass inflows | Real | None | +| pressure | Read only for pressure outflows, where it must be 0 (a non-zero value | Real | None | +| | is not yet implemented). Pressure inflows are not yet implemented. | | | +--------------------+---------------------------------------------------------------------------+-------------+-----------+ diff --git a/Docs/sphinx_documentation/source/RunningProblems.rst b/Docs/sphinx_documentation/source/RunningProblems.rst index 29d93500b..80c307a91 100644 --- a/Docs/sphinx_documentation/source/RunningProblems.rst +++ b/Docs/sphinx_documentation/source/RunningProblems.rst @@ -25,9 +25,16 @@ past it. +----------------------+-----------------------------------------------------------------------+-------------+--------------+ | stop_time | Maximum time to reach | Real | -1.0 | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ -| stop_when_steady | Stop when steady state is reached | Bool | false | +| num_steps | If > 0, take at most this many level-0 steps beyond the step count | Int | -1 | +| | the run starts from. Combined with max_step by taking the smaller | | | +| | of the two. Useful for advancing a restarted run by a fixed number | | | +| | of steps without having to look up its step count. | | | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ -| steady_tol | Specify tolerance to define steady state | Real | 1e-10 | +| strt_time | Simulation time at which to start (ignored when restarting). | Real | 0.0 | +| | Must be non-negative. | | | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| stop_interval | If > 0, stop this much simulation time after the time the run starts | Real | 0.0 | +| | from, i.e. stop_time becomes the restart time plus stop_interval. | | | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ The inputs below must be preceded by "ns." @@ -37,13 +44,18 @@ The inputs below must be preceded by "ns." +======================+=======================================================================+=============+==============+ | fixed_dt | Value of fixed dt if > 0 | Real | -1. | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ -| cfl | CFL constraint (dt < cfl * dx / u) if fixed_dt not > 0 | Real | 0.5 | +| cfl | CFL constraint (dt < cfl * dx / u); must be set if fixed_dt not > 0 | Real | none | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ | init_shrink | Factor by which to shrink the initial time step | Real | 1.0 | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ -| max_change | Factor by which time step can grow in subsequent steps | Real | 1.1 | +| change_max | Factor by which time step can grow in subsequent steps | Real | 1.1 | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ -| dt_cutoff | Time step below which the simulation will abort | Real | 0.0 | +| dt_cutoff | Time step at or below which the simulation stops (writing its final | Real | 0.0 | +| | checkpoint and plotfile) | | | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| stop_when_steady | Stop when steady state is reached | Bool | false | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| steady_tol | Specify tolerance to define steady state | Real | 1e-10 | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ * If you want to fix the dt, simply set :cpp:`ns.fixed_dt = XXX` to set the fluid time @@ -54,10 +66,10 @@ The inputs below must be preceded by "ns." condition, then set :cpp:`ns.cfl = 0.7` for example, and the fluid time step will be computed to be dt = 0.7 * dx / max(vel). - * Note that the cfl defaults to 0.5 so it does not have to be set in the inputs file. If neither - :cpp:`ns.cfl` nor :cpp:`fixed_dt` is set, then default value of cfl will be used. - If :cpp:`ns.fixed_dt` is set, then it will override the cfl option whether - :cpp:`ns.cfl` is set or not. + * Note that :cpp:`ns.cfl` has no default and must always be set in the inputs file, even + when :cpp:`ns.fixed_dt` is used. If :cpp:`ns.fixed_dt` is set, then it overrides the cfl + option: the time step is :cpp:`fixed_dt` (reduced by :cpp:`init_shrink` and then ramped + back up by :cpp:`change_max`) and the advective CFL estimate is not used to limit it. As an example, consider: @@ -169,6 +181,68 @@ and 61 level-0 steps is the first point when simulation time :math:`>=0.2`, etc. +.. _sec:InputsDerived: + +Derived Quantities +~~~~~~~~~~~~~~~~~~ + +The names below may be listed in ``amr.derive_plot_vars`` and used as the ``field_name`` +of a refinement indicator (see :ref:`sec:tagging`). They are registered in ``derive_lst`` +in ``Source/NS_setup.cpp``; state variables themselves (``x_velocity``, ``density``, +``tracer``, ...) are selected with ``amr.plot_vars`` instead. + ++----------------------+-----------------------------------------------------------------------+------------------------------+ +| Name | Description | Available when | ++======================+=======================================================================+==============================+ +| mag_vort | Magnitude of the vorticity | always | ++----------------------+-----------------------------------------------------------------------+------------------------------+ +| energy | Kinetic energy per unit volume, 0.5 rho U.U | always | ++----------------------+-----------------------------------------------------------------------+------------------------------+ +| avg_pressure | Nodal pressure averaged onto cell centres | always | ++----------------------+-----------------------------------------------------------------------+------------------------------+ +| velocity_average | Running time average and RMS fluctuation of each velocity component, | ns.avg_interval > 0 | +| | i.e. x_vel_average, ..., x_vel_rms, ... See | | +| | :ref:`sec:InputsTimeAverage`. | | ++----------------------+-----------------------------------------------------------------------+------------------------------+ +| particle_count | Number of particles in each cell of this level | USE_PARTICLES=TRUE | ++----------------------+-----------------------------------------------------------------------+------------------------------+ +| total_particle_count | Number of particles in each cell of this level and all finer levels | USE_PARTICLES=TRUE | ++----------------------+-----------------------------------------------------------------------+------------------------------+ + + +.. _sec:InputsTimeAverage: + +Time-Averaged Velocity +~~~~~~~~~~~~~~~~~~~~~~ + +IAMR can accumulate a running time average of the velocity field as the simulation +proceeds, and write it into each plotfile via the ``velocity_average`` derived quantity. +The following must be preceded by "ns." + ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| | Description | Type | Default | ++======================+=======================================================================+=============+==============+ +| avg_interval | Accumulate the average every this many level-0 steps. If <= 0, the | Int | 0 | +| | feature is off and the velocity_average derive is not registered. | | | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| compute_fluctuations | Also accumulate the RMS velocity fluctuation about the running mean. | Int | 0 | +| | Best enabled only once the mean itself has converged. | | | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ + +For example + +:: + + ns.avg_interval = 10 + ns.compute_fluctuations = 1 + amr.derive_plot_vars = velocity_average + +Note that the averages are carried in their own state type, so a checkpoint written +before this feature was enabled does not contain them; see ``ns.avg_in_checkpoint`` +in :ref:`sec:InputsCheckpoint`. + + + .. _sec:InputsCheckpoint: Checkpointing and Restarting @@ -199,6 +273,24 @@ The following inputs must be preceded by "amr." and control checkpoint/restart. | | If -1, number of files is always set equal to number of processors. | | | +-------------------------+-----------------------------------------------------------------------+-------------+-----------+ +The following must be preceded by "ns." and are needed only when restarting from a +checkpoint that was written by a build with a different set of state types. + ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| | Description | Type | Default | ++======================+=======================================================================+=============+==============+ +| gradp_in_checkpoint | Is the pressure gradient state present in the checkpoint being | Int | -1 | +| | restarted from? 1 = yes, 0 = no. | | | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| avg_in_checkpoint | Is the time-average state present in the checkpoint being restarted | Int | -1 | +| | from? 1 = yes, 0 = no. | | | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ + +If the checkpoint turns out to be missing a state type, IAMR aborts and asks for both +keys to be set. When in doubt, set both to 0. In particular, when turning on time +averaging (:ref:`sec:InputsTimeAverage`) for a run restarted from a checkpoint written +without it, set ``ns.avg_in_checkpoint = 0`` and ``ns.gradp_in_checkpoint = 1``. + Note: * ``amr.check_per`` will write a checkpoint at the first timestep whose ending time is past an integer multiple of this interval. @@ -247,7 +339,7 @@ The “derived quantity” ``particle_count`` represents the number of particles To visualize the particle locations as represented on the grid, add ``particle_count`` to the list of derived quanties in ``amr.derive_plot_vars =`` in the inputs file. -If ``particles.write_in_plotfile = 1`` in the inputs file, +If ``particles.particles_in_plotfile = 1`` in the inputs file, then the particle positions and velocities will be written in a binary file in each plotfile directory. This allows the use of the AMReX tools such as the particle comparison tool found in ``amrex/Tools/Postprocessing/C_Src/``, and/or ``amrex/Tools/Py_util/amrex_particles_to_vtp`` to generate a vtp file you can open with ParaView. @@ -397,6 +489,25 @@ Note also that amr.ref_ratio, amr.n_error_buf, amr.max_grid_size and amr.blocking_factor can be read in as a single value which is assigned to every level, or as multiple values, one for each level. +IAMR adjusts the tags near outflow faces after the user's tagging criteria have been +applied. The following must be preceded by "ns." + ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| | Description | Type | Default | ++======================+=======================================================================+=============+==============+ +| do_refine_outflow | If 1, and anything in the row of cells adjacent to an outflow face is | Int | 0 | +| | already tagged, refine that entire face. Avoids a coarse/fine | | | +| | boundary running along the outflow. | | | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| do_derefine_outflow | If 1, clear the tags within Nbuf_outflow cells of an outflow face, so | Int | 1 | +| | that fine grids never touch it. May not be set together with | | | +| | do_refine_outflow. | | | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ +| Nbuf_outflow | Number of level-0 cells to leave uncovered at an outflow face when | Int | 1 | +| | do_derefine_outflow = 1; rounded up to a multiple of the blocking | | | +| | factor. Must be > 0 in that case. | | | ++----------------------+-----------------------------------------------------------------------+-------------+--------------+ + As an example, consider: :: @@ -440,12 +551,12 @@ Tiling For details on IAMR's approach to tiling see :ref:`Chap:Parallel`. -The following inputs determine how we create the logical tiles and must be preceded by "fabarray_mfiter." : +The following inputs determine how we create the logical tiles and must be preceded by "fabarray." : +----------------------+-----------------------------------------------------------------------+----------+---------------+ | | Description | Type | Default | +======================+=======================================================================+==========+===============+ -| tile_size | Maximum number of cells in each direction for (logical) tiles. | IntVect | 1024000 | +| mfiter_tile_size | Maximum number of cells in each direction for (logical) tiles. | IntVect | 1024000 | | | (Default for 3D CPU-only) | | (1024000,8,8) | +----------------------+-----------------------------------------------------------------------+----------+---------------+ @@ -464,7 +575,7 @@ Here is some of the more frequently used options: +======================+=======================================================================+=============+==============+ | ns.v | Verbosity in IAMR routines | Int | 0 | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ -| ns.amr | Verbosity in AMR routines | Int | 0 | +| amr.v | Verbosity in AMR routines | Int | 0 | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ | particles.verbose | Verbosity in particle routines | Int | 0 | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ @@ -549,7 +660,7 @@ These control the MAC projection and must be preceded by "mac_proj.": Viscous and Diffusive Solve ~~~~~~~~~~~~~~~~~~~~~~~~~~~ -These control the diffusion solver and must be preceded by "diffusion.": +These control the diffusion solver and must be preceded by "diffuse.": +-------------------------+-----------------------------------------------------------------------+-------------+--------------+ | | Description | Type | Default | @@ -574,9 +685,9 @@ The following inputs must be preceded by "ns." +======================+=======================================================================+=============+==============+ | do_init_proj | Do the initial projections? False is primarily for debugging. | Bool | True | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ -| init_iter | How many pressure iterations before starting the first timestep. | Int | 3 | +| init_iter | How many pressure iterations before starting the first timestep. | Int | 2 | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ -| init_vel_iters | How many projection iterations to ensure the velocity satisfies the | Int | 3 | +| init_vel_iter | How many projection iterations to ensure the velocity satisfies the | Int | 1 | | | constraint. Set = 0 to skip this part of the initialization. | | | +----------------------+-----------------------------------------------------------------------+-------------+--------------+ diff --git a/Docs/sphinx_documentation/source/Software.rst b/Docs/sphinx_documentation/source/Software.rst index 431cea62d..c4fddbc2a 100644 --- a/Docs/sphinx_documentation/source/Software.rst +++ b/Docs/sphinx_documentation/source/Software.rst @@ -571,7 +571,7 @@ more than one grid per MPI rank, and different strategies for distributing the g IAMR relies on AMReX for the implementation. For more information, please see AMReX's documentation, found here: :ref:`amrex:Chap:ManagingGridHierarchy`. -See :ref:`sec:gridCreation` and :ref:`Chap:InputsLoadBalancing` for how grids are created, +See :ref:`sec:gridCreation` and :ref:`sec:InputsLoadBalancing` for how grids are created, i.e. how the :cpp:`BoxArray` on which :cpp:`MultiFabs` will be built is defined at each level. diff --git a/Docs/sphinx_documentation/source/TimeStep.rst b/Docs/sphinx_documentation/source/TimeStep.rst index 20cd89a83..ee911446a 100644 --- a/Docs/sphinx_documentation/source/TimeStep.rst +++ b/Docs/sphinx_documentation/source/TimeStep.rst @@ -2,57 +2,11 @@ To learn how the convective terms are constructed, see `AMReX-Hydro `_ -Time Step -- MOL -~~~~~~~~~~~~~~~~ - -In the predictor - -- Define :math:`U^{MAC,n}`, the face-centered (staggered) MAC velocity which is used for advection, using :math:`U^n` - -- Define an approximation to the new-time state, :math:`(\rho U)^{\ast}` by setting - -.. math:: (\rho U)^{\ast} &= (\rho U)^n - - \Delta t \left( \nabla \cdot (\rho U^{MAC} U) + \nabla {p}^{n-1/2} \right) \\ &+ - \Delta t \left( \nabla \cdot \tau^n + \sum_p \beta_p (V_p - {U}^{\ast}) + \rho g \right) - -- Project :math:`U^{\ast}` by solving - -.. math:: \nabla \cdot \frac{1}{\rho} \nabla \phi = \nabla \cdot \left( \frac{1}{\Delta t} - U^{\ast}+ \frac{1}{\rho} \nabla {p}^{n-1/2} \right) - -then defining - -.. math:: U^{\ast \ast} = U^{\ast} - \frac{\Delta t}{\rho} \nabla \phi - -and - -.. math:: {p}^{n+1/2, \ast} = \phi - - -In the corrector - -- Define :math:`U^{MAC,\ast \ast}` at the "new" time using :math:`U^{\ast \ast}` - -- Define a new approximation to the new-time state, :math:`(\rho U)^{\ast \ast \ast}` by setting - -.. math:: (\rho U)^{\ast \ast \ast} &= (\rho U)^n - \frac{\Delta t}{2} \left( \nabla \cdot (\rho U^{MAC} U)^n + \nabla \cdot (\rho U^{MAC} U)^{\ast \ast}\right) + \\ &+ \frac{\Delta t}{2} \left( \nabla \cdot \tau^n + \nabla \cdot \tau^{\ast \ast \ast} \right) + \Delta t \left( - \nabla {p}^{n+1/2,\ast} + \sum_p \beta_p (V_p - {U}^{\ast \ast \ast}) + \rho g \right) - -- Project :math:`U^{\ast \ast \ast}` by solving - -.. math:: \nabla \cdot \frac{1}{\rho} \nabla \phi = \nabla \cdot \left( \frac{1}{\Delta t} U^{\ast \ast \ast} + \frac{1}{\rho} \nabla {p}^{n+1/2,\ast} \right) - -then defining - -.. math:: U^{n+1} = U^{\ast \ast \ast} - \frac{\Delta t}{\rho} \nabla \phi - -and - -.. math:: {p}^{n+1/2} = \phi - Time Step -- Godunov ~~~~~~~~~~~~~~~~~~~~ -When we use the time-centered Godunov advection, we no longer need the predictor and corrector steps. +With the time-centered Godunov advection (``ns.advection_scheme = Godunov_PLM`` or +``Godunov_PPM``), a single step advances the solution from :math:`t^n` to :math:`t^{n+1}`. - Define the time-centered face-centered (staggered) MAC velocity which is used for advection: :math:`U^{MAC,n+1/2}` diff --git a/Docs/sphinx_documentation/source/Tutorials.rst b/Docs/sphinx_documentation/source/Tutorials.rst index 47bce8abe..738b3140a 100644 --- a/Docs/sphinx_documentation/source/Tutorials.rst +++ b/Docs/sphinx_documentation/source/Tutorials.rst @@ -109,11 +109,6 @@ Non-EB: We use :math:`p_0 = \mu = L = 1`. -* **Euler**: The test case is a "vortex tube" in a constant density fluid - in a triply periodic geometry. The refinement criteria are the - presence of a tracer and the magnitude of vorticity. - - * **TaylorGreen**: This case is an unsteady viscous benchmark for which the exact solution in 2D is @@ -122,7 +117,7 @@ Non-EB: v(x,y,t) &= -&& V_0 Cos(2\pi x) Sin(2\pi y) Cos(2\pi z) \exp(-2 (2\pi)^2 \nu t) \\ p(x,y,t) &= -&& \rho_0 V_0^2 \{Cos(4 \pi x) + Cos(4 \pi y)\} \exp(-4 (2\pi)^2 \nu t) / 4 - In ``TaylorGreen/benchmarks``, there is a tool, ViscBench2d.cpp, that reads a plot file and compares the solution against this exact solution. This benchmark was originally derived by G.I. Taylor (Phil. Mag., Vol. 46, No. 274, pp. 671-674, 1923) and Ethier & Steinman (Intl. J. Num. Meth. Fluids, Vol. 19, pp. 369-375, 1994) give the pressure field. + In ``TaylorGreen/benchmarks``, there is a tool, ViscBench.cpp, that reads a plot file and compares the solution against this exact solution. This benchmark was originally derived by G.I. Taylor (Phil. Mag., Vol. 46, No. 274, pp. 671-674, 1923) and Ethier & Steinman (Intl. J. Num. Meth. Fluids, Vol. 19, pp. 369-375, 1994) give the pressure field. In 3D, the problem is initialized with @@ -133,13 +128,6 @@ Non-EB: p(x,y,t) &= -&& \rho_0 V_0^2 \{2 + Cos(4 \pi z)\}\{Cos(4 \pi x) + Cos(4 \pi y)\} \exp(-4 (2\pi)^2 \nu t) / 16 -* **HIT**: Homogeneous isentropic forced turbulence with constant density. - This demonstrates defining a new forcing function by using a local edited - version of ``NS_getForce.cpp``. IAMR's make system is automatically configured - to select any local versions of files and ignore the corresponding versions in - ``IAMR/Source``. This problem is 3D only. - - * **Particles**: Particles in a double shear layer. Uses 2 levels of refinement and fixed grids. With fixed grids, a grid file (called ``fixed_grids_ml`` here) is used to define the grids for levels >= 1. diff --git a/Docs/sphinx_documentation/source/Visualization.rst b/Docs/sphinx_documentation/source/Visualization.rst index 4d3d76d69..ed4ed7eb5 100644 --- a/Docs/sphinx_documentation/source/Visualization.rst +++ b/Docs/sphinx_documentation/source/Visualization.rst @@ -696,7 +696,7 @@ pc.save(‘profile’) Density/velocity magnitude/kinetic energy phase plot -.. figure:: ./Visualization/Profile2D_1_Density_magveel_kineng.png +.. figure:: ./Visualization/Profile2D_1_Density_magvel_kineng.png :alt: Density/velocity magnitude/kinetic energy phase plot :width: 4.00000in diff --git a/Docs/sphinx_documentation/source/conf.py b/Docs/sphinx_documentation/source/conf.py index e22bba616..549e68d74 100644 --- a/Docs/sphinx_documentation/source/conf.py +++ b/Docs/sphinx_documentation/source/conf.py @@ -22,7 +22,6 @@ import re import sphinx_rtd_theme -import breathe from datetime import datetime def get_IAMR_version(): @@ -43,8 +42,7 @@ def get_IAMR_version(): 'sphinx.ext.githubpages', 'sphinx.ext.viewcode', 'sphinx.ext.intersphinx', - 'sphinx.ext.autosectionlabel', - 'breathe'] + 'sphinx.ext.autosectionlabel'] # bibtex bibtex_bibfiles = ["refs.bib"] @@ -61,7 +59,7 @@ def get_IAMR_version(): } # Add any paths that contain templates here, relative to this directory. -templates_path = ['ytemplates'] +templates_path = [] # The suffix(es) of source filenames. # You can specify multiple suffix as a list of string: @@ -106,22 +104,6 @@ def get_IAMR_version(): numfig = True -# -- breathe options ------------------------------------------------------ - -breathe_projects = { - "IAMR": "../../../out/docs_xml/doxygen/", - } - -breathe_default_project = "IAMR" - -breathe_default_members = ('members', 'undoc-members', 'protected-members', - 'private-members', 'content-only') - -breathe_doxygen_config_options = {'EXTRACT_ALL': 'YES', - 'SHOW_USED_FILES': 'YES', - 'RECURSIVE': 'YES'} - - # -- Options for HTML output ---------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for diff --git a/Exec/Make.IAMR b/Exec/Make.IAMR index 954bcf1b3..e1b8be224 100644 --- a/Exec/Make.IAMR +++ b/Exec/Make.IAMR @@ -61,13 +61,6 @@ ifeq ($(USE_VELOCITY), TRUE) include $(AMREX_HOME)/Src/Extern/amrdata/Make.package endif -ifeq ($(USE_TURBULENT_FORCING), TRUE) - DEFINES += -DAMREX_USE_TURBULENT_FORCING - ifeq ($(USE_FAST_FORCE), TRUE) - DEFINES += -DAMREX_USE_FAST_FORCE - endif -endif - # job_info support CEXE_sources += AMReX_buildInfo.cpp CEXE_headers += $(AMREX_HOME)/Tools/C_scripts/AMReX_buildInfo.H diff --git a/Exec/run3d/Make.package b/Exec/run3d/Make.package deleted file mode 100644 index 4169aa7af..000000000 --- a/Exec/run3d/Make.package +++ /dev/null @@ -1,6 +0,0 @@ -ifeq ($(USE_VELOCITY), TRUE) - CEXE_headers += DataServices.H AmrData.H XYPlotDataList.H AmrvisConstants.H - CEXE_sources += DataServices.cpp AmrData.cpp - FEXE_sources += FABUTIL_$(DIM)D.F -endif - diff --git a/Exec/square_grid_turbulence/Make.package b/Exec/square_grid_turbulence/Make.package deleted file mode 100644 index cdb8f2038..000000000 --- a/Exec/square_grid_turbulence/Make.package +++ /dev/null @@ -1,2 +0,0 @@ -FEXE_headers += probdata.H -F90EXE_sources += PROB_$(DIM)D.F90 diff --git a/Exec/square_grid_turbulence/inputs.3d.square_grid b/Exec/square_grid_turbulence/inputs.3d.square_grid index a9d5cd7eb..210083f65 100644 --- a/Exec/square_grid_turbulence/inputs.3d.square_grid +++ b/Exec/square_grid_turbulence/inputs.3d.square_grid @@ -25,11 +25,10 @@ amr.plot_int = 20 # Steps between plot files amr.plot_per = 0.1 # Steps between plot files amr.check_int = 100 # Steps between checkpoint files #amr.restart = chk13100 # Checkpoint to restart from -amr.probin_file = probin.3d.square_grid #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # PHYSICS # #.......................................# -ns.gravity = 0. 0. 0. # Gravitational force (3D) +ns.gravity = 0. # Gravitational force ns.vel_visc_coef = 1.157407407407407e-05 ns.scal_diff_coefs = 1.0 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# @@ -73,7 +72,7 @@ square_grid.dim_L0 = 0.1 square_grid.ratio_t0_L0_cross = 1.886792452830189e-01 square_grid.ratio_t0_stream_thickness = 0.25 -amr.derive_plot_vars = mag_vort avg_pressure velocity_average +amr.derive_plot_vars = mag_vort avg_pressure amr.plot_vars = x_velocity y_velocity z_velocity #******************************************************************************* diff --git a/Source/Diffusion.cpp b/Source/Diffusion.cpp index 13250590f..39279c73b 100644 --- a/Source/Diffusion.cpp +++ b/Source/Diffusion.cpp @@ -1523,12 +1523,12 @@ Diffusion::computeExtensiveFluxes(MLMG& a_mg, MultiFab& Soln, if ( flags.getType(amrex::grow(bx,0)) == FabType::covered ) { // - // For now, set to very large num so we know if you accidentally use it - // MLMG will set covered fluxes to zero + // Covered fluxes must be zero: they are passed to the viscous + // FluxRegister, which does not know about EB. // - AMREX_D_TERM(AMREX_PARALLEL_FOR_4D(ubx, ncomp, i, j, k, n, {fx(i,j,k,n) = COVERED_VAL;});, - AMREX_PARALLEL_FOR_4D(vbx, ncomp, i, j, k, n, {fy(i,j,k,n) = COVERED_VAL;});, - AMREX_PARALLEL_FOR_4D(wbx, ncomp, i, j, k, n, {fz(i,j,k,n) = COVERED_VAL;});); + AMREX_D_TERM(AMREX_PARALLEL_FOR_4D(ubx, ncomp, i, j, k, n, {fx(i,j,k,n) = 0.0;});, + AMREX_PARALLEL_FOR_4D(vbx, ncomp, i, j, k, n, {fy(i,j,k,n) = 0.0;});, + AMREX_PARALLEL_FOR_4D(wbx, ncomp, i, j, k, n, {fz(i,j,k,n) = 0.0;});); } else if ( flags.getType(amrex::grow(bx,0)) != FabType::regular ) { diff --git a/Source/NS_bcfill.H b/Source/NS_bcfill.H index 42e77b450..b45191173 100644 --- a/Source/NS_bcfill.H +++ b/Source/NS_bcfill.H @@ -319,36 +319,35 @@ void dummy_fill (Box const& bx, FArrayBox& data, } -// struct NodalFillExtDir -// { -// AMREX_GPU_DEVICE -// void operator()( -// const amrex::IntVect& iv, -// amrex::Array4 const& dest, -// const int dcomp, -// const int numcomp, -// amrex::GeometryData const& geom, -// const amrex::Real time, -// const amrex::BCRec* bcr, -// const int bcomp, -// const int orig_comp) const -// { -// // do something for external Dirichlet (BCType::ext_dir) -// } -// }; +// +// Nodal pressure: press_bc never contains ext_dir and GpuBndryFuncFab already +// extrapolates nodal ghost values, so there is nothing left to do here. +// +struct NodalFillExtDir +{ + AMREX_GPU_DEVICE + void operator()( + const amrex::IntVect& /*iv*/, + amrex::Array4 const& /*dest*/, + const int /*dcomp*/, + const int /*numcomp*/, + amrex::GeometryData const& /*geom*/, + const amrex::Real /*time*/, + const amrex::BCRec* /*bcr*/, + const int /*bcomp*/, + const int /*orig_comp*/) const + {} +}; inline -void press_fill (Box const& /*bx*/, FArrayBox& /*data*/, - const int /*dcomp*/, const int /*numcomp*/, - Geometry const& /*geom*/, const Real /*time*/, - const Vector& /*bcr*/, const int /*bcomp*/, - const int /*scomp*/) +void press_fill (Box const& bx, FArrayBox& data, + const int dcomp, const int numcomp, + Geometry const& geom, const Real time, + const Vector& bcr, const int bcomp, + const int scomp) { - amrex::Abort("press_fill: Need to write fill for external Dirichlet (BCType::ext_dir)"); - - // GpuBndryFuncFab gpu_bndry_func(NodalFillExtDir{}); - // gpu_bndry_func(bx,data,dcomp,numcomp,geom,time,bcr,bcomp,scomp); - + GpuBndryFuncFab gpu_bndry_func(NodalFillExtDir{}); + gpu_bndry_func(bx,data,dcomp,numcomp,geom,time,bcr,bcomp,scomp); } #endif diff --git a/Source/NavierStokes.cpp b/Source/NavierStokes.cpp index 2edc390e3..643074fef 100644 --- a/Source/NavierStokes.cpp +++ b/Source/NavierStokes.cpp @@ -1791,20 +1791,12 @@ NavierStokes::reflux () } } -#ifdef AMREX_USE_EB - fr_adv.Reflux(Vsync,*volfrac, 0, 0, AMREX_SPACEDIM); - fr_adv.Reflux(Ssync,*volfrac, AMREX_SPACEDIM, 0, NUM_STATE-AMREX_SPACEDIM); -#else - fr_adv.Reflux(Vsync, 0, 0, AMREX_SPACEDIM); - fr_adv.Reflux(Ssync, AMREX_SPACEDIM, 0,NUM_STATE-AMREX_SPACEDIM); -#endif - const Real scale = 1.0/dt_crse; - Vsync.mult(scale); - Ssync.mult(scale); - const BoxArray& fine_boxes = getLevel(level+1).boxArray(); // // Zero out coarse grid cells which underlie fine grid cells. + // Done before the advective reflux: with EB it re-redistributes part + // of a cut cell's correction into covered cells, which SyncInterp + // must carry to the fine level. // BoxArray baf = fine_boxes; @@ -1839,6 +1831,17 @@ NavierStokes::reflux () } } } + +#ifdef AMREX_USE_EB + fr_adv.Reflux(Vsync,*volfrac, 0, 0, AMREX_SPACEDIM); + fr_adv.Reflux(Ssync,*volfrac, AMREX_SPACEDIM, 0, NUM_STATE-AMREX_SPACEDIM); +#else + fr_adv.Reflux(Vsync, 0, 0, AMREX_SPACEDIM); + fr_adv.Reflux(Ssync, AMREX_SPACEDIM, 0,NUM_STATE-AMREX_SPACEDIM); +#endif + const Real scale = 1.0/dt_crse; + Vsync.mult(scale); + Ssync.mult(scale); } // diff --git a/Source/NavierStokesBase.cpp b/Source/NavierStokesBase.cpp index f1920d7bf..6498f5425 100644 --- a/Source/NavierStokesBase.cpp +++ b/Source/NavierStokesBase.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -23,10 +24,6 @@ #include #endif -#ifdef AMREX_USE_TURBULENT_FORCING -#include -#endif - using namespace amrex; @@ -995,7 +992,19 @@ NavierStokesBase::computeNewDt (int finest_level, for (i = 0; i <= finest_level; i++) { NavierStokesBase& adv_level = getLevel(i); - dt_min[i] = std::min(dt_min[i],adv_level.estTimeStep()); + // + // dt_min[i] arrives holding the CFL-based estimate returned by + // advance(). With fixed_dt that estimate must not undercut fixed_dt; + // only ramp by change_max from an init_shrink-reduced start. Ramp from + // no less than init_shrink*fixed_dt, so a restart after a step cut + // short by stop_time does not crawl back up. + // + if (fixed_dt > 0.0) { + const Real est = adv_level.estTimeStep(); + dt_min[i] = std::min(change_max*std::max(dt_level[i],init_shrink*est),est); + } else { + dt_min[i] = std::min(dt_min[i],adv_level.estTimeStep()); + } } if (fixed_dt <= 0.0) @@ -1252,7 +1261,7 @@ NavierStokesBase::create_umac_grown (int nGrow, for(int jj(-1); jj<=1; jj++) { for(int ii(-1); ii<=1; ii++) { if ( Math::abs(ii)+Math::abs(jj)+Math::abs(kk) == 1 && - (maskarr(i+ii,j+jj,k+kk) == interior || maskarr(i+ii,j+jj,k+kk) == level_mask_covered) ) + (maskarr(i+ii,j+jj,k+kk) == level_mask_interior || maskarr(i+ii,j+jj,k+kk) == level_mask_covered) ) { count++; } @@ -2587,14 +2596,6 @@ NavierStokesBase::post_restart () } } -#ifdef AMREX_USE_TURBULENT_FORCING - // - // Initialize data structures used for homogeneous isentropic forced turbulence. - // Only need to do it once. - if (level == 0) - TurbulentForcing::init_turbulent_forcing(geom.ProbLoArray(),geom.ProbHiArray()); -#endif - #ifdef AMREX_PARTICLES post_restart_particle (); #endif @@ -2671,7 +2672,11 @@ NavierStokesBase::post_timestep (int crse_iteration) BoxArray ba(bx); DistributionMapping dm{ba}; - MultiFab mf(ba, dm, AMREX_SPACEDIM, 0, MFInfo(), Factory()); + // + // ba is not this level's BoxArray, so this level's (EB) factory + // must not be used to build the slab. + // + MultiFab mf(ba, dm, AMREX_SPACEDIM, 0); mf.ParallelCopy(get_new_data(State_Type), Xvel, 0, AMREX_SPACEDIM); @@ -4063,7 +4068,18 @@ NavierStokesBase::ParticleDerive (const std::string& name, } MultiFab* ret = new MultiFab(grids, dmap, ncomp, ngrow, MFInfo(), Factory()); + // + // The MultiFab overload only writes the valid region. Callers such as + // errorEst read the ghost cells (AmrLevel::derive fills them), so + // fill them from neighbouring grids and extrapolate at physical and + // coarse/fine boundaries. + // + ret->setVal(0.); ParticleDerive(name,time,*ret,0); + if (ngrow > 0) { + ret->FillBoundary(geom.periodicity()); + Extrapolater::FirstOrderExtrap(*ret, geom, 0, ncomp, ngrow); + } return std::unique_ptr{ret}; } else { @@ -5100,8 +5116,14 @@ NavierStokesBase::ComputeAofs ( MultiFab& advc, int a_comp, // Advection term "A int as_crse = (fr_as_crse != nullptr); int as_fine = (fr_as_fine != nullptr); - FArrayBox* p_drho_as_crse = (fr_as_crse) ? - fr_as_crse->getCrseData(mfi) : &fab_drho_as_crse; + // + // The register's coarse data holds all NUM_STATE components, + // but the redistribution writes components 0..ncomp-1 of the + // array it is handed, so offset to this call's state range. + // + Array4 const drho_as_crse_arr = (fr_as_crse) ? + fr_as_crse->getCrseData(mfi)->array(state_indx) + : fab_drho_as_crse.array(); const IArrayBox* p_rrflag_as_crse = (fr_as_crse) ? fr_as_crse->getCrseFlag(mfi) : &fab_rrflag_as_crse; @@ -5117,7 +5139,7 @@ NavierStokesBase::ComputeAofs ( MultiFab& advc, int a_comp, // Advection term "A AMREX_D_DECL(apx,apy,apz), vfrac_arr, AMREX_D_DECL(fcx,fcy,fcz), ccent_arr, bcrec_d, geom, dt, redistribution_type, - as_crse, p_drho_as_crse->array(), p_rrflag_as_crse->array(), + as_crse, drho_as_crse_arr, p_rrflag_as_crse->array(), as_fine, dm_as_fine.array(), coarse_fine_mask->const_array(mfi), level_mask_notcovered, /*fac_for_deltaR*/ sync_factor, diff --git a/Source/Projection.cpp b/Source/Projection.cpp index 78f6d1f6b..bb5e43906 100644 --- a/Source/Projection.cpp +++ b/Source/Projection.cpp @@ -1287,9 +1287,14 @@ Projection::scaleVar (MultiFab* sig, } }); + // + // Scale vel on its own ghost cells (the nodal divu stencil reads + // them), not only on sig's, which may have none. + // + const Box& vbx = mfi.growntilebox(1); auto const& velarr = vel->array(mfi); - amrex::ParallelFor(bx, AMREX_SPACEDIM, [=] + amrex::ParallelFor(vbx, AMREX_SPACEDIM, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { // NOTE: cells outside the domain in the axial (j) direction are @@ -1410,9 +1415,14 @@ Projection::rescaleVar (MultiFab* sig, } }); + // + // Scale vel on its own ghost cells (the nodal divu stencil reads + // them), not only on sig's, which may have none. + // + const Box& vbx = mfi.growntilebox(1); auto const& velarr = vel->array(mfi); - amrex::ParallelFor(bx, AMREX_SPACEDIM, [=] + amrex::ParallelFor(vbx, AMREX_SPACEDIM, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { // Mirrors scaleVar: the axial out-of-domain ghosts were scaled diff --git a/Tutorials/HIT/.gitignore b/Tutorials/HIT/.gitignore deleted file mode 100644 index b7beb8714..000000000 --- a/Tutorials/HIT/.gitignore +++ /dev/null @@ -1 +0,0 @@ -*.0.0 diff --git a/Tutorials/HIT/GNUmakefile b/Tutorials/HIT/GNUmakefile deleted file mode 100644 index d0afcd2d0..000000000 --- a/Tutorials/HIT/GNUmakefile +++ /dev/null @@ -1,56 +0,0 @@ -#AMREX_HOME defines the directory in which we will find the AMReX directory -AMREX_HOME ?= ../../../amrex -AMREX_HYDRO_HOME ?= ../../../AMReX-Hydro - -#TOP defines the directory in which we will find Source, Exec, etc. -TOP = ../.. - -# -# Variables for the user to set ... -# - -DIM = 3 -COMP = gcc -DEBUG = FALSE -USE_MPI = TRUE -USE_OMP = FALSE - -USE_CUDA = FALSE - -USE_TURBULENT_FORCING = TRUE -USE_FAST_FORCE = FALSE - -#TEST=TRUE -#USE_ASSERTION=TRUE - -PRECISION = DOUBLE - -USE_HYPRE = FALSE -USE_METIS = FALSE - -USE_VELOCITY = TRUE -USE_VELOCITY = FALSE - -USE_XBLAS = TRUE -USE_XBLAS = FALSE - -USE_SENSEI_INSITU = FALSE - -EBASE = amr - -ifeq (${USE_XBLAS}, TRUE) - XTRADEFS += -DXBLAS - XTRAINCLOC += $(HOME)/tmp/xblas-1.0.248/src - XTRALIBLOC += $(HOME)/tmp/xblas-1.0.248/src - XTRALIBS += -lxblas -endif - -ifeq (${USE_VELOCITY}, TRUE) - #AMRVIS_DIR defines the directory in which we will find pAmrvis (e.g. DataServices, AmrData and FABUTIL) - AMRVIS_DIR = $(AMREX_HOME)/Src/Extern/amrdata -endif - -Blocs := . - -include ./Make.package -include $(TOP)/Exec/Make.IAMR diff --git a/Tutorials/HIT/Make.package b/Tutorials/HIT/Make.package deleted file mode 100644 index 026272b01..000000000 --- a/Tutorials/HIT/Make.package +++ /dev/null @@ -1,10 +0,0 @@ -ifeq ($(USE_VELOCITY), TRUE) - CEXE_headers += DataServices.H AmrData.H XYPlotDataList.H AmrvisConstants.H - CEXE_sources += DataServices.cpp AmrData.cpp - FEXE_sources += FABUTIL_$(DIM)D.F -endif - -ifeq ($(USE_TURBULENT_FORCING), TRUE) - CEXE_headers += TurbulentForcing_params.H TurbulentForcing_def.H depRand.H - CEXE_sources += depRand.cpp -endif diff --git a/Tutorials/HIT/NS_getForce.cpp b/Tutorials/HIT/NS_getForce.cpp deleted file mode 100644 index b25379a76..000000000 --- a/Tutorials/HIT/NS_getForce.cpp +++ /dev/null @@ -1,780 +0,0 @@ - -#include -#ifdef AMREX_USE_TURBULENT_FORCING -#include -#include -#endif - - -using namespace amrex; - -// -// Virtual access function for getting the forcing terms for the -// velocities and scalars. The base version computes a buoyancy. -// -// NOTE: This function returns a rho weighted source term. -// -// For conservative (i.e. do_mom_diff=1, do_cons_trac=1), velocities -// are integrated according to the equation -// -// ui_t + uj ui_j = S_ui ===> tforces = rho S_ui -// -// and scalars psi where (psi = rho q = Scal) as -// -// psi_t + (uj psi)_j = S_psi ===> tforces = S_psi = rho S_q -// -// For non-conservative, this rho-weighted source term will get divided -// by rho in the predict_velocity, velocity_advection, scalar_advection, -// and advection_update routines. -// -// For temperature (which is always non-conservative), we evolve -// -// dT/dt - U dot grad T = [del dot lambda grad T + S_T] / (rho*c_p) -// ===> tforces = S_T/c_p -// -// -// For user-defined forcing, this means -// - For conservative variables, the force term computed here gets used -// as-is -// - For non-conservative variables, the force term computed here is -// divided by rho before use -// - -void -NavierStokesBase::getForce (FArrayBox& force, - const Box& bx, - int scomp, - int ncomp, - const Real time, - const FArrayBox& Vel, - const FArrayBox& Scal, - int scalScomp, - const MFIter& mfi) -{ - - const Real* VelDataPtr = Vel.dataPtr(); - const Real* ScalDataPtr = Scal.dataPtr(scalScomp); - - const Real grav = gravity; - const int* f_lo = force.loVect(); - const int* f_hi = force.hiVect(); - const int* v_lo = Vel.loVect(); - const int* v_hi = Vel.hiVect(); - const int* s_lo = Scal.loVect(); - const int* s_hi = Scal.hiVect(); - - if (ParallelDescriptor::IOProcessor() && getForceVerbose) { - amrex::Print() << "NavierStokesBase::getForce(): Entered..." << '\n' - << "time = " << time << '\n' - << "scomp = " << scomp << '\n' - << "ncomp = " << ncomp << '\n' - << "scalScomp = " << scalScomp << '\n'; - - if (scomp==0) - if (ncomp==3) amrex::Print() << "Doing velocities only" << '\n'; - else amrex::Print() << "Doing all components" << '\n'; - else if (scomp==3) - if (ncomp==1) amrex::Print() << "Doing density only" << '\n'; - else amrex::Print() << "Doing all scalars" << '\n'; - else if (scomp==4) amrex::Print() << "Doing tracer only" << '\n'; - else amrex::Print() << "Doing individual scalar" << '\n'; - - amrex::Print() << "NavierStokesBase::getForce(): Filling Force on box:" - << bx << '\n'; -#if (AMREX_SPACEDIM == 3) - amrex::Print() << "NavierStokesBase::getForce(): Force Domain:" << '\n'; - amrex::Print() << "(" << f_lo[0] << "," << f_lo[1] << "," << f_lo[2] << ") - " - << "(" << f_hi[0] << "," << f_hi[1] << "," << f_hi[2] << ")" << '\n'; - amrex::Print() << "NavierStokesBase::getForce(): Vel Domain:" << '\n'; - amrex::Print() << "(" << v_lo[0] << "," << v_lo[1] << "," << v_lo[2] << ") - " - << "(" << v_hi[0] << "," << v_hi[1] << "," << v_hi[2] << ")" << '\n'; - amrex::Print() << "NavierStokesBase::getForce(): Scal Domain:" << '\n'; - amrex::Print() << "(" << s_lo[0] << "," << s_lo[1] << "," << s_lo[2] << ") - " - << "(" << s_hi[0] << "," << s_hi[1] << "," << s_hi[2] << ")" << '\n'; -#else - amrex::Print() << "NavierStokesBase::getForce(): Force Domain:" << '\n'; - amrex::Print() << "(" << f_lo[0] << "," << f_lo[1] << ") - " - << "(" << f_hi[0] << "," << f_hi[1] << ")" << '\n'; - amrex::Print() << "NavierStokesBase::getForce(): Vel Domain:" << '\n'; - amrex::Print() << "(" << v_lo[0] << "," << v_lo[1] << ") - " - << "(" << v_hi[0] << "," << v_hi[1] << ")" << '\n'; - amrex::Print() << "NavierStokesBase::getForce(): Scal Domain:" << '\n'; - amrex::Print() << "(" << s_lo[0] << "," << s_lo[1] << ") - " - << "(" << s_hi[0] << "," << s_hi[1] << ")" << '\n'; -#endif - - Vector velmin(AMREX_SPACEDIM), velmax(AMREX_SPACEDIM); - Vector scalmin(NUM_SCALARS), scalmax(NUM_SCALARS); - for (int n=0; nvelmax[n]) velmax[n] = v; - } - } - } -#if (AMREX_SPACEDIM == 3) - } -#endif - for (int n=0; nscalmax[n]) scalmax[n] = s; - } - } - } -#if (AMREX_SPACEDIM == 3) - } -#endif - for (int n=0; n=AMREX_SPACEDIM); - } - - if ( scomp==Xvel ){ - // - // TODO: add some switch for user-supplied/problem-dependent forcing - // - auto const& frc = force.array(scomp); - auto const& scal = Scal.array(scalScomp); - - if ( std::abs(grav) > 0.0001) { - amrex::ParallelFor(bx, [frc, scal, grav] - AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - frc(i,j,k,0) = Real(0.0); -#if ( AMREX_SPACEDIM == 2 ) - frc(i,j,k,1) = grav*scal(i,j,k,0); -#elif ( AMREX_SPACEDIM == 3 ) - frc(i,j,k,1) = Real(0.0); - frc(i,j,k,2) = grav*scal(i,j,k,0); -#endif - }); - } - else { - force.setVal(0.0, bx, Xvel, AMREX_SPACEDIM); - } - -#ifdef AMREX_USE_TURBULENT_FORCING - // - // Homogeneous Isotropic Forced Turbulence - // - - // Physical coordinates of the lower left corner of the domain - auto const& problo = geom.ProbLoArray(); - // Physical coordinates of the upper right corner of the domain - auto const& probhi = geom.ProbHiArray(); - - Real Lx = probhi[0]-problo[0]; - Real Ly = probhi[1]-problo[1]; - Real Lz = probhi[2]-problo[2]; - Real Lmin = min(Lx,Ly,Lz); - - - // For now, only works in 3D and without tiling - AMREX_ALWAYS_ASSERT(bx==force.box()); - AMREX_ASSERT(AMREX_SPACEDIM==3); - - int xstep = static_cast(Lx/Lmin+0.5); - int ystep = static_cast(Ly/Lmin+0.5); - int zstep = static_cast(Lz/Lmin+0.5); - - Real kappaMax = TurbulentForcing::nmodes/Lmin + 1.0e-8; - - // force array bounds - // FIXME -- think about how bx (the box we want to fill) may not be the same as the force box! - int ilo = f_lo[0]; - int jlo = f_lo[1]; - int klo = f_lo[2]; - - auto const& dx = geom.CellSizeArray(); - Real hx = dx[0]; - Real hy = dx[1]; - Real hz = dx[2]; - - // Separate out forcing data into individual Array4's - int i_arr = 0; - int fd_ncomp = 1; - int num_elmts=TurbulentForcing::array_size*TurbulentForcing::array_size*TurbulentForcing::array_size; - Dim3 fd_begin{0,0,0}; - Dim3 fd_end{TurbulentForcing::array_size,TurbulentForcing::array_size,TurbulentForcing::array_size}; - - Array4 FTX(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 TAT(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPX(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPY(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPZ(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FAX(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FAY(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FAZ(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPXX(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPXY(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPXZ(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPYX(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPYY(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPYZ(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPZX(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPZY(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - i_arr++; - Array4 FPZZ(&TurbulentForcing::forcedata[i_arr*num_elmts], fd_begin, fd_end, fd_ncomp); - -//fixme -// check the arrays - // { - // int i = 0; - // int j = 0; - // int k = 0; - - // Print()<<"FTX, TAT: "< xlo = {AMREX_D_DECL(loc_lo[0],loc_lo[1],loc_lo[2])}; - - -#ifdef AMREX_USE_FAST_FORCE - // - // Construct force at fewer points and then interpolate. - // This is much faster on CPU. - // - - int ihi = f_hi[0]; - int jhi = f_hi[1]; - int khi = f_hi[2]; - - // coarse cell size - Real ff_hx = hx*TurbulentForcing::ff_factor; - Real ff_hy = hy*TurbulentForcing::ff_factor; - Real ff_hz = hz*TurbulentForcing::ff_factor; - - // coarse bounds (without accounting for ghost cells) - // FIXME -- think about how bx (the box we want to fill) may not be the same as the force box! - int ff_ilo = ilo/TurbulentForcing::ff_factor; - int ff_jlo = jlo/TurbulentForcing::ff_factor; - int ff_klo = klo/TurbulentForcing::ff_factor; - - // +1 so that the (ff_i+1) node used by the interpolation below always exists - int ff_ihi = ihi/TurbulentForcing::ff_factor + 1; - int ff_jhi = jhi/TurbulentForcing::ff_factor + 1; - int ff_khi = khi/TurbulentForcing::ff_factor + 1; - - // adjust for ghost cells - if (ilo < (ff_ilo*TurbulentForcing::ff_factor)) { - ff_ilo=ff_ilo-1; - } - if (jlo < (ff_jlo*TurbulentForcing::ff_factor)) { - ff_jlo=ff_jlo-1; - } - if (klo < (ff_klo*TurbulentForcing::ff_factor)) { - ff_klo=ff_klo-1; - } - - // allocate coarse force array - Box ffbx(IntVect(ff_ilo, ff_jlo, ff_klo), IntVect(ff_ihi, ff_jhi, ff_khi)); - FArrayBox ff_force(ffbx,AMREX_SPACEDIM); - // The Elixir defers freeing ff_force's device memory until the - // asynchronous kernels below that read it have completed. No-op on CPU. - Elixir ff_force_i = ff_force.elixir(); - const auto& ffarr = ff_force.array(); - - // Construct node-based coarse forcing - amrex::ParallelFor(ffbx, [ = ] - AMREX_GPU_DEVICE (int i, int j, int k ) noexcept - { - Real z = xlo[2] + ff_hz*(k-ff_klo); - Real y = xlo[1] + ff_hy*(j-ff_jlo); - Real x = xlo[0] + ff_hx*(i-ff_ilo); - - for (int n = 0; n < AMREX_SPACEDIM; n++) - ffarr(i,j,k,n) = 0.0; - - // forcedata (and Array4) has column-major layout - for (int kz = TurbulentForcing::mode_start*zstep; kz <= TurbulentForcing::nmodes*zstep; kz += zstep) { - for (int ky = TurbulentForcing::mode_start*ystep; ky <= TurbulentForcing::nmodes*ystep; ky += ystep) { - for (int kx = TurbulentForcing::mode_start*xstep; kx <= TurbulentForcing::nmodes*xstep; kx += xstep) - { - Real kappa = sqrt( (kx*kx)/(Lx*Lx) + (ky*ky)/(Ly*Ly) + (kz*kz)/(Lz*Lz) ); - - if (kappa <= kappaMax) - { - Real xT = cos(FTX(kx,ky,kz)*time + TAT(kx,ky,kz)); - - if ( TurbulentForcing::div_free_force ) - { - ffarr(i,j,k,0) += xT * - ( FAZ(kx,ky,kz)*TwoPi*(ky/Ly) - * sin(TwoPi*kx*x/Lx+FPZX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPZY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZZ(kx,ky,kz)) - - FAY(kx,ky,kz)*TwoPi*(kz/Lz) - * sin(TwoPi*kx*x/Lx+FPYX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPYY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPYZ(kx,ky,kz)) ); - - ffarr(i,j,k,1) += xT * - ( FAX(kx,ky,kz)*TwoPi*(kz/Lz) - * sin(TwoPi*kx*x/Lx+FPXX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPXY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPXZ(kx,ky,kz)) - - FAZ(kx,ky,kz)*TwoPi*(kx/Lx) - * cos(TwoPi*kx*x/Lx+FPZX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPZY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZZ(kx,ky,kz)) ); - - ffarr(i,j,k,2) += xT * - ( FAY(kx,ky,kz)*TwoPi*(kx/Lx) - * cos(TwoPi*kx*x/Lx+FPYX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPYY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPYZ(kx,ky,kz)) - - FAX(kx,ky,kz)*TwoPi*(ky/Ly) - * sin(TwoPi*kx*x/Lx+FPXX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPXY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPXZ(kx,ky,kz)) ); - } - else - { - - ffarr(i,j,k,0) += xT*FAX(kx,ky,kz)*cos(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - - ffarr(i,j,k,1) += xT*FAY(kx,ky,kz)*sin(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - - ffarr(i,j,k,2) += xT*FAZ(kx,ky,kz)*sin(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - } - } - } - } - } - - // - // For high aspect ratio domain, add more modes to break symmetry at a low level. - // We assume Lz is longer, Lx = Ly. - // - for ( int kz = 1; kz <= zstep-1; kz++) { - for ( int ky = TurbulentForcing::mode_start; ky <= TurbulentForcing::nmodes*ystep; ky++) { - for ( int kx = TurbulentForcing::mode_start; kx <= TurbulentForcing::nmodes*xstep; kx++) - { - Real kappa = sqrt( (kx*kx)/(Lx*Lx) + (ky*ky)/(Ly*Ly) + (kz*kz)/(Lz*Lz) ); - - if (kappa <= kappaMax) - { - Real xT = cos(FTX(kx,ky,kz)*time + TAT(kx,ky,kz)); - - if ( TurbulentForcing::div_free_force ) - { - ffarr(i,j,k,0) += xT * - ( FAZ(kx,ky,kz)*TwoPi*(ky/Ly) - * sin(TwoPi*kx*x/Lx+FPZX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPZY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZZ(kx,ky,kz)) - - FAY(kx,ky,kz)*TwoPi*(kz/Lz) - * sin(TwoPi*kx*x/Lx+FPYX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPYY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPYZ(kx,ky,kz)) ); - - ffarr(i,j,k,1) += xT * - ( FAX(kx,ky,kz)*TwoPi*(kz/Lz) - * sin(TwoPi*kx*x/Lx+FPXX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPXY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPXZ(kx,ky,kz)) - - FAZ(kx,ky,kz)*TwoPi*(kx/Lx) - * cos(TwoPi*kx*x/Lx+FPZX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPZY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZZ(kx,ky,kz)) ); - - ffarr(i,j,k,2) += xT * - ( FAY(kx,ky,kz)*TwoPi*(kx/Lx) - * cos(TwoPi*kx*x/Lx+FPYX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPYY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPYZ(kx,ky,kz)) - - FAX(kx,ky,kz)*TwoPi*(ky/Ly) - * sin(TwoPi*kx*x/Lx+FPXX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPXY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPXZ(kx,ky,kz)) ); - } - else - { - - ffarr(i,j,k,0) += xT*FAX(kx,ky,kz)*cos(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - - ffarr(i,j,k,1) += xT*FAY(kx,ky,kz)*sin(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - - ffarr(i,j,k,2) += xT*FAZ(kx,ky,kz)*sin(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - } - } - } - } - } - - }); - - // Need all of ffarr filled for next lambda - amrex::Gpu::synchronize(); - - // Now interpolate onto fine grid - // bx is the box we want to fill and may be smaller than force - auto const& dens = Scal.array(scalScomp); - - amrex::ParallelFor(bx, AMREX_SPACEDIM, [ = ] - AMREX_GPU_DEVICE (int i, int j, int k, int n ) noexcept - { - int ff_k = k/TurbulentForcing::ff_factor; - int ff_j = j/TurbulentForcing::ff_factor; - int ff_i = i/TurbulentForcing::ff_factor; - - Real zd = ( hz*(k-klo + 0.5) - ff_hz*(ff_k-ff_klo) )/ff_hz; - Real yd = ( hy*(j-jlo + 0.5) - ff_hy*(ff_j-ff_jlo) )/ff_hy; - Real xd = ( hx*(i-ilo + 0.5) - ff_hx*(ff_i-ff_ilo) )/ff_hx; - - Real ff00 = ffarr(ff_i ,ff_j ,ff_k ,n) * (1. - xd) - + ffarr(ff_i+1,ff_j ,ff_k ,n) * xd; - Real ff01 = ffarr(ff_i ,ff_j ,ff_k+1,n) * (1. - xd) - + ffarr(ff_i+1,ff_j ,ff_k+1,n) * xd; - Real ff10 = ffarr(ff_i ,ff_j+1,ff_k ,n) * (1. - xd) - + ffarr(ff_i+1,ff_j+1,ff_k ,n) * xd; - Real ff11 = ffarr(ff_i ,ff_j+1,ff_k+1,n) * (1. - xd) - + ffarr(ff_i+1,ff_j+1,ff_k+1,n) * xd; - - Real ff = ( ff00*(1.-yd)+ff10*yd ) * (1. - zd) - + ( ff01*(1.-yd)+ff11*yd ) * zd; - - frc(i,j,k,n) += dens(i,j,k,0) * ff; - - // if ( i==0 && j==0 && k==0 ){ - // printf("%15.13e %15.13e %15.13e %15.13e %15.13e\n", - // FTX(1,0,0), TAT(0,0,0),FAY(9,8,7), FPX(3,4,5), FPZZ(12,4,7) ); - // } - }); - -#else - - // - // Original implementation using all 33 points in k-space. - // May be fast enough on GPU. - // - - auto const& dens = Scal.array(scalScomp); - - // Construct cell-centered forcing - amrex::ParallelFor(bx, [ = ] - AMREX_GPU_DEVICE (int i, int j, int k ) noexcept - { - Real z = xlo[2] + hz*(k-klo + 0.5); - Real y = xlo[1] + hy*(j-jlo + 0.5); - Real x = xlo[0] + hx*(i-ilo + 0.5); - - Real f1 = 0; - Real f2 = 0; - Real f3 = 0; - - // forcedata (and Array4) has column-major layout - for (int kz = TurbulentForcing::mode_start*zstep; kz <= TurbulentForcing::nmodes*zstep; kz += zstep) { - for (int ky = TurbulentForcing::mode_start*ystep; ky <= TurbulentForcing::nmodes*ystep; ky += ystep) { - for (int kx = TurbulentForcing::mode_start*xstep; kx <= TurbulentForcing::nmodes*xstep; kx += xstep) - { - Real kappa = sqrt( (kx*kx)/(Lx*Lx) + (ky*ky)/(Ly*Ly) + (kz*kz)/(Lz*Lz) ); - - if (kappa <= kappaMax) - { - Real xT = cos(FTX(kx,ky,kz)*time + TAT(kx,ky,kz)); - - // if ( i==0 && j==0 && k==0 && kx==0 && ky==0 && kx==0){ - // printf("(0,0,0) : xT : %15.13e %15.13e %15.13e %15.13e\n", - // FTX(kx,ky,kz), time, TAT(kx,ky,kz), xT); - // Abort(); - // } - - - if ( TurbulentForcing::div_free_force ) - { - f1 += xT * - ( FAZ(kx,ky,kz)*TwoPi*(ky/Ly) - * sin(TwoPi*kx*x/Lx+FPZX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPZY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZZ(kx,ky,kz)) - - FAY(kx,ky,kz)*TwoPi*(kz/Lz) - * sin(TwoPi*kx*x/Lx+FPYX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPYY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPYZ(kx,ky,kz)) ); - - f2 += xT * - ( FAX(kx,ky,kz)*TwoPi*(kz/Lz) - * sin(TwoPi*kx*x/Lx+FPXX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPXY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPXZ(kx,ky,kz)) - - FAZ(kx,ky,kz)*TwoPi*(kx/Lx) - * cos(TwoPi*kx*x/Lx+FPZX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPZY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZZ(kx,ky,kz)) ); - - f3 += xT * - ( FAY(kx,ky,kz)*TwoPi*(kx/Lx) - * cos(TwoPi*kx*x/Lx+FPYX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPYY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPYZ(kx,ky,kz)) - - FAX(kx,ky,kz)*TwoPi*(ky/Ly) - * sin(TwoPi*kx*x/Lx+FPXX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPXY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPXZ(kx,ky,kz)) ); - } - else - { - - f1 += xT*FAX(kx,ky,kz)*cos(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - - f2 += xT*FAY(kx,ky,kz)*sin(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - - f3 += xT*FAZ(kx,ky,kz)*sin(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - } - } - } - } - } - - // - // For high aspect ratio domain, add more modes to break symmetry at a low level. - // We assume Lz is longer, Lx = Ly. - // - for ( int kz = 1; kz <= zstep-1; kz++) { - for ( int ky = TurbulentForcing::mode_start; ky <= TurbulentForcing::nmodes*ystep; ky++) { - for ( int kx = TurbulentForcing::mode_start; kx <= TurbulentForcing::nmodes*xstep; kx++) - { - Real kappa = sqrt( (kx*kx)/(Lx*Lx) + (ky*ky)/(Ly*Ly) + (kz*kz)/(Lz*Lz) ); - - if (kappa <= kappaMax) - { - Real xT = cos(FTX(kx,ky,kz)*time + TAT(kx,ky,kz)); - - if ( TurbulentForcing::div_free_force ) - { - f1 += xT * - ( FAZ(kx,ky,kz)*TwoPi*(ky/Ly) - * sin(TwoPi*kx*x/Lx+FPZX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPZY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZZ(kx,ky,kz)) - - FAY(kx,ky,kz)*TwoPi*(kz/Lz) - * sin(TwoPi*kx*x/Lx+FPYX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPYY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPYZ(kx,ky,kz)) ); - - f2 += xT * - ( FAX(kx,ky,kz)*TwoPi*(kz/Lz) - * sin(TwoPi*kx*x/Lx+FPXX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPXY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPXZ(kx,ky,kz)) - - FAZ(kx,ky,kz)*TwoPi*(kx/Lx) - * cos(TwoPi*kx*x/Lx+FPZX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPZY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZZ(kx,ky,kz)) ); - - f3 += xT * - ( FAY(kx,ky,kz)*TwoPi*(kx/Lx) - * cos(TwoPi*kx*x/Lx+FPYX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPYY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPYZ(kx,ky,kz)) - - FAX(kx,ky,kz)*TwoPi*(ky/Ly) - * sin(TwoPi*kx*x/Lx+FPXX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPXY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPXZ(kx,ky,kz)) ); - } - else - { - - f1 += xT*FAX(kx,ky,kz)*cos(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - - f2 += xT*FAY(kx,ky,kz)*sin(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * cos(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * sin(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - - f3 += xT*FAZ(kx,ky,kz)*sin(TwoPi*kx*x/Lx+FPX(kx,ky,kz)) - * sin(TwoPi*ky*y/Ly+FPY(kx,ky,kz)) - * cos(TwoPi*kz*z/Lz+FPZ(kx,ky,kz)); - } - } - } - } - } - - frc(i,j,k,0) += dens(i,j,k,0) * f1; - frc(i,j,k,1) += dens(i,j,k,0) * f2; - frc(i,j,k,2) += dens(i,j,k,0) * f3; - - // if ( i==0 && j==0 && k==0 ){ - // printf("(0,0,0) : %15.13e %15.13e %15.13e\n", - // f1, f2, f3); - // } - // if ( i==16 && j==16 && k==16 ){ - // printf("(16,16,16) : %15.13e %15.13e %15.13e\n", - // f1, f2, f3); - // } - // if ( i==25 && j==12 && k==3 ){ - // printf("(25,12,3) : %15.13e %15.13e %15.13e\n", - // f1, f2, f3); - // } - - }); -#endif // Fast Force -#endif // Turbulent forcing - - } - - - // - // Scalar forcing - // - if ( scomp >= AMREX_SPACEDIM ) { - // Doing only scalars - force.setVal(0.0, bx, 0, ncomp); - - // - // Or create user-defined forcing. - // Recall we compute a density-weighted forcing term. - // - // auto const& frc = force.array(); - // amrex::ParallelFor(bx, ncomp, [frc] - // AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept - // { - // frc(i,j,k,n) = ; - // frc(i,j,k,n) *= rho; - // }); - } - else if ( scomp+ncomp > AMREX_SPACEDIM) { - // Doing scalars with vel - force.setVal(0.0, bx, Density, ncomp-Density); - - // - // Or create user-defined forcing. - // Recall we compute a density-weighted forcing term. - // - // auto const& frc = force.array(Density); - // amrex::ParallelFor(bx, ncomp-Density, [frc] - // AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept - // { - // frc(i,j,k,n) = ; - // frc(i,j,k,n) *= rho; - // }); - } - - if (ParallelDescriptor::IOProcessor() && getForceVerbose) { - Vector forcemin(ncomp); - Vector forcemax(ncomp); - for (int n=0; nforcemax[n]) forcemax[n] = f; - } - } - } -#if (AMREX_SPACEDIM == 3) - } -#endif - for (int n=0; n -#include -#include -#include -#include - -// factor by which to reduce sampling for faster performance -AMREX_GPU_MANAGED int TurbulentForcing::ff_factor; -// make the forcing divergence free? -AMREX_GPU_MANAGED bool TurbulentForcing::div_free_force; -// how many modes to use -AMREX_GPU_MANAGED int TurbulentForcing::nmodes; -// don't use any modes below mode_start. We probably don't need this -AMREX_GPU_MANAGED int TurbulentForcing::mode_start; -// Diagnostic print outs -AMREX_GPU_MANAGED int TurbulentForcing::verbose; - -amrex::Real* TurbulentForcing::forcedata; - - -void -TurbulentForcing::init_turbulent_forcing (const amrex::GpuArray& problo, const amrex::GpuArray& probhi) -{ - using namespace amrex; - - // Start with checks. - // This is for 3D only. - AMREX_ALWAYS_ASSERT(AMREX_SPACEDIM==3); - - // Forcing requires that Lx==Ly, Lz can be longer - Real Lx = probhi[0]-problo[0]; - Real Ly = probhi[1]-problo[1]; - Real Lz = probhi[2]-problo[2]; - AMREX_ALWAYS_ASSERT(Lx==Ly); - - // Read in parameters - ParmParse pp("turb"); - - nmodes = 4; - pp.get("nmodes", nmodes); - - div_free_force = true; - pp.query("div_free_force", div_free_force); - - ff_factor = 4; - pp.query("ff_factor", ff_factor); - - mode_start = 0; - pp.query("mode_start", mode_start); - - verbose = 0; - pp.query("verbose", verbose); - - // Inputs not yet defined. Could make runtime parameters if desired. - int hack_lz(0), spectrum_type(2), moderate_zero_modes(1); - Real forcing_time_scale_min(0.5), forcing_time_scale_max(1.0), force_scale(1.0); - - // tmp CPU storage that holds everything in one flat array - const int num_elmts=array_size*array_size*array_size; - const int tmp_size = num_fdarray*num_elmts; - Real tmp[tmp_size]; - - // Separate out forcing data into individual Array4's - int i_arr = 0; - int fd_ncomp = 1; - Dim3 fd_begin{0,0,0}; - Dim3 fd_end{array_size,array_size,array_size}; - - Array4 FTX(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 TAT(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPX(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPY(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPZ(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FAX(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FAY(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FAZ(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPXX(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPXY(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPXZ(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPYX(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPYY(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPYZ(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPZX(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPZY(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - Array4 FPZZ(&tmp[(i_arr++)*num_elmts], fd_begin, fd_end, fd_ncomp); - - if (hack_lz>0) { - if (hack_lz==1) { - Lz = Lz/2.0; - } else { - Lz = Lz/hack_lz; - } - } - - if (verbose) - { - Print() << "Lx = " << Lx << '\n'; - Print() << "Ly = " << Ly << '\n'; - Print() << "Lz = " << Lz << '\n'; - } - - Real Lmin = std::min(Lx,std::min(Ly,Lz)); - Real kappaMax = ((Real)nmodes)/Lmin + 1.0e-8; - int nxmodes = nmodes*(int)(0.5+Lx/Lmin); - int nymodes = nmodes*(int)(0.5+Ly/Lmin); - int nzmodes = nmodes*(int)(0.5+Lz/Lmin); - - // Mode indices run 0:n?modes, so they must fit in the 0:array_size-1 - // extent of the forcing arrays. - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::max(nxmodes,std::max(nymodes,nzmodes)) < array_size, - "TurbulentForcing: too many modes; reduce turb.nmodes or increase array_size"); - - if (verbose) - { - Print() << "Lmin = " << Lmin << '\n'; - Print() << "kappaMax = " << kappaMax << '\n'; - Print() << "nxmodes = " << nxmodes << '\n'; - Print() << "nymodes = " << nymodes << '\n'; - Print() << "nzmodes = " << nzmodes << '\n'; - } - - Real freqMin = 1.0/forcing_time_scale_max; - Real freqMax = 1.0/forcing_time_scale_min; - Real freqDiff= freqMax-freqMin; - - if (verbose) - { - Print() << "forcing_time_scale_min = " << forcing_time_scale_min << '\n'; - Print() << "forcing_time_scale_max = " << forcing_time_scale_max << '\n'; - Print() << "freqMin = " << freqMin << '\n'; - Print() << "freqMax = " << freqMax << '\n'; - Print() << "freqDiff = " << freqDiff << '\n'; - } - - // initiate the magic - DepRand::InitRandom((unsigned long)111397); - - int mode_count = 0; - - int xstep = (int)(Lx/Lmin+0.5); - int ystep = (int)(Ly/Lmin+0.5); - int zstep = (int)(Lz/Lmin+0.5); - - if (verbose) - Print() << "Mode step = " << xstep << " " << ystep << " " << zstep << '\n'; - - for (int kz = mode_start*zstep; kz <= nzmodes; kz += zstep ) { - Real kzd = (Real)kz; - for (int ky = mode_start*ystep; ky <= nymodes; ky += ystep ) { - Real kyd = (Real)ky; - for (int kx = mode_start*xstep; kx <= nxmodes; kx += xstep ) { - Real kxd = (Real)kx; - - Real kappa = sqrt( (kxd*kxd)/(Lx*Lx) + (kyd*kyd)/(Ly*Ly) + (kzd*kzd)/(Lz*Lz) ); - - if (kappa<=kappaMax) { - FTX(kx,ky,kz) = (freqMin + freqDiff*DepRand::Random() )*TwoPi; - // Translation angles, theta=0..2Pi and phi=0..Pi - TAT(kx,ky,kz) = DepRand::Random()*TwoPi; - // Phases - FPX(kx,ky,kz) = DepRand::Random()*TwoPi; - FPY(kx,ky,kz) = DepRand::Random()*TwoPi; - FPZ(kx,ky,kz) = DepRand::Random()*TwoPi; - if (div_free_force==1) { - FPXX(kx,ky,kz) = DepRand::Random()*TwoPi; - FPYX(kx,ky,kz) = DepRand::Random()*TwoPi; - FPZX(kx,ky,kz) = DepRand::Random()*TwoPi; - FPXY(kx,ky,kz) = DepRand::Random()*TwoPi; - FPYY(kx,ky,kz) = DepRand::Random()*TwoPi; - FPZY(kx,ky,kz) = DepRand::Random()*TwoPi; - FPXZ(kx,ky,kz) = DepRand::Random()*TwoPi; - FPYZ(kx,ky,kz) = DepRand::Random()*TwoPi; - FPZZ(kx,ky,kz) = DepRand::Random()*TwoPi; - } - // Amplitudes (alpha) - Real thetaTmp = DepRand::Random()*TwoPi; - Real cosThetaTmp = cos(thetaTmp); - Real sinThetaTmp = sin(thetaTmp); - - Real phiTmp = DepRand::Random()*Pi; - Real cosPhiTmp = cos(phiTmp); - Real sinPhiTmp = sin(phiTmp); - - Real px = cosThetaTmp * sinPhiTmp; - Real py = sinThetaTmp * sinPhiTmp; - Real pz = cosPhiTmp; - - Real mp2 = px*px + py*py + pz*pz; - if (kappa < 0.000001) { - Print() << "ZERO AMPLITUDE MODE " << kx << ky << kz << '\n'; - FAX(kx,ky,kz) = 0.; - FAY(kx,ky,kz) = 0.; - FAZ(kx,ky,kz) = 0.; - } else { - // Count modes that contribute - mode_count++; - // Set amplitudes - Real Ekh; - if (spectrum_type==1) { - Ekh = 1. / kappa; - } else if (spectrum_type==2) { - Ekh = 1. / (kappa*kappa); - } else { - Ekh = 1.; - } - if (div_free_force==1) { - Ekh /= kappa; - } - if (moderate_zero_modes==1) { - if (kx==0) Ekh /= 2.; - if (ky==0) Ekh /= 2.; - if (kz==0) Ekh /= 2.; - } - if (force_scale>0.) { - FAX(kx,ky,kz) = force_scale * px * Ekh / mp2; - FAY(kx,ky,kz) = force_scale * py * Ekh / mp2; - FAZ(kx,ky,kz) = force_scale * pz * Ekh / mp2; - } else { - FAX(kx,ky,kz) = px * Ekh / mp2; - FAY(kx,ky,kz) = py * Ekh / mp2; - FAZ(kx,ky,kz) = pz * Ekh / mp2; - } - - if (verbose) - { - Print() << "Mode"; - Print() << "kappa = " << kx << " " << ky << " " << kz << " " << kappa << " " - << sqrt(FAX(kx,ky,kz)*FAX(kx,ky,kz)+FAY(kx,ky,kz)*FAY(kx,ky,kz)+FAZ(kx,ky,kz)*FAZ(kx,ky,kz)) << '\n'; - Print() << "Amplitudes - A" << '\n'; - Print() << FAX(kx,ky,kz) << " " << FAY(kx,ky,kz) << " " << FAZ(kx,ky,kz) << '\n'; - Print() << "Frequencies" << '\n'; - Print() << FTX(kx,ky,kz) << '\n'; - Print() << "TAT" << '\n'; - Print() << TAT(kx,ky,kz) << '\n'; - Print() << "Amplitudes - AA" << '\n'; - Print() << FPXX(kx,ky,kz) << " " << FPYX(kx,ky,kz) << " " << FPZX(kx,ky,kz) << '\n'; - Print() << FPXY(kx,ky,kz) << " " << FPYY(kx,ky,kz) << " " << FPZY(kx,ky,kz) << '\n'; - Print() << FPXZ(kx,ky,kz) << " " << FPYZ(kx,ky,kz) << " " << FPZZ(kx,ky,kz) << '\n'; - } - } - } - } - } - } - - // Now let's break symmetry, have to assume high aspect ratio in z for now - int reduced_mode_count = 0; - - for (int kz = 1; kz < zstep; kz++ ) { - Real kzd = (Real)kz; - for (int ky = mode_start; ky <= nymodes; ky += ystep ) { - Real kyd = (Real)ky; - for (int kx = mode_start; kx <= nxmodes; kx += xstep ) { - Real kxd = (Real)kx; - - Real kappa = sqrt( (kxd*kxd)/(Lx*Lx) + (kyd*kyd)/(Ly*Ly) + (kzd*kzd)/(Lz*Lz) ); - - if (kappa<=kappaMax) { - FTX(kx,ky,kz) = (freqMin + freqDiff*DepRand::Random() )*TwoPi; - // Translation angles, theta=0..2Pi and phi=0..Pi - TAT(kx,ky,kz) = DepRand::Random()*TwoPi; - // Phases - FPX(kx,ky,kz) = DepRand::Random()*TwoPi; - FPY(kx,ky,kz) = DepRand::Random()*TwoPi; - FPZ(kx,ky,kz) = DepRand::Random()*TwoPi; - if (div_free_force==1) { - FPXX(kx,ky,kz) = DepRand::Random()*TwoPi; - FPYX(kx,ky,kz) = DepRand::Random()*TwoPi; - FPZX(kx,ky,kz) = DepRand::Random()*TwoPi; - FPXY(kx,ky,kz) = DepRand::Random()*TwoPi; - FPYY(kx,ky,kz) = DepRand::Random()*TwoPi; - FPZY(kx,ky,kz) = DepRand::Random()*TwoPi; - FPXZ(kx,ky,kz) = DepRand::Random()*TwoPi; - FPYZ(kx,ky,kz) = DepRand::Random()*TwoPi; - FPZZ(kx,ky,kz) = DepRand::Random()*TwoPi; - } - // Amplitudes (alpha) - Real thetaTmp = DepRand::Random()*TwoPi; - Real cosThetaTmp = cos(thetaTmp); - Real sinThetaTmp = sin(thetaTmp); - - Real phiTmp = DepRand::Random()*Pi; - Real cosPhiTmp = cos(phiTmp); - Real sinPhiTmp = sin(phiTmp); - - Real px = cosThetaTmp * sinPhiTmp; - Real py = sinThetaTmp * sinPhiTmp; - Real pz = cosPhiTmp; - - Real mp2 = px*px + py*py + pz*pz; - if (kappa < 0.000001) { - Print() << "ZERO AMPLITUDE MODE " << kx << ky << kz << '\n'; - FAX(kx,ky,kz) = 0.; - FAY(kx,ky,kz) = 0.; - FAZ(kx,ky,kz) = 0.; - } else { - // Count modes that contribute - reduced_mode_count++; - // Set amplitudes - Real Ekh; - if (spectrum_type==1) { - Ekh = 1. / kappa; - } else if (spectrum_type==2) { - Ekh = 1. / (kappa*kappa); - } else { - Ekh = 1.; - } - if (div_free_force==1) { - Ekh /= kappa; - } - if (moderate_zero_modes==1) { - if (kx==0) Ekh /= 2.; - if (ky==0) Ekh /= 2.; - if (kz==0) Ekh /= 2.; - } - if (force_scale>0.) { - FAX(kx,ky,kz) = force_scale * px * Ekh / mp2; - FAY(kx,ky,kz) = force_scale * py * Ekh / mp2; - FAZ(kx,ky,kz) = force_scale * pz * Ekh / mp2; - } else { - FAX(kx,ky,kz) = px * Ekh / mp2; - FAY(kx,ky,kz) = py * Ekh / mp2; - FAZ(kx,ky,kz) = pz * Ekh / mp2; - } - - if (verbose) - { - Print() << "Mode"; - Print() << "kappa = " << kx << " " << ky << " " << kz << " " << kappa << " " - << sqrt(FAX(kx,ky,kz)*FAX(kx,ky,kz)+FAY(kx,ky,kz)*FAY(kx,ky,kz)+FAZ(kx,ky,kz)*FAZ(kx,ky,kz)) << '\n'; - Print() << "Amplitudes - A" << '\n'; - Print() << FAX(kx,ky,kz) << " " << FAY(kx,ky,kz) << " " << FAZ(kx,ky,kz) << '\n'; - Print() << "Frequencies" << '\n'; - Print() << FTX(kx,ky,kz) << '\n'; - Print() << "TAT" << '\n'; - Print() << TAT(kx,ky,kz) << '\n'; - Print() << "Amplitudes - AA" << '\n'; - Print() << FPXX(kx,ky,kz) << " " << FPYX(kx,ky,kz) << " " << FPZX(kx,ky,kz) << '\n'; - Print() << FPXY(kx,ky,kz) << " " << FPYY(kx,ky,kz) << " " << FPZY(kx,ky,kz) << '\n'; - Print() << FPXZ(kx,ky,kz) << " " << FPYZ(kx,ky,kz) << " " << FPZZ(kx,ky,kz) << '\n'; - } - } - } - } - } - } - - Print() << "mode_count = " << mode_count << '\n'; - Print() << "reduced_mode_count = " << reduced_mode_count << '\n'; - if (spectrum_type==1) { - Print() << "Spectrum type 1" << '\n'; - } else if (spectrum_type==2) { - Print() << "Spectrum type 2" << '\n'; - } else { - Print() << "Spectrum type OTHER" << '\n'; - } - -// Now allocate forcedata and copy in tmp array. -#ifdef AMREX_USE_GPU - if (Gpu::inLaunchRegion()) - { - forcedata = static_cast(The_Arena()->alloc(tmp_size*sizeof(Real))); - // Synchronous copy: tmp is a stack array that dies when this - // function returns, so the copy must complete before then. - Gpu::htod_memcpy(forcedata, tmp, tmp_size*sizeof(Real)); - } - else -#endif - { - forcedata = static_cast(The_Pinned_Arena()->alloc(tmp_size*sizeof(Real))); - std::memcpy(forcedata, tmp, tmp_size*sizeof(Real)); - } -} diff --git a/Tutorials/HIT/TurbulentForcing_params.H b/Tutorials/HIT/TurbulentForcing_params.H deleted file mode 100644 index cdd1d5fd7..000000000 --- a/Tutorials/HIT/TurbulentForcing_params.H +++ /dev/null @@ -1,24 +0,0 @@ -#ifndef IAMR_TurbulentForcing_params_H_ -#define IAMR_TurbulentForcing_params_H_ - -namespace TurbulentForcing { - - // function to generate the parameters - void init_turbulent_forcing (const amrex::GpuArray& problo, const amrex::GpuArray& probhi); - - extern AMREX_GPU_MANAGED int verbose; - // factor by which to reduce sampling for faster performance - extern AMREX_GPU_MANAGED int ff_factor; - // make the forcing divergence free? - extern AMREX_GPU_MANAGED bool div_free_force; - // how many modes to use - extern AMREX_GPU_MANAGED int nmodes; - // don't use any modes below mode_start. We probably don't need this - extern AMREX_GPU_MANAGED int mode_start; - - constexpr int array_size = 33; - constexpr int num_fdarray = 17; - // forcedata will contain num_fdarray arrays of size (0,0,0)(array_size-1,array_size-1,array_size-1) - extern amrex::Real* forcedata; -} -#endif diff --git a/Tutorials/HIT/depRand.H b/Tutorials/HIT/depRand.H deleted file mode 100644 index fd58b65e3..000000000 --- a/Tutorials/HIT/depRand.H +++ /dev/null @@ -1,113 +0,0 @@ - -#ifndef IAMR_depRand_H_ -#define IAMR_depRand_H_ - -#include - -namespace DepRand -{ - // - // The Mersenne twistor : - // - class mt19937 - { - public: - typedef unsigned long seed_type; - - explicit mt19937 (seed_type seed = 4357UL); - mt19937 (seed_type seed, int numprocs); - mt19937 (seed_type array[], int array_len); - void rewind(); - void reset(unsigned long seed); - - double d_value (); // [0,1] random numbers - double d1_value (); // [0,1) random numbers - double d2_value (); // (0,1) random numbers - - long l_value (); // [0,2^31-1] random numbers - unsigned long u_value (); // [0,2^32-1] random numbers - - void save (amrex::Vector& state) const; - int RNGstatesize() const; - void restore (const amrex::Vector& state); - private: - enum { N = 624 }; - static unsigned long init_seed; - static unsigned long mt[N]; // the array for the state vector - static int mti; // mti==N+1 means mt[N] is not initialized -#ifdef _OPENMP -#pragma omp threadprivate(init_seed,mt,mti) -#endif - private: - void sgenrand (unsigned long seed); - void sgenrand (seed_type seed_array[], int len); - unsigned long igenrand (); - void reload (); - }; - - /* - Mersenne Twister pseudo-random number generator. - - Generates one pseudorandom real number (double) which is - uniformly distributed on [0,1]-interval for each call. - - Accepts any 32-bit integer as a seed -- uses 4357 as the default. - - Has a period of 2**19937. - - Mersenne Twister Home Page: http://www.math.keio.ac.jp/matumoto/emt.html - - There is also an entry point for Fortran callable as: - - REAL_T rn - call blutilrand(rn) - - Internally, blutilrand() calls a static Mersenne Twister object (the - same one used by AMReX::Random()) to get a value in [0,1] and then - sets "rn" to that value. - */ - double Random (); // [0,1] - double Random1 (); // [0,1) - double Random2 (); // (0,1) - unsigned long Random_int(unsigned long n); // [0,n-1], where n<=2^32-1 - /* Set the seed of the random number generator. - - There is also an entry point for Fortran callable as: - - INTEGER seed - call blutilinitrand(seed) - - or - - INTEGER seed - call blinitrand(seed) - */ - void InitRandom (unsigned long seed); - void InitRandom (unsigned long seed, int numprocs); - - void ResetRandomSeed(unsigned long seed); - // - // Save and restore random state. - // - // state.size() == 626 on return from Save & on entry to Restore. - // - void SaveRandomState (amrex::Vector& state); - - int sizeofRandomState (); - - void RestoreRandomState (const amrex::Vector& state); - // - // Create a unique subset of random numbers from a pool - // of integers in the range [0, poolSize - 1] - // the set will be in the order they are found - // setSize must be <= poolSize - // uSet will be resized to setSize - // if you want all processors to have the same set, - // call this on one processor and broadcast the array - // - void UniqueRandomSubset (amrex::Vector &uSet, int setSize, int poolSize, - bool printSet = false); - -} - -#endif diff --git a/Tutorials/HIT/depRand.cpp b/Tutorials/HIT/depRand.cpp deleted file mode 100644 index e994731d0..000000000 --- a/Tutorials/HIT/depRand.cpp +++ /dev/null @@ -1,371 +0,0 @@ -#include -#include - -#include - -using namespace amrex; - -// -// AMReX Interface to Mersenne Twistor -// - -/* A C-program for MT19937: Real number version (1999/10/28) */ -/* genrand() generates one pseudorandom real number (double) */ -/* which is uniformly distributed on [0,1]-interval, for each */ -/* call. sgenrand(seed) sets initial values to the working area */ -/* of 624 words. Before genrand(), sgenrand(seed) must be */ -/* called once. (seed is any 32-bit integer.) */ -/* Integer generator is obtained by modifying two lines. */ -/* Coded by Takuji Nishimura, considering the suggestions by */ -/* Topher Cooper and Marc Rieffel in July-Aug. 1997. */ - -/* This library is free software under the Artistic license: */ -/* see the file COPYING distributed together with this code. */ -/* For the verification of the code, its output sequence file */ -/* mt19937-1.out is attached (2001/4/2) */ - -/* Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura. */ -/* Any feedback is very welcome. For any question, comments, */ -/* see http://www.math.keio.ac.jp/matumoto/emt.html or email */ -/* matumoto@math.keio.ac.jp */ - -/* REFERENCE */ -/* M. Matsumoto and T. Nishimura, */ -/* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform */ -/* Pseudo-Random Number Generator", */ -/* ACM Transactions on Modeling and Computer Simulation, */ -/* Vol. 8, No. 1, January 1998, pp 3--30. */ - -unsigned long DepRand::mt19937::init_seed; -unsigned long DepRand::mt19937::mt[DepRand::mt19937::N]; -int DepRand::mt19937::mti; - -// -// initializing with a NONZERO seed. -// -void -DepRand::mt19937::sgenrand(unsigned long seed) -{ - mt[0]= seed & 0xffffffffUL; - for ( mti=1; mti> 30L)) + mti); - /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ - /* In the previous versions, MSBs of the seed affect */ - /* only MSBs of the array mt[]. */ - /* 2002/01/09 modified by Makoto Matsumoto */ - mt[mti] &= 0xffffffffUL; /* for >32 bit machines */ - } -} - -/* initialize by an array with array-length */ -/* init_key is the array for initializing keys */ -/* key_length is its length */ -void -DepRand::mt19937::sgenrand(unsigned long init_key[], int key_length) -{ - int i, j, k; - sgenrand(19650218UL); - i=1; j=0; - k = (N>key_length ? N : key_length); - for ( ; k; k-- ) - { - mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + init_key[j] + j; /* non linear */ - mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ - i++; j++; - if (i>=N) { mt[0] = mt[N-1]; i=1; } - if (j>=key_length) j=0; - } - for ( k=N-1; k; k-- ) - { - mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - i; /* non linear */ - mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ - i++; - if (i>=N) { mt[0] = mt[N-1]; i=1; } - } - - mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ -} - -void -DepRand::mt19937::reload() -{ - unsigned long y; - int kk; - - const int M = 397; - -#define MATRIX_A 0x9908B0DFUL // Constant vector a -#define UPPER_MASK 0x80000000UL // Most significant w-r bits -#define LOWER_MASK 0x7FFFFFFFUL // least significant r bits - // - // mag01[x] = x * MATRIX_A for x=0,1 - // - static unsigned long mag01[2]={0x0UL, MATRIX_A}; - for ( kk=0; kk> 1L) ^ mag01[y & 0x1UL]; - } - for ( ; kk> 1L) ^ mag01[y & 0x1UL]; - } - y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); - mt[N-1] = mt[M-1] ^ (y >> 1L) ^ mag01[y & 0x1UL]; - - mti = 0; - -#undef MATRIX_A -#undef UPPER_MASK -#undef LOWER_MASK -} - -unsigned long -DepRand::mt19937::igenrand() -{ - // - // Generate N words at one time. - // - if ( mti >= N ) reload(); - - unsigned long y = mt[mti++]; - - /* Tempering */ - y ^= (y >> 11); - y ^= (y << 7) & 0x9d2c5680UL; - y ^= (y << 15) & 0xefc60000UL; - y ^= (y >> 18); - - return y; -} - -DepRand::mt19937::mt19937(unsigned long seed) -{ - init_seed = seed; - mti = N; - sgenrand(seed); -} - -DepRand::mt19937::mt19937(unsigned long seed, int numprocs) -{ -#ifdef _OPENMP -#pragma omp parallel - { - init_seed = seed + omp_get_thread_num() * numprocs; - mti = N; - sgenrand(init_seed); - } -#else - init_seed = seed; - mti = N; - sgenrand(init_seed); -#endif -} - -DepRand::mt19937::mt19937 (unsigned long seed_array[], int len) -{ - sgenrand(seed_array, len); -} - -void -DepRand::mt19937::rewind() -{ - sgenrand(init_seed); -} - -void -DepRand::mt19937::reset(unsigned long seed) -{ - sgenrand(seed); -} - -// -// [0,1] random numbers -// -double -DepRand::mt19937::d_value() -{ - return double(igenrand()) * (1.0/4294967295.0); // divided by 2^32-1 -} - -// -// [0,1) random numbers -// -double -DepRand::mt19937::d1_value() -{ - return double(igenrand()) * (1.0/4294967296.0); // divided by 2^32 -} - -// -// (0,1) random numbers -// -double -DepRand::mt19937::d2_value() -{ - return (double(igenrand()) + .5) * (1.0/4294967296.0); // divided by 2^32 -} - -long -DepRand::mt19937::l_value() -{ - return (long)(igenrand()>>1); -} - -unsigned long -DepRand::mt19937::u_value() -{ - return igenrand(); -} - -void -DepRand::mt19937::save (Vector& state) const -{ - state.resize(N+2); - state[0] = init_seed; - for (int i = 0; i < N; i++) - state[i+1] = mt[i]; - state[N+1] = mti; -} - -int -DepRand::mt19937::RNGstatesize () const -{ - return N+2; -} - -void -DepRand::mt19937::restore (const Vector& state) -{ - if (state.size() != N+2) - Error("mt19937::restore(): incorrectly sized state vector"); - - init_seed = state[0]; - for (int i = 0; i < N; i++) - mt[i] = state[i+1]; - mti = state[N+1]; - - if (mti < 0 || mti > N) - Error("mt19937::restore(): mti out-of-bounds"); -} - -namespace -{ - DepRand::mt19937 the_generator; -} - -void -DepRand::InitRandom (unsigned long seed) -{ - the_generator = mt19937(seed); -} - -void -DepRand::InitRandom (unsigned long seed, int numprocs) -{ - the_generator = mt19937(seed, numprocs); -} - -void DepRand::ResetRandomSeed(unsigned long seed) -{ - the_generator.reset(seed); -} - -double -DepRand::Random () -{ - return the_generator.d_value(); -} - -double -DepRand::Random1 () -{ - return the_generator.d1_value(); -} - -double -DepRand::Random2 () -{ - return the_generator.d2_value(); -} - -unsigned long -DepRand::Random_int(unsigned long n) -{ - const unsigned long umax = 4294967295UL; // 2^32-1 - BL_ASSERT( n > 0 && n <= umax ); - unsigned long scale = umax/n; - unsigned long r; - do { - r = the_generator.u_value() / scale; - } while (r >= n); - return r; -} - -void -DepRand::SaveRandomState (Vector& state) -{ - the_generator.save(state); -} - -int -DepRand::sizeofRandomState () -{ - return the_generator.RNGstatesize(); -} - -void -DepRand::RestoreRandomState (const Vector& state) -{ - the_generator.restore(state); -} - -void -DepRand::UniqueRandomSubset (Vector &uSet, int setSize, int poolSize, - bool printSet) -{ - if(setSize > poolSize) { - Abort("**** Error in UniqueRandomSubset: setSize > poolSize."); - } - std::set copySet; - Vector uSetTemp; - while(copySet.size() < (unsigned long)setSize) { - int r(DepRand::Random_int(poolSize)); - if(copySet.find(r) == copySet.end()) { - copySet.insert(r); - uSetTemp.push_back(r); - } - } - uSet = uSetTemp; - if(printSet) { - for(int i(0); i < uSet.size(); ++i) { - std::cout << "uSet[" << i << "] = " << uSet[i] << '\n'; - } - } -} - - -#if 0 -// -// Fortran entry points for DepRand::Random(). -// - -BL_FORT_PROC_DECL(BLUTILINITRAND,blutilinitrand)(const int* sd) -{ - unsigned long seed = *sd; - DepRand::InitRandom(seed); -} - -BL_FORT_PROC_DECL(BLINITRAND,blinitrand)(const int* sd) -{ - unsigned long seed = *sd; - DepRand::InitRandom(seed); -} - -BL_FORT_PROC_DECL(BLUTILRAND,blutilrand)(Real* rn) -{ - *rn = DepRand::Random(); -} -#endif diff --git a/Tutorials/HIT/derivespect-inputs b/Tutorials/HIT/derivespect-inputs deleted file mode 100644 index 0d9fdfaab..000000000 --- a/Tutorials/HIT/derivespect-inputs +++ /dev/null @@ -1,42 +0,0 @@ -##----------------------------------------------------------- -## INPUT PARAMETERS FOR AMRDERIVESPECTRUM PROBLEM -##----------------------------------------------------------- - -# Verbose flag -verbose = 1 - -# List of input plotfiles to process -infile = tmp_plotfile - -# Variables to read from the input plotfiles -vars = density x_vort y_vort z_vort x_velocity y_velocity z_velocity divu - -# Flag determining whether div_free -div_free = 0 - -# Flag determining whether or not to transpose spatial indices -transpose_dp = 0 - -# Flag sets whether or not to weight by density^(1/3) -density_weighting = 0 - -# Density field -density = density - -# Flag sets whether or not to apply cutoff density -use_cutoff_density = 0 - -# Cutoff density below which to zero the velocities -cutoff_density = 0.000 - -# Flag determining whether or not to read a list of -# wavenumbers to filter on. -do_filter = 0 - -# List of wavenumbers to filter on -# filterWN = 0.0 0.0 0.0 - -# This sets the finest level to do the FFT -finestLevel = 2 - - diff --git a/Tutorials/HIT/gen_hit_ic.py b/Tutorials/HIT/gen_hit_ic.py deleted file mode 100755 index 3d8472373..000000000 --- a/Tutorials/HIT/gen_hit_ic.py +++ /dev/null @@ -1,463 +0,0 @@ -#!/usr/bin/env python -# -# Generate a table of the velocity fluctuations for the homogeneous -# isotropic turbulence case at a specific k0 (default to 4) -# -# Order of operations: -# 1. velocity fluctuations generated on a 512^3 grid in wavenumber space -# 2. Coefficients associated to wavenumbers that cannot be represented on -# the desired grid are set to 0 (sharp wavenumber cutoff) -# 3. inverse Fourier transform of the velocity fluctuations (512^3 grid) -# 4. velocity fluctuations resampled on the desired grid (N^3) -# -# The velocity fluctuations are normalized by urms0 so to get the -# actual velocity fluctuations, one must multiply these velocities by -# the appropriate urms0. -# -# - -# ======================================================================== -# -# Imports -# -# ======================================================================== -import argparse -import sys -import time -from datetime import timedelta -import numpy as np -import scipy.interpolate as spi -import matplotlib as mpl - -mpl.use("Agg") -import matplotlib.pyplot as plt - - -# ======================================================================== -# -# Parse arguments -# -# ======================================================================== -parser = argparse.ArgumentParser( - description="Generate the velocity fluctuations for the HIT IC" -) -parser.add_argument( - "-k0", help="Wave number containing highest energy", type=float, default=4.0 -) -parser.add_argument("-N", help="Resolution", type=int, default=16) -parser.add_argument( - "-s", "--seed", help="Random number generator seed", type=int, default=42 -) -parser.add_argument( - "-p", "--plot", help="Save a plot of the x-velocity", action="store_true" -) -args = parser.parse_args() - -# =============================================================================== -# -# Some defaults variables -# -# =============================================================================== -plt.rc("text", usetex=True) -plt.rc("font", family="serif", serif="Times") -cmap_med = [ - "#F15A60", - "#7AC36A", - "#5A9BD4", - "#FAA75B", - "#9E67AB", - "#CE7058", - "#D77FB4", - "#737373", -] -cmap = [ - "#EE2E2F", - "#008C48", - "#185AA9", - "#F47D23", - "#662C91", - "#A21D21", - "#B43894", - "#010202", -] -dashseq = [ - (None, None), - [10, 5], - [10, 4, 3, 4], - [3, 3], - [10, 4, 3, 4, 3, 4], - [3, 3], - [3, 3], -] -markertype = ["s", "d", "o", "p", "h"] - -# ======================================================================== -# -# Function definitions -# -# ======================================================================== -def div0(a, b): - """ Ignore division by 0, just replace it by 0, - - From: http://stackoverflow.com/questions/26248654/numpy-return-0-with-divide-by-zero - e.g. div0( [-1, 0, 1], 0 ) -> [0, 0, 0] - """ - with np.errstate(divide="ignore", invalid="ignore"): - c = np.true_divide(a, b) - c[~np.isfinite(c)] = 0 # -inf inf NaN - return c - - -# ======================================================================== -def abs2(x): - """This is equivalent to np.abs(x)**2 or x*np.conj(x) - - To make it faster, add this right before the function definition - import numba - @numba.vectorize([numba.float64(numba.complex128),numba.float32(numba.complex64)]) - """ - return x.real ** 2 + x.imag ** 2 - - -# ======================================================================== -# -# Main -# -# ======================================================================== - -# Timer -start = time.time() - -# ======================================================================== -# 1. velocity fluctuations generated on a 512^3 grid in wavenumber space - -# Dimension of the large cube -N = 512 -halfN = int(round(0.5 * N)) -xs = 0 -xe = 2.0 * np.pi -L = xe - xs -dx = L / N - -# Only work if N and args.N are even -if not ((args.N % 2 == 0) and N % 2 == 0): - print("N or args.N is not even. Exiting") - sys.exit(1) - -# Get cell centered values and meshed grid -x = np.linspace(xs, xe, N + 1) -xc = (x[1:] + x[:-1]) / 2 # get cell center coordinates -X, Y, Z = np.meshgrid(xc, xc, xc, indexing="ij") - -# Get the wave numbers and associated quantities -k = np.concatenate((np.arange(halfN), np.arange(-halfN, 0, 1)), axis=0) -khalf = np.arange(halfN + 1) -k1, k2, k3 = np.meshgrid(k, k, khalf, indexing="ij") -kmag = np.sqrt(k1 ** 2 + k2 ** 2 + k3 ** 2) -k12 = np.sqrt(k1 ** 2 + k2 ** 2) -k1k12 = div0(k1, k12) -k2k12 = div0(k2, k12) -k3kmag = div0(k3, kmag) -k12kmag = div0(k12, kmag) - -# Generate data - -# # Toy Fourier data corresponding to uo = cos(X) * cos(2*Y) * cos(3*Z) -# uo = np.cos(X) * np.cos(2*Y) * np.cos(3*Z) -# uf = np.fft.rfftn(uo) -# vf = np.copy(uf) -# wf = np.copy(uf) - -# Energy spectrum -Ek = ( - 16.0 - * np.sqrt(2.0 / np.pi) - * (kmag ** 4) - / (args.k0 ** 5) - * np.exp(-2.0 * (kmag ** 2) / (args.k0 ** 2)) -) - -# Draw random numbers -np.random.seed(args.seed) -phi1 = np.random.uniform(0, 2 * np.pi, np.shape(kmag)) -phi2 = np.random.uniform(0, 2 * np.pi, np.shape(kmag)) -phi3 = np.random.uniform(0, 2 * np.pi, np.shape(kmag)) - -# the random quantities -prefix = np.sqrt(2.0 * div0(Ek, 4.0 * np.pi * (kmag ** 2))) -a = prefix * np.exp(1j * phi1) * np.cos(phi3) -b = prefix * np.exp(1j * phi2) * np.sin(phi3) - -# the random velocities -uf = k2k12 * a + k1k12 * k3kmag * b -vf = k2k12 * k3kmag * b - k1k12 * a -wf = -k12kmag * b - -# Impose the 3D spherical symmetry (to ensure we have a real signal) -# equiv: uf[-l,-m,0] = np.conj(uf[ l, m,0]) for l=0..N/2 and m=0..N/2 -uf[N:halfN:-1, N:halfN:-1, 0] = np.conj(uf[1:halfN, 1:halfN, 0]) -# symmetry on first column -uf[N:halfN:-1, 0, 0] = np.conj(uf[1:halfN, 0, 0]) -# symmetry on first row -uf[0, N:halfN:-1, 0] = np.conj(uf[0, 1:halfN, 0]) -# symmetry about the (halfN,halfN) element -uf[halfN - 1 : 0 : -1, N : halfN - 1 : -1, 0] = np.conj( - uf[halfN + 1 : N, 1 : halfN + 1, 0] -) - -vf[N:halfN:-1, N:halfN:-1, 0] = np.conj(vf[1:halfN, 1:halfN, 0]) -vf[halfN - 1 : 0 : -1, N : halfN - 1 : -1, 0] = np.conj( - vf[halfN + 1 : N, 1 : halfN + 1, 0] -) -vf[N:halfN:-1, 0, 0] = np.conj(vf[1:halfN, 0, 0]) -vf[0, N:halfN:-1, 0] = np.conj(vf[0, 1:halfN, 0]) - -wf[N:halfN:-1, N:halfN:-1, 0] = np.conj(wf[1:halfN, 1:halfN, 0]) -wf[halfN - 1 : 0 : -1, N : halfN - 1 : -1, 0] = np.conj( - wf[halfN + 1 : N, 1 : halfN + 1, 0] -) -wf[N:halfN:-1, 0, 0] = np.conj(wf[1:halfN, 0, 0]) -wf[0, N:halfN:-1, 0] = np.conj(wf[0, 1:halfN, 0]) - -# Normalize. Because we are generating the data in wavenumber space, -# we have to multiply by N**3 because in the definition of the numpy -# ifftn there is a 1/N**n. -uf = uf * N ** 3 -vf = vf * N ** 3 -wf = wf * N ** 3 - -# # Quick check on energy content (make sure you add both the current -# # contribution and the one we are neglecting because we are assuming -# # real input data) -# print('Energy = int E(k) dk = 0.5 * int (uf**2 + vf**2 wf**2) dk1 dk2 dk3 = {0:.10f} ~= 3/2'.format( -# (np.sum(abs2(uf ) + abs2(vf ) + abs2(wf )) + -# np.sum(abs2(uf[:,:,1:-1]) + abs2(vf[:,:,1:-1]) + abs2(wf[:,:,1:-1]))) -# * 0.5 / N**6)) - -# if plotting, save the original field (before filtering) -if args.plot: - uo = np.fft.irfftn(uf) - Eko = ( - 16.0 - * np.sqrt(2.0 / np.pi) - * (khalf ** 4) - / (args.k0 ** 5) - * np.exp(-2.0 * (khalf ** 2) / (args.k0 ** 2)) - ) - - # Get the spectrum from 3D velocity field - kbins = np.arange(1, halfN + 1) - Nbins = len(kbins) - whichbin = np.digitize(kmag.flat, kbins) - ncount = np.bincount(whichbin) - - KI = (abs2(uf) + abs2(vf) + abs2(wf)) * 0.5 / N ** 6 - KI[:, :, 1:-1] += ( - (abs2(uf[:, :, 1:-1]) + abs2(vf[:, :, 1:-1]) + abs2(wf[:, :, 1:-1])) - * 0.5 - / N ** 6 - ) - - Eku = np.zeros(len(ncount) - 1) - for n in range(1, len(ncount)): - Eku[n - 1] = np.sum(KI.flat[whichbin == n]) - - ku = 0.5 * (kbins[0 : Nbins - 1] + kbins[1:Nbins]) + 1 - Eku = Eku[1:Nbins] - - -# ======================================================================== -# 2. Coefficients associated to wavenumbers that cannot be represented -# on the desired grid are set to 0 (sharp wavenumber cutoff) -kmagc = 0.5 * args.N -uf[kmag > kmagc] = 0.0 -vf[kmag > kmagc] = 0.0 -wf[kmag > kmagc] = 0.0 - - -# ======================================================================== -# 3. inverse Fourier transform of the velocity fluctuations (512^3 grid) -u = np.fft.irfftn(uf, s=(N, N, N)) -v = np.fft.irfftn(vf, s=(N, N, N)) -w = np.fft.irfftn(wf, s=(N, N, N)) - -# Another energy content check -print( - "Energy = 1/V * int E(x,y,z) dV = 0.5/V * int (u**2 + v**2 + w**2) dx dy dz = {0:.10f} ~= 3/2".format( - np.sum(u ** 2 + v ** 2 + w ** 2) * 0.5 * (dx / L) ** 3 - ) -) - -# # Enstrophy check -# _, dudy, dudz = np.gradient(u, dx) -# dvdx, _, dvdz = np.gradient(v, dx) -# dwdx, dwdy, _ = np.gradient(w, dx) -# wx = dwdy-dvdz -# wy = dudz-dwdx -# wz = dvdx-dudy -# lambda0 = 2.0/args.k0 -# print('Enstrophy = 0.5/V * int (wx**2 + wy**2 + wz**2) dx dy dz= -# {0:.10f} ~= '.format(np.sum(wx**2+wy**2+wz**2) * 0.5 * (dx/L)**3 * -# lambda0**2)) - -# ======================================================================== -# 4. velocity fluctuations re-sampled on the desired grid (N^3) -xr = np.linspace(xs, xe, args.N + 1) -xrc = (xr[1:] + xr[:-1]) / 2 -Xr, Yr, Zr = np.meshgrid(xrc, xrc, xrc, indexing="ij") - -Xr = Xr.reshape(-1, order="F") -Yr = Yr.reshape(-1, order="F") -Zr = Zr.reshape(-1, order="F") - -ur = spi.interpn((xc, xc, xc), u, (Xr, Yr, Zr), method="linear") -vr = spi.interpn((xc, xc, xc), v, (Xr, Yr, Zr), method="linear") -wr = spi.interpn((xc, xc, xc), w, (Xr, Yr, Zr), method="linear") - - -# ======================================================================== -# Save the data in Fortran ordering -fname = "hit_ic_{0:d}_{1:d}.dat".format(int(args.k0), args.N) -data = np.vstack((Xr, Yr, Zr, ur, vr, wr)).T -np.savetxt(fname, data, fmt="%.18e", delimiter=",", header="x, y, z, u, v, w") - - -# ======================================================================== -# plot (only u fluctuations) -if args.plot: - import matplotlib as mpl - - mpl.use("Agg") - import matplotlib.pyplot as plt - - datmin = u.min() - datmax = u.max() - # print("min/max u =",datmin,datmax) - - # Original data - # transpose and origin change bc I used meshgrid with ij and not xy - fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(14, 14)) - ax[0, 0].imshow( - uo[:, :, 0].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[0, 0].set_title("Original data (x,y)") - ax[0, 1].imshow( - uo[:, 0, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[0, 1].set_title("Original data (x,z)") - ax[0, 2].imshow( - uo[0, :, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[0, 2].set_title("Original data (y,z)") - - # Filtered original data - ax[1, 0].imshow( - u[:, :, 0].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[1, 0].set_title("Filtered original data (x,y)") - ax[1, 1].imshow( - u[:, 0, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[1, 1].set_title("Filtered original data (x,z)") - ax[1, 2].imshow( - u[0, :, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[1, 2].set_title("Filtered original data (y,z)") - - # Downsampled filtered data - ur = ur.reshape(args.N, args.N, args.N, order="F") - ax[2, 0].imshow( - ur[:, :, 0].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[2, 0].set_title("Downsampled data (x,y)") - ax[2, 1].imshow( - ur[:, 0, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[2, 1].set_title("Downsampled data (x,z)") - ax[2, 2].imshow( - ur[0, :, :].T, - origin="lower", - extent=[xs, xe, xs, xe], - cmap="RdBu_r", - vmin=datmin, - vmax=datmax, - ) - ax[2, 2].set_title("Downsampled data (y,z)") - - plt.savefig("hit_ic_u_{0:d}_{1:d}.png".format(int(args.k0), args.N), format="png") - - # Fourier coefficients of original data - fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(14, 8)) - ax[0, 0].imshow(np.real(uf[:, :, 0].T), origin="lower", cmap="RdBu_r") - ax[0, 0].set_title("Real Fourier coefficients (x,y)") - ax[0, 1].imshow(np.real(uf[:, 0, :].T), origin="lower", cmap="RdBu_r") - ax[0, 1].set_title("Real Fourier coefficients (x,z)") - ax[0, 2].imshow(np.real(uf[0, :, :].T), origin="lower", cmap="RdBu_r") - ax[0, 2].set_title("Real Fourier coefficients (y,z)") - ax[1, 0].imshow(np.imag(uf[:, :, 0].T), origin="lower", cmap="RdBu_r") - ax[1, 0].set_title("Imag Fourier coefficients (x,y)") - ax[1, 1].imshow(np.imag(uf[:, 0, :].T), origin="lower", cmap="RdBu_r") - ax[1, 1].set_title("Imag Fourier coefficients (x,z)") - ax[1, 2].imshow(np.imag(uf[0, :, :].T), origin="lower", cmap="RdBu_r") - ax[1, 2].set_title("Imag Fourier coefficients (y,z)") - plt.savefig("hit_ic_uf_{0:d}_{1:d}.png".format(int(args.k0), args.N), format="png") - - # Spectrum - plt.figure(20) - ax = plt.gca() - p = plt.loglog(khalf, Eko, color=cmap[-1], lw=2) - p[0].set_dashes(dashseq[0]) - p = plt.loglog(ku, Eku, color=cmap[0], lw=2) - p[0].set_dashes(dashseq[1]) - plt.ylim([1e-16, 10]) - plt.xlabel(r"$k$", fontsize=22, fontweight="bold") - plt.ylabel(r"$E(k)$", fontsize=22, fontweight="bold") - plt.setp(ax.get_xmajorticklabels(), fontsize=18, fontweight="bold") - plt.setp(ax.get_ymajorticklabels(), fontsize=18, fontweight="bold") - plt.savefig( - "hit_ic_spectrum_{0:d}_{1:d}.png".format(int(args.k0), args.N), format="png" - ) - -# output timer -end = time.time() - start -print("Elapsed time " + str(timedelta(seconds=end)) + " (or {0:f} seconds)".format(end)) diff --git a/Tutorials/HIT/inputs.3d.forced b/Tutorials/HIT/inputs.3d.forced deleted file mode 100644 index 4ff3a799e..000000000 --- a/Tutorials/HIT/inputs.3d.forced +++ /dev/null @@ -1,129 +0,0 @@ - -#******************************************************************************* - -#NOTE: You may set *either* max_step or stop_time, or you may set them both. - -# Maximum number of coarse grid timesteps to be taken, if stop_time is -# not reached first. -max_step = 2000 -stop_time = 8 -#******************************************************************************* - -# Number of cells in each coordinate direction at the coarsest level -amr.n_cell = 128 128 128 - -#******************************************************************************* - -# Maximum level (defaults to 0 for single level calculation) -amr.max_level = 0 # maximum number of levels of refinement - -#******************************************************************************* - -# Interval (in number of level l timesteps) between regridding -amr.regrid_int = 2 2 2 2 2 2 2 - -#******************************************************************************* - -# Refinement ratio as a function of level -amr.ref_ratio = 2 2 2 2 - -#******************************************************************************* - -# Sets the "NavierStokes" code to be verbose -ns.v = 1 - -#******************************************************************************* - -# Sets the "amr" code to be verbose -amr.v = 1 - -#******************************************************************************* - -# Interval (in number of coarse timesteps) between checkpoint(restart) files -amr.check_int = 100 -amr.check_file = chk - -#******************************************************************************* - -# Interval (in number of coarse timesteps) between plot files -amr.plot_int = 50 -amr.plot_file = plt - -#******************************************************************************* - -# CFL number to be used in calculating the time step : dt = dx / max(velocity) -ns.cfl = 0.7 # CFL number used to set dt - -#******************************************************************************* - -# Factor by which the first time is shrunk relative to CFL constraint -ns.init_shrink = 1.0 # factor which multiplies the very first time step - -#******************************************************************************* - -# Viscosity coefficient -ns.vel_visc_coef = 1.e-4 - -#******************************************************************************* - -# Diffusion coefficient for first scalar -ns.scal_diff_coefs = 0.0 - -#******************************************************************************* - -# Set to 0 if x-y coordinate system, set to 1 if r-z (in 2-d). -geometry.coord_sys = 0 - -#******************************************************************************* - -# Physical dimensions of the low end of the domain. -geometry.prob_lo = -0.5 -0.5 -0.5 - -# Physical dimensions of the high end of the domain. -geometry.prob_hi = 0.5 0.5 0.5 - -#******************************************************************************* - -#Set to 1 if periodic in that direction -geometry.is_periodic = 1 1 1 - -#******************************************************************************* - -# Boundary conditions on the low end of the domain. -ns.lo_bc = 0 0 0 - -# Boundary conditions on the high end of the domain. -ns.hi_bc = 0 0 0 - -# 0 = Interior/Periodic 3 = Symmetry -# 1 = Inflow 4 = SlipWall -# 2 = Outflow 5 = NoSlipWall - -#******************************************************************************* - -# Problem parameters -prob.probtype = 100 - -#******************************************************************************* - -# Turbulent forcing parameters -turb.nmodes = 4 -turb.force_file = forcedata.dat - -# Turn off tiling. Turbulent forcing doesn't work with tiling for now -fabarray.mfiter_tile_size = 1024 1024 1024 - - -#******************************************************************************* - -# Factor by which grids must be coarsenable. -#amr.blocking_factor = 8 - -#******************************************************************************* - -# Add vorticity to the variables in the plot files. -#amr.derive_plot_vars = NONE -amr.derive_plot_vars = mag_vort diveru avg_pressure - -#******************************************************************************* -nodal_proj.proj_tol = 1.e-10 diff --git a/Tutorials/HIT/prob_init.H b/Tutorials/HIT/prob_init.H deleted file mode 100644 index 3c6fc1430..000000000 --- a/Tutorials/HIT/prob_init.H +++ /dev/null @@ -1,41 +0,0 @@ -#ifndef PROB_INIT_H_ -#define PROB_INIT_H_ - - -// This header is included by NavierStokes.H. These are members of NavierStokes - -// -// struct to hold initial conditions parameters -// -struct InitialConditions -{ - // - // For initializing with random combination of cosine waves, used with forced - // turbulence, where ICs are less important as the forcing takes over with time. - // - amrex::Real turb_scale = 1.0; - amrex::Real density = 1.0; -}; - -// -// Problem initialization functions -// -void prob_initData(); - -void init_forced (amrex::Box const& vbx, - /* amrex::Array4 const& press, */ - amrex::Array4 const& vel, - amrex::Array4 const& scal, - const int nscal, - amrex::Box const& domain, - amrex::GpuArray const& dx, - amrex::GpuArray const& problo, - amrex::GpuArray const& probhi, - InitialConditions IC); - -// -// Problems parameters, to be read from inputs file -// -static int probtype; - -#endif diff --git a/Tutorials/HIT/prob_init.cpp b/Tutorials/HIT/prob_init.cpp deleted file mode 100644 index 2b8130ab9..000000000 --- a/Tutorials/HIT/prob_init.cpp +++ /dev/null @@ -1,134 +0,0 @@ -#include -#include -#include - -#ifdef AMREX_USE_TURBULENT_FORCING -#include -#endif - -using namespace amrex; - -int NavierStokes::probtype = -1; - - -// -// Initialize state and pressure with problem-specific data -// -void NavierStokes::prob_initData () -{ - // - // Fill state and, optionally, pressure - // - MultiFab& P_new = get_new_data(Press_Type); - MultiFab& S_new = get_new_data(State_Type); - const int nscal = NUM_STATE-Density; - - S_new.setVal(0.0); - P_new.setVal(0.0); - - // Integer indices of the lower left and upper right corners of the - // valid region of the entire domain. - Box const& domain = geom.Domain(); - auto const& dx = geom.CellSizeArray(); - // Physical coordinates of the lower left corner of the domain - auto const& problo = geom.ProbLoArray(); - auto const& probhi = geom.ProbHiArray(); - -// FIXME - should remove ifdef and use runtime parameter instead... -#ifdef AMREX_USE_TURBULENT_FORCING - // - // Initialize data structures used for homogeneous isentropic forced turbulence. - // Only need to do it once. - if (level == 0) - TurbulentForcing::init_turbulent_forcing(problo,probhi); -#endif - - // - // Create struct to hold initial conditions parameters - // - InitialConditions IC; - - // - // Read problem parameters from inputs file - // - ParmParse pp("prob"); - - pp.query("probtype",probtype); - - if ( probtype == 100 ) - { - // - // Random combination of cosine waves to be used with forced turbulence, - // where ICs are less important as the forcing takes over with time. - // - - pp.query("turb_scale",IC.turb_scale); - pp.query("density_ic",IC.density); - -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(S_new,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& vbx = mfi.tilebox(); - - init_forced(vbx, /*P_new.array(mfi),*/ S_new.array(mfi, Xvel), - S_new.array(mfi, Density), nscal, - domain, dx, problo, probhi, IC); - } - } - else - { - amrex::Abort("NavierStokes::prob_init: unknown probtype"); - } -} - -void NavierStokes::init_forced (Box const& vbx, - /* Array4 const& press, */ - Array4 const& vel, - Array4 const& scal, - const int nscal, - Box const& domain, - GpuArray const& dx, - GpuArray const& problo, - GpuArray const& probhi, - InitialConditions IC) -{ - const auto domlo = amrex::lbound(domain); - - amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real x = problo[0] + (i - domlo.x + 0.5)*dx[0]; - Real y = problo[1] + (j - domlo.y + 0.5)*dx[1]; -#if (AMREX_SPACEDIM == 3) - Real z = problo[2] + (k - domlo.z + 0.5)*dx[2]; -#else - constexpr Real z = 0.0; -#endif - - const Real Lx = (probhi[0] - problo[0]); - const Real Ly = (probhi[1] - problo[1]); -#if (AMREX_SPACEDIM == 3) - const Real Lz = (probhi[2] - problo[2]); -#else - const Real Lz = 1.0; -#endif - // - // Fill Velocity - // - AMREX_D_TERM(vel(i,j,k,0) = IC.turb_scale * std::cos(TwoPi*y/Ly) * std::cos(TwoPi*z/Lz);, - vel(i,j,k,1) = IC.turb_scale * std::cos(TwoPi*x/Lx) * std::cos(TwoPi*z/Lz);, - vel(i,j,k,2) = IC.turb_scale * std::cos(TwoPi*x/Lx) * std::cos(TwoPi*y/Ly);); - - // - // Scalars, ordered as Density, Tracer(s) - // - scal(i,j,k,0) = IC.density; - - // Tracers - for ( int nt=1; nt #include #include +#include #include #include #include @@ -60,8 +61,7 @@ int max_grid_size(4096); const std::string CheckPointVersion = "CheckPointVersion_1.0"; std::string interp_kind; -Real avg_time; -Real avg_time_fluct; +std::string TimeAverageContents; bool TimeAverageFile_exist = false; int flag_eb = 0; @@ -337,10 +337,12 @@ static void ReadCheckpointFile(const std::string& fileName) { if( is_avg.good()) { TimeAverageFile_exist = true; - std::string first_line_avg; - std::getline(is_avg,first_line_avg); - is_avg >> avg_time; - is_avg >> avg_time_fluct; + // Copy the file verbatim. IAMR writes a title line and then one + // (time_avg, time_avg_fluct, dt_avg) triple per level; none of these + // change when the grids are refined or coarsened. + std::ostringstream oss; + oss << is_avg.rdbuf(); + TimeAverageContents = oss.str(); } } @@ -571,9 +573,7 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string amrex::FileOpenFailed(HeaderFileName); } - HeaderFile << "Writing time_average to checkpoint" << '\n' - << avg_time << '\n' - << avg_time_fluct << '\n'; + HeaderFile << TimeAverageContents; } } diff --git a/Util/run_scripts/run.corigpu b/Util/run_scripts/run.corigpu deleted file mode 100755 index 1198390a4..000000000 --- a/Util/run_scripts/run.corigpu +++ /dev/null @@ -1,117 +0,0 @@ -#!/bin/bash -l -#SBATCH -C gpu -#SBATCH -t 00:05:00 -#SBATCH -J AMREX_GPU -#SBATCH -o AMREX_GPU.o%j -#SBATCH -A m3406 -#SBATCH -N 1 -#SBATCH -n 8 -#SBATCH -c 10 -#SBATCH --gres=gpu:8 -#SBATCH --ntasks-per-node=8 - -# Note: Given exclusive configuration mode, -# you MUST specify your desired resources up top like this. -# Cannot put it in the srun line alone. -# (You can force lower than your full request in the srun line, -# or put the configuration again for safety, but shouldn't be needed.) -# ============ -# -N = nodes -# -n = tasks (MPI ranks) -# -c = CPU per task (full coriGPU node, c*n <= 80) -# --gres=gpu: = GPUs per node (full coriGPU node, 8) -# --ntasks-per-node = number of tasks (MPI ranks) per node (full node, 8) -# - -# For one node: -N 1, -n 8, -c 10, --gres=gpu:8 --ntasks-per-node 8 -# For two nodes: -N 2, -n 16, -c 10, --gres=gpu:8 --ntasks-per-node 8 - -# salloc commands: -# (Make sure the appropriate module is loaded to see the gpu partition, esslurm.) -# ================ -# Single GPU. (If you don't require an independent node, please use a shared node.) -# salloc -N 1 -t 2:00:00 -c 10 -C gpu --gres=gpu:1 -A m3406 -# Single node: -# salloc -N 1 -t 2:00:00 -c 80 -C gpu --exclusive --gres=gpu:8 -A m3406 -# Multi node: -# salloc -N 2 -t 2:00:00 -c 80 -C gpu --exclusive --gres=gpu:8 -A m3406 - - -EXE=./main3d.gnu.TPROF.MPI.CUDA.ex -INPUTS=inputs - -# Basic job submissions: -# ============================= -# Run inside the current salloc session using available resources. -# Change parameters to match available resources & run with "./run.corigpu" -# srun -n 8 -c 10 --gres=gpu:8 ${EXE} ${INPUTS} - - -# Submit with the SBATCH configuration above to the gpu queue: "sbatch run.corigpu" -# Can also be ran with "./run.corigpu" to run with 1 CPU and 1 GPU. -srun ${EXE} ${INPUTS} - - - -# NSight Systems -# ============== - -# @@ Simple Example: -#srun nsys profile -o nsys_out.%q{SLURM_PROCID}.%q{SLURM_JOBID} ${EXE} ${INPUTS} - -# @@ Recommended Example: -#srun nsys profile -c nvtx -p "@*" -e NSYS_NVTX_PROFILER_REGISTER_ONLY=0 -o nsys_out.%q{SLURM_PROCID}.%q{SLURM_JOBID} ${EXE} ${INPUTS} - -# @@ Discussion: -# This will run nsys profile and store performance data in a qdrep file named after '-o' -# Open using nsight-sys $(pwd)/nsys_out.#.######.qdrep - -# To capture the NVTX ranges, included in TINY_PROFILE objects, use: -# "-e NSYS_NVTX_PROFILER_REGISTER_ONLY=0" -# (TINY_PROFILE's NVTX regions do not use registered strings at this time.) - -# Nsight systems creates a timeline over a single, contiguous block of time. -# The start of the timeline can be selected using TINY_PROFILER's NVTX markers with: -# -c nvtx -p "region_name@*" -# This will turn on the profiling analysis at the first instance of the TINY_PROFILER region -# and run to the end of the program. To stop the analysis at the end of the same region, add: -# -x true -# Note: This will only analyze the first instance of the region, so "-x true" should be used -# for specific analyses, or on more inclusive timers, e.g. a timer around a full timestep. - -# @@ Documentation: -# For NSight System profiling flags: -# https://docs.nvidia.com/nsight-systems/profiling/index.html#cli-profile-command-switch-options -# For NSight examples to launch profiling, including region limiting: -# https://docs.nvidia.com/nsight-systems/profiling/index.html#example-interactive-cli-command-sequences - -# Running NSight Systems on multiple ranks -# ======================================== - -# Run Nsight Systems only profiling on $PROFILE_RANK rank on a multi-rank job -# **** Preferred for most basic use cases -#srun ./profile_1rank.sh ${EXE} ${INPUTS} - -# Uncomment and copy the following lines into profile_1rank.sh -# Adjust the nsys command line as needed for your test case. -# #!/bin/bash -# PROFILE_RANK=0 -# if [ $SLURM_PROCID == $PROFILE_RANK ]; then -# nsys profile -o nsys_out.%q{SLURM_PROCID}.%q{SLURM_JOBID} "$@" -# else -# "$@" -# fi - -# NSight Compute -# ============== - -# Run Nsight Compute: -# **** This will do a A LOT of analysis. Unless you want the entire job ran 7 times -# **** with full profiling, limit the kernels profiled with additional flags: -# For filtering examples, see: -# https://docs.nvidia.com/nsight-compute/NsightComputeCli/index.html#nvtx-filtering -# For full list of profile options, see: -# https://docs.nvidia.com/nsight-compute/NsightComputeCli/index.html#command-line-options-profile -# Recommended: limit kernels tested within a given BL_PROFILER timer with "--nvtx-include " -# Note: Must use TINY_PROFILE=TRUE and nvtx region names are equal to BL_PROFILER timer names. -#srun nv-nsight-cu-cli -o cucli_out.%q{SLURM_PROCID}.%q{SLURM_JOBID} ${EXE} ${INPUTS} diff --git a/Util/run_scripts/run.perlmutter b/Util/run_scripts/run.perlmutter deleted file mode 100755 index 050029547..000000000 --- a/Util/run_scripts/run.perlmutter +++ /dev/null @@ -1,49 +0,0 @@ -#!/bin/bash - -# SLURM syntax for Cori, probably same for Perlmutter: -# -# -p = partition/queue ("regular", "debug", etc.) -# -N = number of nodes -# --ntasks-per-node = number of tasks per node (if you don't want to pack the node) -# -t = time -# -J = job name -# -o = STDOUT file -# -e = STDERR file (merges with STDOUT if -e is not given explicitly) -# -A = account to charge -# --mail-type = events to e-mail user about (ALL is short for BEGIN,END,FAIL,REQUEUE,STAGE_OUT) -# and also have TIME_LIMIT_50,TIME_LIMIT_90,TIME_LIMIT -# --mail-user = NIM username of user to notify -# -# SLURM by default will cd to your working directory when you submit the job, - -# 1 node, 1 task, 1 GPU -#SBATCH -C gpu -#SBATCH -q regular -#SBATCH -t 1:00:00 -#SBATCH -n 1 -#SBATCH --ntasks-per-node=1 -#SBATCH -c 128 -#SBATCH --gpus-per-task=1 -#SBATCH -J HIT - -EXE=./amr3d.gnu.MPI.CUDA.ex -# For profiling. Also must change srun command -#EXE=./amr3d.gnu.PROF.MPI.CUDA.ex -INPUTS=inputs.3d.HIT -IC_FILE=hit_ic_4_32.dat - -WORKDIR=$SCRATCH/$SLURM_JOB_NAME.$SLURM_JOB_ID -[[ ! -d "$WORKDIR" ]] && mkdir -p "$WORKDIR" - -cp $EXE $INPUTS $IC_FILE $WORKDIR - -cd $WORKDIR - -export SLURM_CPU_BIND="cores" -srun ${EXE} ${INPUTS} - -# For profiling -#srun nsys profile -o nsys_out.fullprof.%q{SLURM_PROCID}.%q{SLURM_JOBID} ${EXE} ${INPUTS} - -# to get interactive node -#salloc --nodes 1 --qos interactive --time 01:00:00 --constraint gpu --gpus 4