mirror/MUSIC
1
0
Fork 0
mirror of https://github.com/cosmo-sims/MUSIC.git synced 2024-09-19 17:03:46 +02:00
Code Releases Activity

merge

Browse source
This commit is contained in:
Oliver Hahn 2013-10-13 00:14:26 +02:00
commit 09c8251d06
8 changed files with 1672 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

26
LICENSE Normal file
View file

@ -0,0 +1,26 @@
MUSIC - Multi-scale Initial Conditions for Cosmological Simulations
Copyright(c) 2011 by Oliver Hahn. All rights reserved.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to use
the Software under the following conditions:
1. Redistributions of the Software must contain the above copyright
notice, this list of conditions and the following disclaimers.
2. Scientific publications that make use of results obtained with
the Software must properly reference that the 'MUSIC' software
has been used and include a reference to Hahn & Abel (2011),
published in MNRAS Vol 415(3), the paper describing the
algorithms used in the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH
THE SOFTWARE.

13
Makefile
View file

@ -3,18 +3,31 @@
FFTW3 = yes
MULTITHREADFFTW = yes
SINGLEPRECISION = no
<<<<<<< local
HAVEHDF5 = yes
=======
HAVEHDF5 = no
>>>>>>> other
HAVEBOXLIB = no
BOXLIB_HOME = ${HOME}/nyx_tot_sterben/BoxLib
##############################################################################
### compiler and path settings
CC = g++
<<<<<<< local
OPT = -Wall -Wno-unknown-pragmas -O3 -g -mtune=native
=======
OPT = -Wall -Wno-unknown-pragmas -O3 -g -msse2
>>>>>>> other
CFLAGS =
LFLAGS = -lgsl -lgslcblas
<<<<<<< local
CPATHS = -I. -I/opt/local/include -I/usr/local/include
LPATHS = -L/opt/local/lib -L/usr/local/lib
=======
CPATHS = -I. -I$(HOME)/local/include -I/opt/local/include -I/usr/local/include
LPATHS = -L$(HOME)/local/lib -L/opt/local/lib -L/usr/local/lib
>>>>>>> other
##############################################################################
# if you have FFTW 2.1.5 or 3.x with multi-thread support, you can enable the

74
ics_example.conf Normal file
View file

@ -0,0 +1,74 @@
[setup]
boxlength = 100
zstart = 50
levelmin = 7
levelmin_TF = 8
levelmax = 9
padding = 8
overlap = 4
ref_center = 0.5, 0.5, 0.5
ref_extent = 0.2, 0.2, 0.2
align_top = no
baryons = no
use_2LPT = no
use_LLA = no
periodic_TF = yes
[cosmology]
Omega_m = 0.276
Omega_L = 0.724
Omega_b = 0.045
H0 = 70.3
sigma_8 = 0.811
nspec = 0.961
transfer = eisenstein
[random]
seed[7] = 12345
seed[8] = 23456
seed[9] = 34567
seed[10] = 45678
seed[11] = 56789
seed[12] = 67890
[output]
##generic MUSIC data format (used for testing)
##requires HDF5 installation and HDF5 enabled in Makefile
#format = generic
#filename = debug.hdf5
##ENZO - also outputs the settings for the parameter file
##requires HDF5 installation and HDF5 enabled in Makefile
#format = enzo
#filename = ic.enzo
##Gadget-2 (type=1: high-res particles, type=5: rest)
#format = gadget2
#filename = ics_gadget.dat
##Grafic2 compatible format for use with RAMSES
##option 'ramses_nml'=yes writes out a startup nml file
#format = grafic2
#filename = ics_ramses
#ramses_nml = yes
##TIPSY compatible with PKDgrav and Gasoline
#format = tipsy
#filename = ics_tipsy.dat
## NYX compatible output format
##requires boxlib installation and boxlib enabled in Makefile
#format = nyx
#filename = init
[poisson]
fft_fine = yes
accuracy = 1e-5
pre_smooth = 3
post_smooth = 3
smoother = gs
laplace_order = 6
grad_order = 6

71
plugins/nyx_plugin/GNUmakefile Normal file
View file

@ -0,0 +1,71 @@
TOP = ${PWD}
PRECISION = DOUBLE
DEBUG = FALSE
DIM = 3
COMP = Intel
FCOMP = Intel
USE_MPI = FALSE
USE_OMP = FALSE
IN_MUSIC = NO
#hmm....awful (needed for general.hh)
ifeq ($(FFTW3), yes)
DEFINES += -DFFTW3
endif
ifeq ($(SINGLE), yes)
PRECISION = SINGLE
endif
include $(BOXLIB_HOME)/Tools/C_mk/Make.defs
NYX = TRUE
VERBOSE = TRUE
DEFINES += -DHAVE_BOXLIB
#These are the directories in Nyx
Bpack += $(TOP)/Make.package
Blocs += $(TOP)
#include $(TOP)/Make.package
include $(Bpack)
INCLUDE_LOCATIONS += $(Blocs)
VPATH_LOCATIONS += $(Blocs)
#These are the directories in BoxLib
Pdirs := C_BaseLib
Ppack += $(foreach dir, $(Pdirs), $(BOXLIB_HOME)/Src/$(dir)/Make.package)
Plocs += $(foreach dir, $(Pdirs), $(BOXLIB_HOME)/Src/$(dir))
include $(Ppack)
INCLUDE_LOCATIONS += $(Plocs)
VPATH_LOCATIONS += $(Plocs)
INCLUDE_LOCATIONS += $(BOXLIB_HOME)/Src/F_BaseLib
VPATH_LOCATIONS += $(BOXLIB_HOME)/Src/F_BaseLib
vpath %.c . $(VPATH_LOCATIONS)
vpath %.cpp . $(VPATH_LOCATIONS)
vpath %.h . $(VPATH_LOCATIONS)
vpath %.H . $(VPATH_LOCATIONS)
vpath %.F . $(VPATH_LOCATIONS)
vpath %.f90 . $(VPATH_LOCATIONS)
vpath %.f . $(VPATH_LOCATIONS)
vpath %.fi . $(VPATH_LOCATIONS)
ifeq ($(IN_MUSIC), NO)
all: $(objForExecs)
@echo BoxLib compiled ...
@touch ../../output.cc
include $(BOXLIB_HOME)/Tools/C_mk/Make.rules
endif

55
plugins/nyx_plugin/Make.ic Normal file
View file

@ -0,0 +1,55 @@
COMP = Intel
FCOMP = Intel
DEBUG = FALSE
include $(BOXLIB_HOME)/Tools/C_mk/Make.defs
NYX = TRUE
VERBOSE = TRUE
DEFINES += -DHAVE_BOXLIB
#These are the directories in Nyx
Bpack += $(TOP)/Make.package
Blocs += $(TOP)
#include $(TOP)/Make.package
include $(Bpack)
INCLUDE_LOCATIONS += $(Blocs)
VPATH_LOCATIONS += $(Blocs)
#These are the directories in BoxLib
Pdirs := C_BaseLib
Ppack += $(foreach dir, $(Pdirs), $(BOXLIB_HOME)/Src/$(dir)/Make.package)
Plocs += $(foreach dir, $(Pdirs), $(BOXLIB_HOME)/Src/$(dir))
include $(Ppack)
INCLUDE_LOCATIONS += $(Plocs)
VPATH_LOCATIONS += $(Plocs)
INCLUDE_LOCATIONS += $(BOXLIB_HOME)/Src/F_BaseLib
VPATH_LOCATIONS += $(BOXLIB_HOME)/Src/F_BaseLib
vpath %.c . $(VPATH_LOCATIONS)
vpath %.cpp . $(VPATH_LOCATIONS)
vpath %.h . $(VPATH_LOCATIONS)
vpath %.H . $(VPATH_LOCATIONS)
vpath %.F . $(VPATH_LOCATIONS)
vpath %.f90 . $(VPATH_LOCATIONS)
vpath %.f . $(VPATH_LOCATIONS)
vpath %.fi . $(VPATH_LOCATIONS)
ifeq ($(IN_MUSIC), NO)
all: $(objForExecs)
@echo BoxLib compiled ...
@touch ../../output.cc
include $(BOXLIB_HOME)/Tools/C_mk/Make.rules
endif

2
plugins/nyx_plugin/Make.package Normal file
View file

@ -0,0 +1,2 @@
CEXE_sources += output_nyx.cpp

654
plugins/nyx_plugin/output_nyx.cpp Normal file
View file

@ -0,0 +1,654 @@
/*
output_nyx.cc - This file is part of MUSIC -
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010 Oliver Hahn
Copyright (C) 2012 Jan Frederik Engels
*/
#ifdef HAVE_BOXLIB
#include "../../output.hh"
#include <VisMF.H>
#include <Box.H>
#include <RealBox.H>
#include <ParallelDescriptor.H>
#include <Utility.H>
#include <PArray.H>
#define MAX_GRID_SIZE 32
#define BL_SPACEDIM 3
class nyx_output_plugin : public output_plugin
{
protected:
struct patch_header{
int component_rank;
size_t component_size;
std::vector<int> dimensions;
int rank;
std::vector<int> top_grid_dims;
std::vector<int> top_grid_end;
std::vector<int> top_grid_start;
};
struct sim_header{
std::vector<int> dimensions;
std::vector<int> offset;
float a_start;
float dx;
float h0;
float omega_b;
float omega_m;
float omega_v;
float vfact;
float boxlength;
int particle_idx;
};
// struct grid_on_one_level{
// IntVect lo;
// IntVect hi;
// };
int n_data_items;
std::vector<std::string> field_name;
int f_lev;
int gridp;
PArray<MultiFab> mfs;
// std::vector<grid_on_one_level> grids;
std::vector<BoxArray> boxarrays;
std::vector<Box> boxes;
sim_header the_sim_header;
void dump_grid_data(int comp, std::string fieldname, const grid_hierarchy& gh, double factor = 1.0, double add = 0.0 )
{
std::cout << fieldname << " is dumped... to mf index " << comp << std::endl;
//FIXME adapt for multiple levels!
for(int mlevel=levelmin_; mlevel<=levelmax_; ++mlevel )
{
int blevel = mlevel-levelmin_;
std::vector<int> ng;
ng.push_back( gh.get_grid(mlevel)->size(0) );
ng.push_back( gh.get_grid(mlevel)->size(1) );
ng.push_back( gh.get_grid(mlevel)->size(2) );
std::cout << ng[0] << " " << ng[1] << " " << ng[2] << std::endl;
//write data to mf
for(MFIter mfi(mfs[blevel]); mfi.isValid(); ++mfi) {
FArrayBox &myFab = mfs[blevel][mfi];
const int *fab_lo = mfi.validbox().loVect();
const int *fab_hi = mfi.validbox().hiVect();
int mk = fab_lo[2] - boxes[blevel].smallEnd()[2];
#ifdef OMP
#pragma omp parallel for default(shared)
#endif
for (int k = fab_lo[2]; k <= fab_hi[2]; k++, mk++) {
int mj = fab_lo[1] - boxes[blevel].smallEnd()[1];
for (int j = fab_lo[1]; j <= fab_hi[1]; j++, mj++) {
int mi = fab_lo[0] - boxes[blevel].smallEnd()[0];
for (int i = fab_lo[0]; i <= fab_hi[0]; i++, mi++) {
if (mi>=ng[0])
std::cout << "mi (" << mi << ") too large " << ng[0] << std::endl;
if (mj>=ng[1])
std::cout << "mj (" << mj << ") too large " << ng[1] << std::endl;
if (mk>=ng[2])
std::cout << "mk (" << mk << ") too large " << ng[2] << std::endl;
IntVect iv(i,j,k);
double data = ( add + (*gh.get_grid(mlevel))(mi,mj,mk) )*factor;
int idx = myFab.box().index(iv);
myFab.dataPtr(comp)[idx] = data;
// mi++;
}
// mj++;
}
// mk++;
}
} // MFI
}
// char nyxname[256], filename[256];
//
// for(unsigned ilevel=levelmin_; ilevel<=levelmax_; ++ilevel )
// {
}
public:
nyx_output_plugin( config_file& cf )
: output_plugin( cf )
{
int argc=1;
char **argv;
BoxLib::Initialize(argc,argv);
bool bhave_hydro = cf_.getValue<bool>("setup","baryons");
if (bhave_hydro)
n_data_items = 10;
else
n_data_items = 6;
field_name.resize(n_data_items);
if (bhave_hydro)
{
field_name[0] = "baryon_density";
field_name[1] = "baryon_vel_x";
field_name[2] = "baryon_vel_y";
field_name[3] = "baryon_vel_z";
field_name[4] = "dm_pos_x";
field_name[5] = "dm_pos_y";
field_name[6] = "dm_pos_z";
field_name[7] = "dm_vel_x";
field_name[8] = "dm_vel_y";
field_name[9] = "dm_vel_z";
the_sim_header.particle_idx = 4;
}
else
{
field_name[0] = "dm_pos_x";
field_name[1] = "dm_pos_y";
field_name[2] = "dm_pos_z";
field_name[3] = "dm_vel_x";
field_name[4] = "dm_vel_y";
field_name[5] = "dm_vel_z";
the_sim_header.particle_idx = 0;
}
f_lev = levelmax_-levelmin_;
std::cout << f_lev+1 << " level" << std::endl;
mfs.resize(f_lev+1);
Array<int> pmap(2);
pmap[0]=0;
pmap[1]=0;
gridp = 1<<levelmin_;
double off[] = {0, 0, 0};
//at first we do this only for the topgrid...
for(int lev = 0; lev <= f_lev; lev++)
{
BoxArray domainBoxArray(1);
int mlev = lev+levelmin_;
int fac = (1<<lev);
// off[0] += fac*offx_[lev];
// off[1] += fac*offy_[lev];
// off[2] += fac*offz_[lev];
off[0] += offx_[lev];
off[1] += offy_[lev];
off[2] += offz_[lev];
for (int asdf = 0; asdf < 3; asdf++)
off[asdf] *= 2;
IntVect pdLo(off[0],
off[1],
off[2]);
IntVect pdHi(off[0]+sizex_[lev]-1,
off[1]+sizey_[lev]-1,
off[2]+sizez_[lev]-1);
// pdLo *= (1<<lev);
// pdHi *= (1<<lev);
std::cout << pdLo << std::endl;
std::cout << pdHi << std::endl;
// Start with a probDomain
// IntVect pdLo(0,0,0);
// IntVect pdHi(gridp-1,gridp-1,gridp-1);
Box probDomain(pdLo,pdHi);
// We just have one box since we don't use mpi.
domainBoxArray.set(0, probDomain);
domainBoxArray.maxSize(32);
pmap.resize(domainBoxArray.size(),0);
DistributionMapping domainDistMap(pmap);
boxarrays.push_back(domainBoxArray);
boxes.push_back(probDomain);
int ngrow(0);
MultiFab *mf = new MultiFab;
mfs.set(lev,mf);
mfs[lev].define(domainBoxArray, n_data_items, ngrow, domainDistMap, Fab_allocate);
}
// if( mkdir( fname_.c_str(), 0777 ) )
// {
// perror( fname_.c_str() );
// throw std::runtime_error("Error in nyx_output_plugin!");
// }
bool haveblockingfactor = cf.containsKey( "setup", "blocking_factor");
if( !haveblockingfactor )
{
LOGERR("nyx output plug-in requires that \'blocking_factor\' is set!");
throw std::runtime_error("nyx output plug-in requires that \'blocking_factor\' is set!");
}
the_sim_header.dimensions.push_back( 1<<levelmin_ );
the_sim_header.dimensions.push_back( 1<<levelmin_ );
the_sim_header.dimensions.push_back( 1<<levelmin_ );
the_sim_header.offset.push_back( 0 );
the_sim_header.offset.push_back( 0 );
the_sim_header.offset.push_back( 0 );
the_sim_header.a_start = 1.0/(1.0+cf.getValue<double>("setup","zstart"));
the_sim_header.dx = cf.getValue<double>("setup","boxlength")/the_sim_header.dimensions[0]/(cf.getValue<double>("cosmology","H0")*0.01); // not sure?!?
the_sim_header.boxlength=cf.getValue<double>("setup","boxlength");
the_sim_header.h0 = cf.getValue<double>("cosmology","H0")*0.01;
if( bhave_hydro )
the_sim_header.omega_b = cf.getValue<double>("cosmology","Omega_b");
else
the_sim_header.omega_b = 0.0;
the_sim_header.omega_m = cf.getValue<double>("cosmology","Omega_m");
the_sim_header.omega_v = cf.getValue<double>("cosmology","Omega_L");
the_sim_header.vfact = cf.getValue<double>("cosmology","vfact")*the_sim_header.h0; //.. need to multiply by h, nyx wants this factor for non h-1 units
std::cout << "creating output object" << std::endl;
}
~nyx_output_plugin()
{
std::string FullPath = fname_;
if (!BoxLib::UtilCreateDirectory(FullPath, 0755))
BoxLib::CreateDirectoryFailed(FullPath);
if (!FullPath.empty() && FullPath[FullPath.size()-1] != '/')
FullPath += '/';
FullPath += "Header";
std::ofstream Header(FullPath.c_str());
for(int lev=0; lev <= f_lev; lev++)
{
writeLevelPlotFile ( fname_,
Header,
VisMF::OneFilePerCPU,
lev);
//FIXME I would prefer VisMF::NFiles
}
Header.close();
writeGridsFile(fname_);
std::cout << "destroying output object" << std::endl;
}
void write_dm_mass( const grid_hierarchy& gh )
{ /* do nothing, not needed */ }
void write_dm_density( const grid_hierarchy& gh )
{ /* write the parameter file data */
// It's very useful to write a parameter file, but WHY here?
}
void write_dm_velocity( int coord, const grid_hierarchy& gh )
{
char nyxname[256];
sprintf( nyxname, "ParticleVelocities_%c", (char)('x'+coord) );
double vunit = 1.0/(1.225e2*sqrt(the_sim_header.omega_m/the_sim_header.a_start));
dump_grid_data(the_sim_header.particle_idx+3+coord, nyxname, gh);
}
void write_dm_position( int coord, const grid_hierarchy& gh )
{
char nyxname[256];
sprintf( nyxname, "ParticleDisplacements_%c", (char)('x'+coord) );
//dump_grid_data( nyxname, gh );
dump_grid_data(the_sim_header.particle_idx+coord, nyxname, gh);
}
void write_dm_potential( const grid_hierarchy& gh )
{ }
void write_gas_potential( const grid_hierarchy& gh )
{ }
void write_gas_velocity( int coord, const grid_hierarchy& gh )
{
double vunit = 1.0/(1.225e2*sqrt(the_sim_header.omega_m/the_sim_header.a_start));
char nyxname[256];
sprintf( nyxname, "GridVelocities_%c", (char)('x'+coord) );
dump_grid_data(coord+1, nyxname, gh);
}
void write_gas_position( int coord, const grid_hierarchy& gh )
{
/* do nothing, not needed */
}
void write_gas_density( const grid_hierarchy& gh )
{
char nyxname[256];
sprintf( nyxname, "density" );
//FIXME factor and add have to be adjusted to the
//corresponding nyx units...
dump_grid_data(0, nyxname, gh);
}
void finalize( void )
{
//
//before finalizing we write out an inputs and a probin file for Nyx.
//
std::ofstream inputs("inputs");
std::ofstream probin("probin");
//at first the fortran stuff...
probin << "&fortin" << std::endl;
probin << " comoving_OmM = " << the_sim_header.omega_m << "d0" << std::endl;
probin << " comoving_OmB = " << the_sim_header.omega_b << "d0" << std::endl;
probin << " comoving_OmL = " << the_sim_header.omega_v << "d0" << std::endl;
probin << " comoving_h = " << the_sim_header.h0 << "d0" << std::endl;
probin << "/" << std::endl;
probin << std::endl;
//afterwards the cpp stuff...(for which we will need a template, which is read in by the code...)
inputs << "nyx.final_a = 1.0 " << std::endl;
inputs << "max_step = 100000 " << std::endl;
inputs << "nyx.small_dens = 1e-4" << std::endl;
inputs << "nyx.small_temp = 10" << std::endl;
inputs << "nyx.cfl = 0.9 # cfl number for hyperbolic system" << std::endl;
inputs << "nyx.init_shrink = 1.0 # scale back initial timestep" << std::endl;
inputs << "nyx.change_max = 1.05 # scale back initial timestep" << std::endl;
inputs << "nyx.dt_cutoff = 5.e-20 # level 0 timestep below which we halt" << std::endl;
inputs << "nyx.sum_interval = 1 # timesteps between computing mass" << std::endl;
inputs << "nyx.v = 1 # verbosity in Castro.cpp" << std::endl;
inputs << "gravity.v = 1 # verbosity in Gravity.cpp" << std::endl;
inputs << "amr.v = 1 # verbosity in Amr.cpp" << std::endl;
inputs << "mg.v = 0 # verbosity in Amr.cpp" << std::endl;
inputs << "particles.v = 1 # verbosity in Particle class" << std::endl;
inputs << "amr.ref_ratio = 2 2 2 2 2 2 2 2 " << std::endl;
inputs << "amr.regrid_int = 2 2 2 2 2 2 2 2 " << std::endl;
inputs << "amr.initial_grid_file = init/grids_file" << std::endl;
inputs << "amr.useFixedCoarseGrids = 1" << std::endl;
inputs << "amr.check_file = chk " << std::endl;
inputs << "amr.check_int = 10 " << std::endl;
inputs << "amr.plot_file = plt " << std::endl;
inputs << "amr.plot_int = 10 " << std::endl;
inputs << "amr.derive_plot_vars = particle_count particle_mass_density pressure" << std::endl;
inputs << "amr.plot_vars = ALL" << std::endl;
inputs << "nyx.add_ext_src = 0" << std::endl;
inputs << "gravity.gravity_type = PoissonGrav " << std::endl;
inputs << "gravity.no_sync = 1 " << std::endl;
inputs << "gravity.no_composite = 1 " << std::endl;
inputs << "mg.bottom_solver = 1 " << std::endl;
inputs << "geometry.is_periodic = 1 1 1 " << std::endl;
inputs << "geometry.coord_sys = 0 " << std::endl;
inputs << "amr.max_grid_size = 32 " << std::endl;
inputs << "nyx.lo_bc = 0 0 0 " << std::endl;
inputs << "nyx.hi_bc = 0 0 0 " << std::endl;
inputs << "nyx.do_grav = 1 " << std::endl;
inputs << "nyx.do_dm_particles = 1 " << std::endl;
inputs << "nyx.particle_init_type = Cosmological " << std::endl;
inputs << "nyx.print_fortran_warnings = 0" << std::endl;
inputs << "cosmo.initDirName = init " << std::endl;
inputs << "nyx.particle_move_type = Gravitational" << std::endl;
inputs << "amr.probin_file = probin " << std::endl;
inputs << "cosmo.ic-source = MUSIC " << std::endl;
inputs << "amr.blocking_factor = " << cf_.getValue<double>("setup","blocking_factor") << std::endl;
inputs << "nyx.do_hydro = "<< (the_sim_header.omega_b>0?1:0) << std::endl;
inputs << "amr.max_level = " << levelmax_-levelmin_ << std::endl;
inputs << "nyx.initial_z = " << 1/the_sim_header.a_start-1 << std::endl;
inputs << "amr.n_cell = " << sizex_[0] << " " << sizey_[0] << " " << sizez_[0] << std::endl;
inputs << "nyx.n_particles = " << sizex_[0] << " " << sizey_[0] << " " << sizez_[0] << std::endl;
inputs << "geometry.prob_lo = 0 0 0" << std::endl;
//double dx = the_sim_header.dx/the_sim_header.h0;
double bl = the_sim_header.boxlength/the_sim_header.h0;
inputs << "geometry.prob_hi = " << bl << " " << bl << " " << bl << std::endl;
probin.close();
inputs.close();
std::cout << "finalizing..." << std::endl;
}
void writeLevelPlotFile (const std::string& dir,
std::ostream& os,
VisMF::How how,
int level)
{
int i, n;
const Real cur_time = 0.0;
std::cout << "in writeLevelPlotFile" << std::endl;
double h0 = cf_.getValue<double>("cosmology", "H0")*0.01;
// for (MFIter mfi(mf); mfi.isValid(); ++mfi)
// {
// std::cout << "bla" << std::endl;
// std::cout << mf[mfi] << std::endl;
// }
if (level == 0)
{
//
// The first thing we write out is the plotfile type.
//
os << "MUSIC_for_Nyx_v0.1" << '\n';
os << n_data_items << '\n';
for (i = 0; i < n_data_items; i++)
os << field_name[i] << '\n';
os << 3 << '\n';
os << 0 << '\n';
os << f_lev << '\n';
for (i = 0; i < BL_SPACEDIM; i++)
os << 0 << ' '; //ProbLo
os << '\n';
double boxlength = cf_.getValue<double>("setup","boxlength");
for (i = 0; i < BL_SPACEDIM; i++)
os << boxlength/h0 << ' '; //ProbHi
os << '\n';
for (i = 0; i < f_lev; i++)
os << 2 << ' '; //refinement factor
os << '\n';
IntVect pdLo(0,0,0);
IntVect pdHi(gridp-1,gridp-1,gridp-1);
// Box probDomain(pdLo,pdHi);
for (i = 0; i <= f_lev; i++) //Geom(i).Domain()
{
// IntVect pdLo(offx_[i], offy_[i], offz_[i]);
// IntVect pdHi(offx_[i]+sizex_[i], offy_[i]+sizey_[i], offz_[i]+sizez_[i]);
Box probDomain(pdLo,pdHi);
os << probDomain << ' ';
pdHi *= 2;
pdHi += 1;
}
os << '\n';
for (i = 0; i <= f_lev; i++) //level steps
os << 0 << ' ';
os << '\n';
double dx = cf_.getValue<double>("setup","boxlength")/gridp/h0;
for (i = 0; i <= f_lev; i++)
{
for (int k = 0; k < BL_SPACEDIM; k++)
os << dx << ' ';
os << '\n';
dx = dx/2.;
}
os << 0 << '\n';
os << "0\n"; // Write bndry data.
}
//
// Build the directory to hold the MultiFab at this level.
// The name is relative to the directory containing the Header file.
//
static const std::string BaseName = "/Cell";
std::string Level = BoxLib::Concatenate("Level_", level, 1);
//
// Now for the full pathname of that directory.
//
std::string FullPath = dir;
if (!FullPath.empty() && FullPath[FullPath.size()-1] != '/')
FullPath += '/';
FullPath += Level;
//
// Only the I/O processor makes the directory if it doesn't already exist.
//
if (!BoxLib::UtilCreateDirectory(FullPath, 0755))
BoxLib::CreateDirectoryFailed(FullPath);
os << level << ' ' << boxarrays[level].size() << ' ' << 0 << '\n';
os << 0 << '\n';
double cellsize[3];
double dx = cf_.getValue<double>("setup","boxlength")/gridp/h0;
for (n = 0; n < BL_SPACEDIM; n++)
{
cellsize[n] = dx;
}
for (i = 0; i < level; i++)
{
for (n = 0; n < BL_SPACEDIM; n++)
{
cellsize[n] /= 2.;
}
}
std::cout << cellsize[0] << std::endl;
for (i = 0; i < boxarrays[level].size(); ++i)
{
double problo[] = {0,0,0};
std::cout << boxarrays[level][i] << std::endl;
RealBox gridloc = RealBox(boxarrays[level][i], cellsize, problo);
for (n = 0; n < BL_SPACEDIM; n++)
os << gridloc.lo(n) << ' ' << gridloc.hi(n) << '\n';
}
//
// The full relative pathname of the MultiFabs at this level.
// The name is relative to the Header file containing this name.
// It's the name that gets written into the Header.
//
std::string PathNameInHeader = Level;
PathNameInHeader += BaseName;
os << PathNameInHeader << '\n';
//
// Use the Full pathname when naming the MultiFab.
//
std::string TheFullPath = FullPath;
TheFullPath += BaseName;
VisMF::Write(mfs[level],TheFullPath,how,true);
}
void writeGridsFile (const std::string& dir)
{
int i, n;
std::cout << "in writeGridsFile" << std::endl;
std::string myFname = dir;
if (!myFname.empty() && myFname[myFname.size()-1] != '/')
myFname += '/';
myFname += "grids_file";
std::ofstream os(myFname.c_str());
os << f_lev << '\n';
for (int lev = 1; lev <= f_lev; lev++)
{
os << boxarrays[lev].size() << '\n';
boxarrays[lev].coarsen(2);
for (i=0; i < boxarrays[lev].size(); i++)
os << boxarrays[lev][i] << "\n";
}
os.close();
}
// void get_grids (const grid_hierarchy &gh)
// {
// for(unsigned ilevel=levelmin_; ilevel<=levelmax_; ++ilevel )
// {
// grid_on_one_level gool;
//
// int xs = gh.get_grid(ilevel)->size(0);
// int ys = gh.get_grid(ilevel)->size(1);
// int zs = gh.get_grid(ilevel)->size(2);
//
// int xo = gh.get_grid(ilevel)->offset(0);
// int yo = gh.get_grid(ilevel)->offset(1);
// int zo = gh.get_grid(ilevel)->offset(2);
//
// IntVect gridLo(xo,yo,zo);
// gool.lo = gridLo;
//
// IntVect gridHi(xo+xs,yo+zs,zo+zs);
// gool.hi = gridHi;
//
// grids.push_back(gool);
// }
//
// }
};
namespace{
output_plugin_creator_concrete<nyx_output_plugin> creator("nyx");
}
#endif //HAVE_BOXLIB

777
tools/ramses/get_music_refmask.f90 Normal file
View file

@ -0,0 +1,777 @@
program get_music_refmask
!--------------------------------------------------------------------------
! Ce programme calcule la carte de densite surfacique projetee
! des particules de matiere noire d'une simulation RAMSES.
! Version F90 par R. Teyssier le 01/04/01.
!--------------------------------------------------------------------------
implicit none
integer::ncpu,ndim,npart,ngrid,n,i,j,k,icpu,ipos,nstar,nstart,inull,ico
integer::ncpu2,npart2,ndim2,levelmin,levelmax,ilevel,ndark,ismooth
integer::nx=0,ny=0,ix,iy,iz,ixp1,iyp1,idim,jdim,ncpu_read
real(KIND=8)::mtot,ddx,ddy,dex,dey,t,soft,poty,mass,btime,unit_l,aexp,unit_t
real(KIND=8)::xmin=0,xmax=1,ymin=0,ymax=1,zmin=0,zmax=1,r,xc=0.5,yc=0.5,zc=0.5,rad=-1
integer::imin,imax,jmin,jmax,kmin,kmax,lmin,ipart
real(KIND=8)::xxmin,xxmax,yymin,yymax,dy,deltax,fakeage
real(KIND=4),dimension(:,:),allocatable::toto
real(KIND=8),dimension(:,:),allocatable::map
real(KIND=8),dimension(:) ,allocatable::x
real(KIND=8),dimension(:) ,allocatable::y
real(KIND=8),dimension(:) ,allocatable::z
real(KIND=8),dimension(:) ,allocatable::vx
real(KIND=8),dimension(:) ,allocatable::vy
real(KIND=8),dimension(:) ,allocatable::vz
real(KIND=8),dimension(:) ,allocatable::m
real(KIND=8),dimension(:) ,allocatable::temp
real(KIND=8),dimension(:) ,allocatable::bt
real(KIND=8),dimension(:) ,allocatable::tempx,tempy,tempz,tempvx,tempvy,tempvz,tempm,tempbt,tempr
integer ,allocatable,dimension(:)::tempid,temp2,indtempid
integer ,allocatable,dimension(:)::id
integer ,allocatable,dimension(:)::idpart,indidpart
integer::outputmode
real ,allocatable,dimension(:,:,:)::imark
character(LEN=1)::proj='z'
character(LEN=5)::nchar,ncharcpu
character(LEN=80)::ordering,format_grille
character(LEN=80)::GMGM
character(LEN=128)::nomfich,repository,filetype='bin',grafic
character(LEN=128)::nomfich2,repository2, outputname
logical::ok,ok_part,periodic=.false.,star=.false.,okerode=.false.
integer::impi,ndom,bit_length,maxdom,maxid,idd
integer,dimension(1:8)::idom,jdom,kdom,cpu_min,cpu_max
real(KIND=8),dimension(1:8)::bounding_min,bounding_max
real(KIND=8)::dkey,order_min,dmax,vfact
real(kind=8),dimension(:),allocatable::bound_key
logical,dimension(:),allocatable::cpu_read
integer,dimension(:),allocatable::cpu_list
integer(kind=4)::np1,np2,np3
real::dx,dx2,x1o,x2o,x3o,astart,omegam,omegav,h0,x1or,x2or,x3or,dxor,omegak
call read_params
!-----------------------------------------------
! Lecture du fichier particules au format RAMSES
!-----------------------------------------------
ipos=INDEX(repository,'output_')
nchar=repository(ipos+7:ipos+13)
nomfich=TRIM(repository)//'/part_'//TRIM(nchar)//'.out00001'
inquire(file=nomfich, exist=ok) ! verify input file
if ( .not. ok ) then
print *,TRIM(nomfich)//' not found.'
stop
endif
nomfich=TRIM(repository)//'/info_'//TRIM(nchar)//'.txt'
inquire(file=nomfich, exist=ok) ! verify input file
if ( .not. ok ) then
print *,TRIM(nomfich)//' not found.'
stop
endif
open(unit=10,file=nomfich,form='formatted',status='old')
read(10,'("ncpu =",I11)')ncpu
read(10,'("ndim =",I11)')ndim
read(10,'("levelmin =",I11)')levelmin
read(10,'("levelmax =",I11)')levelmax
read(10,*)
read(10,*)
read(10,*)
write(*,*)ncpu,ndim,levelmin,levelmax
read(10,*)
read(10,'("time =",E23.15)')t
read(10,'("aexp =",E23.15)')aexp
read(10,*)
read(10,*)
read(10,*)
read(10,*)
read(10,*)
read(10,'("unit_l =",E23.15)')unit_l
read(10,*)
read(10,'("unit_t =",E23.15)')unit_t
read(10,*)
read(10,'("ordering type=",A80)'),ordering
read(10,*)
write(*,'(" ordering type=",A20)'),TRIM(ordering)
allocate(cpu_list(1:ncpu))
if(TRIM(ordering).eq.'hilbert')then
allocate(bound_key(0:ncpu))
allocate(cpu_read(1:ncpu))
cpu_read=.false.
do impi=1,ncpu
read(10,'(I8,1X,E23.15,1X,E23.15)')i,bound_key(impi-1),bound_key(impi)
end do
endif
close(10)
if(TRIM(ordering).eq.'hilbert')then
dmax=max(xmax-xmin,ymax-ymin,zmax-zmin)
do ilevel=1,levelmax
deltax=0.5d0**ilevel
if(deltax.lt.dmax)exit
end do
lmin=ilevel
bit_length=lmin-1
maxdom=2**bit_length
imin=0; imax=0; jmin=0; jmax=0; kmin=0; kmax=0
if(bit_length>0)then
imin=int(xmin*dble(maxdom))
imax=imin+1
jmin=int(ymin*dble(maxdom))
jmax=jmin+1
kmin=int(zmin*dble(maxdom))
kmax=kmin+1
endif
dkey=(dble(2**(levelmax+1)/dble(maxdom)))**ndim
ndom=1
if(bit_length>0)ndom=8
idom(1)=imin; idom(2)=imax
idom(3)=imin; idom(4)=imax
idom(5)=imin; idom(6)=imax
idom(7)=imin; idom(8)=imax
jdom(1)=jmin; jdom(2)=jmin
jdom(3)=jmax; jdom(4)=jmax
jdom(5)=jmin; jdom(6)=jmin
jdom(7)=jmax; jdom(8)=jmax
kdom(1)=kmin; kdom(2)=kmin
kdom(3)=kmin; kdom(4)=kmin
kdom(5)=kmax; kdom(6)=kmax
kdom(7)=kmax; kdom(8)=kmax
do i=1,ndom
if(bit_length>0)then
call hilbert3d(idom(i),jdom(i),kdom(i),order_min,bit_length,1)
else
order_min=0.0d0
endif
bounding_min(i)=(order_min)*dkey
bounding_max(i)=(order_min+1.0D0)*dkey
end do
cpu_min=0; cpu_max=0
do impi=1,ncpu
do i=1,ndom
if ( bound_key(impi-1).le.bounding_min(i).and.&
& bound_key(impi ).gt.bounding_min(i))then
cpu_min(i)=impi
endif
if ( bound_key(impi-1).lt.bounding_max(i).and.&
& bound_key(impi ).ge.bounding_max(i))then
cpu_max(i)=impi
endif
end do
end do
ncpu_read=0
do i=1,ndom
do j=cpu_min(i),cpu_max(i)
if(.not. cpu_read(j))then
ncpu_read=ncpu_read+1
cpu_list(ncpu_read)=j
cpu_read(j)=.true.
endif
enddo
enddo
else
ncpu_read=ncpu
do j=1,ncpu
cpu_list(j)=j
end do
end if
npart=0
do k=1,ncpu_read
write(*,*) 'CPU=',k
icpu=cpu_list(k)
call title(icpu,ncharcpu)
nomfich=TRIM(repository)//'/part_'//TRIM(nchar)//'.out'//TRIM(ncharcpu)
open(unit=1,file=nomfich,status='old',form='unformatted')
read(1)ncpu2
read(1)ndim2
read(1)npart2
read(1)
read(1)nstar
close(1)
npart=npart+npart2
end do
write(*,*) npart,' particles in the region'
allocate(m(1:npart))
allocate(x(1:npart))
allocate(y(1:npart))
allocate(z(1:npart))
allocate(vx(1:npart))
allocate(vy(1:npart))
allocate(vz(1:npart))
allocate(id(1:npart))
if(nstar>0) then
allocate(bt(1:npart))
endif
!-----------------------------------------------
! Compute projected mass using CIC smoothing
!----------------------------------------------
mtot=0.0d0
nstart=1
do k=1,ncpu_read
icpu=cpu_list(k)
call title(icpu,ncharcpu)
nomfich=TRIM(repository)//'/part_'//TRIM(nchar)//'.out'//TRIM(ncharcpu)
open(unit=1,file=nomfich,status='old',form='unformatted')
write(*,*)'Processing file '//TRIM(nomfich)
read(1)ncpu2
read(1)ndim2
read(1)npart2
read(1)
read(1)
read(1)
read(1)
read(1)
allocate(temp(1:npart2))
allocate(temp2(1:npart2))
! Read positions
read(1)temp
x(nstart:nstart+npart2-1)=temp
read(1)temp
y(nstart:nstart+npart2-1)=temp
read(1)temp
z(nstart:nstart+npart2-1)=temp
! Read velocity
read(1)temp
vx(nstart:nstart+npart2-1)=temp
read(1)temp
vy(nstart:nstart+npart2-1)=temp
read(1)temp
vz(nstart:nstart+npart2-1)=temp
! Read mass
read(1)temp
m(nstart:nstart+npart2-1)=temp
!Read identity
read(1)temp2
id(nstart:nstart+npart2-1)=temp2
!Read level
read(1)temp2
if(nstar>0) then
! Read BT
read(1)temp
bt(nstart:nstart+npart2-1)=temp
endif
! ----------------------------
nstart=nstart+npart2 !Fill up the next set
deallocate(temp)
deallocate(temp2)
enddo
40 format(3e16.8)
50 format(2I16)
!Outputs IDs of selected particles
ipart=0
mass=0.0
write(*,*) 'Getting IDs...'
open(18,file='partID.dat',form='formatted')
maxid=0
do i=1,npart !To get maximum identity of the particle
if(nstar.eq.0) then !Only DM particles
btime=0
else
btime=bt(i)
endif
if(btime.eq.0) then
ok_part=(x(i)>=xmin.and.x(i)<=xmax.and. &
& y(i)>=ymin.and.y(i)<=ymax.and. &
& z(i)>=zmin.and.z(i)<=zmax)
if(rad>0) then
r=(x(i)-xc)**2+(y(i)-yc)**2+(z(i)-zc)**2
ok_part=(sqrt(r)<=rad)
endif
if(ok_part) then
maxid=max(maxid,id(i))
ipart=ipart+1
mass=mass+m(i)
endif
endif
enddo
write(*,*) 'We have',ipart,' particles in selected region'
write(*,*) 'Total mass =', mass
30 format(i16)
write(18,50) ipart,npart,maxid
allocate(idpart(1:ipart))
j=1
do i=1,npart !Start finding the IDs
if(nstar.eq.0) then !Only DM particles
btime=0
else
btime=bt(i)
endif
if(btime.eq.0) then
ok_part=(x(i)>=xmin.and.x(i)<=xmax.and. &
& y(i)>=ymin.and.y(i)<=ymax.and. &
& z(i)>=zmin.and.z(i)<=zmax)
if(rad>0) then
r=(x(i)-xc)**2+(y(i)-yc)**2+(z(i)-zc)**2
ok_part=sqrt(r)<=rad
endif
if(ok_part) then
write(18,30) id(i) !Write IDs
idpart(j) = id(i)
j = j+1
endif
endif
enddo
npart = ipart
allocate(indidpart(1:npart))
CALL quick_sort(idpart, indidpart, npart)
!------- read IC data and match -----------
deallocate(m)
deallocate(x)
deallocate(y)
deallocate(z)
deallocate(vx)
deallocate(vy)
deallocate(vz)
deallocate(id)
if(nstar>0) then
deallocate(bt)
endif
allocate(m(1:npart))
allocate(x(1:npart))
allocate(y(1:npart))
allocate(z(1:npart))
allocate(vx(1:npart))
allocate(vy(1:npart))
allocate(vz(1:npart))
allocate(id(1:npart))
allocate(bt(1:npart))
!-----------------------------------------------
! Compute projected mass using CIC smoothing
!----------------------------------------------
mtot=0.0d0
!nstart=1
ipart=0
ipos=INDEX(repository2,'output_')
nchar=repository2(ipos+7:ipos+13)
do k=1,ncpu_read
icpu=cpu_list(k)
call title(icpu,ncharcpu)
nomfich=TRIM(repository2)//'/part_'//TRIM(nchar)//'.out'//TRIM(ncharcpu)
open(unit=1,file=nomfich,status='old',form='unformatted')
write(*,*)'Processing file '//TRIM(nomfich)
read(1)ncpu2
read(1)ndim2
read(1)npart2
read(1)
read(1)
read(1)
read(1)
read(1)
allocate(tempx(1:npart2))
allocate(tempy(1:npart2))
allocate(tempz(1:npart2))
allocate(tempvx(1:npart2))
allocate(tempvy(1:npart2))
allocate(tempvz(1:npart2))
allocate(tempm(1:npart2))
allocate(tempid(1:npart2))
allocate(indtempid(1:npart2))
allocate(temp2(1:npart2))
allocate(tempbt(1:npart2))
! Read positions
read(1)tempx
read(1)tempy
read(1)tempz
! Read velocity
read(1)tempvx
read(1)tempvy
read(1)tempvz
! Read mass
read(1)tempm
! Read identity
read(1)tempid
! Read level
read(1)temp2
if(nstar.gt.0) then
! Read BT
read(1)tempbt
else
tempbt=0
endif
close(1)
call quick_sort(tempid, indtempid, npart2)
ico=1
i=1
do while (i.le.npart2.and.ico.le.npart)
if(tempid(i).lt.idpart(ico))then
i=i+1
else
if(tempid(i).gt.idpart(ico))then
ico=ico+1
else
if(tempid(i).eq.idpart(ico).and.tempbt(i).ne.0)then
i=i+1
else
if(tempid(i).eq.idpart(ico).and.tempbt(i).eq.0)then
ipart=ipart+1
m(ipart)=tempm(indtempid(i))
x(ipart)=tempx(indtempid(i))
y(ipart)=tempy(indtempid(i))
z(ipart)=tempz(indtempid(i))
vx(ipart)=tempvx(indtempid(i))
vy(ipart)=tempvy(indtempid(i))
vz(ipart)=tempvz(indtempid(i))
id(ipart)=tempid(i)
bt(ipart)=tempbt(indtempid(i))
i=i+1
ico=ico+1
end if
end if
end if
end if
end do
deallocate(tempx)
deallocate(tempy)
deallocate(tempz)
deallocate(tempvx)
deallocate(tempvy)
deallocate(tempvz)
deallocate(tempm)
deallocate(tempid)
deallocate(indtempid)
deallocate(temp2)
deallocate(tempbt)
end do
open(20,file=TRIM(outputname),form='formatted')
if( outputmode .eq. 1 ) then
do i=1,ipart
write(20,1001) x(i), y(i), z(i), vx(i), vy(i), vz(i)
end do
else
do i=1,ipart
write(20,1002) x(i), y(i), z(i)
end do
end if
close(20)
write(*,*)'Wrote data to file ',trim(outputname)
1001 format (f16.8,f16.8,f16.8,e16.7,e16.7,e16.7)
1002 format (f16.8,f16.8,f16.8)
contains
subroutine read_params
implicit none
integer :: i,n
integer :: iargc
character(len=4) :: opt
character(len=128) :: arg
LOGICAL :: bad, ok
n = iargc()
if (n < 4) then
print *, 'usage: geticref -inf input_dir_final_snapshot'
print *, ' -ini input_dir_first_snapshot'
print *, ' [-out output_name] '
print *, ' [-xc xc] '
print *, ' [-yc yc] '
print *, ' [-zc zc] '
print *, ' [-rad rad] '
print *, ' [-vel 0|1] output also velocity data'
print *, 'ex: geticref -inf output_00010 -ini output_00001 -xc 0.5 -yc 0.5 -zc 0.5 -rad 0.1'
stop
end if
outputname = 'music_region_file.txt'
outputmode = 0
do i = 1,n,2
call getarg(i,opt)
if (i == n) then
print '("option ",a2," has no argument")', opt
stop 2
end if
call getarg(i+1,arg)
select case (opt)
case ('-inf')
repository = trim(arg)
case ('-ini')
repository2 = trim(arg)
case ('-dir')
proj = trim(arg)
case ('-xmi')
read (arg,*) xmin
case ('-xma')
read (arg,*) xmax
case ('-ymi')
read (arg,*) ymin
case ('-yma')
read (arg,*) ymax
case ('-zmi')
read (arg,*) zmin
case ('-zma')
read (arg,*) zmax
case ('-xc')
read (arg,*) xc
case ('-yc')
read (arg,*) yc
case ('-zc')
read (arg,*) zc
case ('-rad')
read (arg,*) rad
case ('-out')
outputname = trim(arg)
case ('-vel')
read (arg,*) outputmode
case ('-per')
read (arg,*) periodic
case ('-gfc')
grafic = trim(arg)
case ('-fil')
filetype = trim(arg)
case default
print '("unknown option ",a2," ignored")', opt
end select
end do
return
end subroutine read_params
end program get_music_refmask
!=======================================================================
subroutine title(n,nchar)
!=======================================================================
implicit none
integer::n
character*5::nchar
character*1::nchar1
character*2::nchar2
character*3::nchar3
character*4::nchar4
character*5::nchar5
if(n.ge.10000)then
write(nchar5,'(i5)') n
nchar = nchar5
elseif(n.ge.1000)then
write(nchar4,'(i4)') n
nchar = '0'//nchar4
elseif(n.ge.100)then
write(nchar3,'(i3)') n
nchar = '00'//nchar3
elseif(n.ge.10)then
write(nchar2,'(i2)') n
nchar = '000'//nchar2
else
write(nchar1,'(i1)') n
nchar = '0000'//nchar1
endif
end subroutine title
!================================================================
!================================================================
!================================================================
!================================================================
subroutine hilbert3d(x,y,z,order,bit_length,npoint)
implicit none
integer ,INTENT(IN) ::bit_length,npoint
integer ,INTENT(IN) ,dimension(1:npoint)::x,y,z
real(kind=8),INTENT(OUT),dimension(1:npoint)::order
logical,dimension(0:3*bit_length-1)::i_bit_mask
logical,dimension(0:1*bit_length-1)::x_bit_mask,y_bit_mask,z_bit_mask
integer,dimension(0:7,0:1,0:11)::state_diagram
integer::i,ip,cstate,nstate,b0,b1,b2,sdigit,hdigit
if(bit_length>bit_size(bit_length))then
write(*,*)'Maximum bit length=',bit_size(bit_length)
write(*,*)'stop in hilbert3d'
stop
endif
state_diagram = RESHAPE( (/ 1, 2, 3, 2, 4, 5, 3, 5,&
& 0, 1, 3, 2, 7, 6, 4, 5,&
& 2, 6, 0, 7, 8, 8, 0, 7,&
& 0, 7, 1, 6, 3, 4, 2, 5,&
& 0, 9,10, 9, 1, 1,11,11,&
& 0, 3, 7, 4, 1, 2, 6, 5,&
& 6, 0, 6,11, 9, 0, 9, 8,&
& 2, 3, 1, 0, 5, 4, 6, 7,&
& 11,11, 0, 7, 5, 9, 0, 7,&
& 4, 3, 5, 2, 7, 0, 6, 1,&
& 4, 4, 8, 8, 0, 6,10, 6,&
& 6, 5, 1, 2, 7, 4, 0, 3,&
& 5, 7, 5, 3, 1, 1,11,11,&
& 4, 7, 3, 0, 5, 6, 2, 1,&
& 6, 1, 6,10, 9, 4, 9,10,&
& 6, 7, 5, 4, 1, 0, 2, 3,&
& 10, 3, 1, 1,10, 3, 5, 9,&
& 2, 5, 3, 4, 1, 6, 0, 7,&
& 4, 4, 8, 8, 2, 7, 2, 3,&
& 2, 1, 5, 6, 3, 0, 4, 7,&
& 7, 2,11, 2, 7, 5, 8, 5,&
& 4, 5, 7, 6, 3, 2, 0, 1,&
& 10, 3, 2, 6,10, 3, 4, 4,&
& 6, 1, 7, 0, 5, 2, 4, 3 /), &
& (/8 ,2, 12 /) )
do ip=1,npoint
! convert to binary
do i=0,bit_length-1
x_bit_mask(i)=btest(x(ip),i)
y_bit_mask(i)=btest(y(ip),i)
z_bit_mask(i)=btest(z(ip),i)
enddo
! interleave bits
do i=0,bit_length-1
i_bit_mask(3*i+2)=x_bit_mask(i)
i_bit_mask(3*i+1)=y_bit_mask(i)
i_bit_mask(3*i )=z_bit_mask(i)
end do
! build Hilbert ordering using state diagram
cstate=0
do i=bit_length-1,0,-1
b2=0 ; if(i_bit_mask(3*i+2))b2=1
b1=0 ; if(i_bit_mask(3*i+1))b1=1
b0=0 ; if(i_bit_mask(3*i ))b0=1
sdigit=b2*4+b1*2+b0
nstate=state_diagram(sdigit,0,cstate)
hdigit=state_diagram(sdigit,1,cstate)
i_bit_mask(3*i+2)=btest(hdigit,2)
i_bit_mask(3*i+1)=btest(hdigit,1)
i_bit_mask(3*i )=btest(hdigit,0)
cstate=nstate
enddo
! save Hilbert key as double precision real
order(ip)=0.
do i=0,3*bit_length-1
b0=0 ; if(i_bit_mask(i))b0=1
order(ip)=order(ip)+dble(b0)*dble(2)**i
end do
end do
end subroutine hilbert3d
SUBROUTINE quick_sort(list,order,n)
! Quick sort routine from:
! Brainerd, W.S., Goldberg, C.H. & Adams, J.C. (1990) "Programmer's Guide to
! Fortran 90", McGraw-Hill ISBN 0-07-000248-7, pages 149-150.
! Modified by Alan Miller to include an associated integer array which gives
! the positions of the elements in the original order.
IMPLICIT NONE
INTEGER :: n
INTEGER, DIMENSION (1:n), INTENT(INOUT) :: list
INTEGER, DIMENSION (1:n), INTENT(OUT) :: order
! Local variable
INTEGER :: i
DO i = 1, n
order(i) = i
END DO
CALL quick_sort_1(1, n)
CONTAINS
RECURSIVE SUBROUTINE quick_sort_1(left_end, right_end)
INTEGER, INTENT(IN) :: left_end, right_end
! Local variables
INTEGER :: i, j, itemp
INTEGER :: reference, temp
INTEGER, PARAMETER :: max_simple_sort_size = 6
IF (right_end < left_end + max_simple_sort_size) THEN
! Use interchange sort for small lists
CALL interchange_sort(left_end, right_end)
ELSE
! Use partition ("quick") sort
reference = list((left_end + right_end)/2)
i = left_end - 1; j = right_end + 1
DO
! Scan list from left end until element >= reference is found
DO
i = i + 1
IF (list(i) >= reference) EXIT
END DO
! Scan list from right end until element <= reference is found
DO
j = j - 1
IF (list(j) <= reference) EXIT
END DO
IF (i < j) THEN
! Swap two out-of-order elements
temp = list(i); list(i) = list(j); list(j) = temp
itemp = order(i); order(i) = order(j); order(j) = itemp
ELSE IF (i == j) THEN
i = i + 1
EXIT
ELSE
EXIT
END IF
END DO
IF (left_end < j) CALL quick_sort_1(left_end, j)
IF (i < right_end) CALL quick_sort_1(i, right_end)
END IF
END SUBROUTINE quick_sort_1
SUBROUTINE interchange_sort(left_end, right_end)
INTEGER, INTENT(IN) :: left_end, right_end
! Local variables
INTEGER :: i, j, itemp
INTEGER :: temp
DO i = left_end, right_end - 1
DO j = i+1, right_end
IF (list(i) > list(j)) THEN
temp = list(i); list(i) = list(j); list(j) = temp
itemp = order(i); order(i) = order(j); order(j) = itemp
END IF
END DO
END DO
END SUBROUTINE interchange_sort
END SUBROUTINE quick_sort