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

wip commit

Browse source
This commit is contained in:
Oliver Hahn 2024-02-26 05:18:53 +08:00
parent c3b80ce0c7
commit a616b2ebf6
3 changed files with 82 additions and 35 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

34
src/densities.cc
View file

@ -17,6 +17,8 @@
#include <cstring>
#include "math/special.hh"
#include "densities.hh"
#include "random.hh"
#include "convolution_kernel.hh"
@ -24,20 +26,6 @@
//TODO: this should be a larger number by default, just to maintain consistency with old default
#define DEF_RAN_CUBE_SIZE 32
double Meyer_scaling_function( double k, double kmax )
{
constexpr double twopithirds{2.0*M_PI/3.0};
constexpr double fourpithirds{4.0*M_PI/3.0};
auto nu = []( double x ){ return x<0.0?0.0:(x<1.0?x:1.0); };
k = std::abs(k)/kmax * fourpithirds;
if( k < twopithirds ) return 1.0;
else if( k< fourpithirds ){
return std::cos( 0.5*M_PI * nu(3*k/(2*M_PI)-1.0) );
}
return 0.0;
}
template <typename m1, typename m2>
void fft_coarsen(m1 &v, m2 &V)
@ -131,6 +119,7 @@ void fft_interpolate(m1 &V, m2 &v, int margin, bool from_basegrid = false)
size_t nxf = v.size(0), nyf = v.size(1), nzf = v.size(2), nzfp = nzf + 2;
size_t mxf = v.margin(0), myf = v.margin(1), mzf = v.margin(2);
// adjust offsets to respect margins, all grids have 'margins' except basegrid (which is periodic)
if (!from_basegrid)
{
oxf += mxf/2;
@ -217,22 +206,21 @@ void fft_interpolate(m1 &V, m2 &v, int margin, bool from_basegrid = false)
double kz = (k <= (int)nzc / 2) ? (double)k : (double)(k - (int)nzc);
double phase = -0.5 * M_PI * (kx / nxc + ky / nyc + kz / nzc);
// double phase = -0.5 * M_PI * (kx / nxc + ky / nyc + kz / nzc);
std::complex<double> val_phas(cos(phase), sin(phase));
std::complex<double> val(RE(ccoarse[qc]), IM(ccoarse[qc]));
val *= val_phas * 8.0;
if(i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){
double blend_coarse_x = Meyer_scaling_function(kx, nxc / 2);
double blend_coarse_y = Meyer_scaling_function(ky, nyc / 2);
double blend_coarse_z = Meyer_scaling_function(kz, nzc / 2);
double blend_coarse = blend_coarse_x*blend_coarse_y*blend_coarse_z;
double blend_fine = 1.0-blend_coarse;
double blend_coarse_x = Meyer_scaling_function(kx, nxc / 4);
double blend_coarse_y = Meyer_scaling_function(ky, nyc / 4);
double blend_coarse_z = Meyer_scaling_function(kz, nzc / 4);
double blend_coarse = blend_coarse_x*blend_coarse_y*blend_coarse_z;
double blend_fine = 1.0-blend_coarse;
RE(cfine[qf]) = blend_fine * RE(cfine[qf]) + blend_coarse * val.real();
IM(cfine[qf]) = blend_fine * IM(cfine[qf]) + blend_coarse * val.imag();
}
RE(cfine[qf]) = blend_fine * RE(cfine[qf]) + blend_coarse * val.real();
IM(cfine[qf]) = blend_fine * IM(cfine[qf]) + blend_coarse * val.imag();
}
delete[] rcoarse;

42
src/math/special.hh Normal file
View file

@ -0,0 +1,42 @@
// This file is part of monofonIC (MUSIC2)
// A software package to generate ICs for cosmological simulations
// Copyright (C) 2020 by Oliver Hahn
//
// monofonIC is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// monofonIC is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#pragma once
#include <cmath>
inline double Meyer_scaling_function( double k, double kmax )
{
constexpr double twopithirds{2.0*M_PI/3.0};
constexpr double fourpithirds{4.0*M_PI/3.0};
auto nu = []( double x ){ return x<0.0?0.0:(x<1.0?x:1.0); };
// k = std::abs(k)/kmax * fourpithirds;
k = std::abs(k)/kmax * 2 * M_PI;
if( k < twopithirds ) return 1.0;
else if( k< fourpithirds ){
return std::cos( 0.5*M_PI * nu(3*k/(2*M_PI)-1.0) );
}
// if( k < fourpithirds ) return 1.0;
return 0.0;
}
inline double Shannon_scaling_function( double k, double kmax )
{
if( std::abs(k) < kmax ) return 1.0;
return 0.0;
}

41
src/plugins/random_music_wnoise_generator.cc
View file

@ -4,6 +4,8 @@
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "math/special.hh"
#include "random.hh"
#include "random_music_wnoise_generator.hh"
@ -627,18 +629,33 @@ music_wnoise_generator<T>::music_wnoise_generator(music_wnoise_generator<T> &rc,
val *= val_phas * sqrt8;
// if (x0_ == NULL || lx_ == NULL){
if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2)
{
RE(cfine[qf]) = val.real();
IM(cfine[qf]) = val.imag();
}
else
{
RE(cfine[qf]) = val.real();
IM(cfine[qf]) = 0.0;
}
// }
// // if (x0_ == NULL || lx_ == NULL){
// if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2)
// {
// RE(cfine[qf]) = val.real();
// IM(cfine[qf]) = val.imag();
// }
// else
// {
// RE(cfine[qf]) = val.real();
// IM(cfine[qf]) = 0.0;
// }
// // }
if(i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){
double blend_coarse_x = Meyer_scaling_function(kx, nxc / 2);
double blend_coarse_y = Meyer_scaling_function(ky, nyc / 2);
double blend_coarse_z = Meyer_scaling_function(kz, nzc / 2);
// double blend_coarse_x = Shannon_scaling_function(kx, nxc / 2);
// double blend_coarse_y = Shannon_scaling_function(ky, nyc / 2);
// double blend_coarse_z = Shannon_scaling_function(kz, nzc / 2);
double blend_coarse = blend_coarse_x*blend_coarse_y*blend_coarse_z;
double blend_fine = std::sqrt(1.0-blend_coarse*blend_coarse);
RE(cfine[qf]) = blend_fine * RE(cfine[qf]) + blend_coarse * val.real();
IM(cfine[qf]) = blend_fine * IM(cfine[qf]) + blend_coarse * val.imag();
}
}
delete[] rcoarse;