2020-08-21 17:03:33 +02:00
|
|
|
// 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/>.
|
2019-11-03 15:54:17 +01:00
|
|
|
#pragma once
|
|
|
|
|
|
|
|
#include <general.hh>
|
|
|
|
#include <unistd.h> // for unlink
|
|
|
|
|
|
|
|
#include <iostream>
|
|
|
|
#include <fstream>
|
|
|
|
|
|
|
|
#include <random>
|
2019-12-05 05:43:49 +01:00
|
|
|
#include <map>
|
2019-11-03 15:54:17 +01:00
|
|
|
|
2020-02-28 16:15:37 +01:00
|
|
|
#include <cassert>
|
|
|
|
|
2019-11-06 17:55:09 +01:00
|
|
|
#include <particle_generator.hh>
|
2019-11-04 00:25:45 +01:00
|
|
|
#include <grid_fft.hh>
|
2020-04-04 23:59:13 +02:00
|
|
|
#include <math/mat3.hh>
|
2019-11-03 15:54:17 +01:00
|
|
|
|
2020-02-28 16:15:37 +01:00
|
|
|
#include <gsl/gsl_sf_hyperg.h>
|
|
|
|
inline double Hypergeometric2F1( double a, double b, double c, double x )
|
|
|
|
{
|
|
|
|
return gsl_sf_hyperg_2F1( a, b, c, x);
|
|
|
|
}
|
|
|
|
|
2019-12-04 14:26:42 +01:00
|
|
|
#define PRODUCTION
|
2019-11-14 15:36:39 +01:00
|
|
|
|
2019-11-03 15:54:17 +01:00
|
|
|
namespace particle{
|
2019-12-05 05:43:49 +01:00
|
|
|
//! implement Joyce, Marcos et al. PLT calculation
|
2019-11-03 15:54:17 +01:00
|
|
|
|
2019-11-05 19:14:14 +01:00
|
|
|
class lattice_gradient{
|
|
|
|
private:
|
2020-03-01 17:21:17 +01:00
|
|
|
const real_t boxlen_, aini_;
|
2019-11-05 19:14:14 +01:00
|
|
|
const size_t ngmapto_, ngrid_, ngrid32_;
|
2020-03-01 17:21:17 +01:00
|
|
|
const real_t mapratio_, XmL_;
|
2019-11-14 15:36:39 +01:00
|
|
|
Grid_FFT<real_t,false> D_xx_, D_xy_, D_xz_, D_yy_, D_yz_, D_zz_;
|
|
|
|
Grid_FFT<real_t,false> grad_x_, grad_y_, grad_z_;
|
2020-03-29 14:45:43 +02:00
|
|
|
std::vector<vec3_t<real_t>> vectk_;
|
|
|
|
std::vector<vec3_t<int>> ico_, vecitk_;
|
2019-11-05 19:14:14 +01:00
|
|
|
|
2019-12-04 14:26:42 +01:00
|
|
|
bool is_even( int i ){ return (i%2)==0; }
|
|
|
|
|
2020-03-29 14:45:43 +02:00
|
|
|
bool is_in( int i, int j, int k, const mat3_t<int>& M ){
|
|
|
|
vec3_t<int> v({i,j,k});
|
2019-12-04 14:26:42 +01:00
|
|
|
auto vv = M * v;
|
|
|
|
return is_even(vv.x)&&is_even(vv.y)&&is_even(vv.z);
|
|
|
|
}
|
|
|
|
|
2019-12-02 01:04:03 +01:00
|
|
|
void init_D( lattice lattice_type )
|
2019-11-05 19:14:14 +01:00
|
|
|
{
|
2019-12-01 15:40:24 +01:00
|
|
|
constexpr real_t pi = M_PI;
|
|
|
|
constexpr real_t twopi = 2.0*M_PI;
|
2019-12-01 15:04:25 +01:00
|
|
|
constexpr real_t fourpi = 4.0*M_PI;
|
2019-12-01 15:40:24 +01:00
|
|
|
const real_t sqrtpi = std::sqrt(M_PI);
|
|
|
|
const real_t pi32 = std::pow(M_PI,1.5);
|
2019-11-27 16:23:43 +01:00
|
|
|
|
2019-12-01 15:04:25 +01:00
|
|
|
//! === vectors, reciprocals and normals for the SC lattice ===
|
|
|
|
const int charge_fac_sc = 1;
|
2020-03-29 14:45:43 +02:00
|
|
|
const mat3_t<real_t> mat_bravais_sc{
|
2019-12-01 15:04:25 +01:00
|
|
|
1.0, 0.0, 0.0,
|
|
|
|
0.0, 1.0, 0.0,
|
|
|
|
0.0, 0.0, 1.0,
|
|
|
|
};
|
2020-03-29 14:45:43 +02:00
|
|
|
const mat3_t<real_t> mat_reciprocal_sc{
|
2019-12-01 15:04:25 +01:00
|
|
|
twopi, 0.0, 0.0,
|
|
|
|
0.0, twopi, 0.0,
|
|
|
|
0.0, 0.0, twopi,
|
|
|
|
};
|
2020-03-29 14:45:43 +02:00
|
|
|
const mat3_t<int> mat_invrecip_sc{
|
2019-12-04 14:26:42 +01:00
|
|
|
2, 0, 0,
|
|
|
|
0, 2, 0,
|
|
|
|
0, 0, 2,
|
|
|
|
};
|
2020-03-29 14:45:43 +02:00
|
|
|
const std::vector<vec3_t<real_t>> normals_sc{
|
2019-12-01 15:04:25 +01:00
|
|
|
{pi,0.,0.},{-pi,0.,0.},
|
|
|
|
{0.,pi,0.},{0.,-pi,0.},
|
|
|
|
{0.,0.,pi},{0.,0.,-pi},
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
//! === vectors, reciprocals and normals for the BCC lattice ===
|
|
|
|
const int charge_fac_bcc = 2;
|
2020-03-29 14:45:43 +02:00
|
|
|
const mat3_t<real_t> mat_bravais_bcc{
|
2019-12-01 14:34:28 +01:00
|
|
|
1.0, 0.0, 0.5,
|
|
|
|
0.0, 1.0, 0.5,
|
|
|
|
0.0, 0.0, 0.5,
|
2019-11-27 16:23:43 +01:00
|
|
|
};
|
2020-03-29 14:45:43 +02:00
|
|
|
const mat3_t<real_t> mat_reciprocal_bcc{
|
2019-12-01 14:34:28 +01:00
|
|
|
twopi, 0.0, 0.0,
|
|
|
|
0.0, twopi, 0.0,
|
2019-12-01 15:04:25 +01:00
|
|
|
-twopi, -twopi, fourpi,
|
2019-11-27 16:23:43 +01:00
|
|
|
};
|
2020-03-29 14:45:43 +02:00
|
|
|
const mat3_t<int> mat_invrecip_bcc{
|
2019-12-04 14:26:42 +01:00
|
|
|
2, 0, 0,
|
|
|
|
0, 2, 0,
|
|
|
|
1, 1, 1,
|
|
|
|
};
|
2020-03-29 14:45:43 +02:00
|
|
|
const std::vector<vec3_t<real_t>> normals_bcc{
|
2019-11-27 16:23:43 +01:00
|
|
|
{0.,pi,pi},{0.,-pi,pi},{0.,pi,-pi},{0.,-pi,-pi},
|
|
|
|
{pi,0.,pi},{-pi,0.,pi},{pi,0.,-pi},{-pi,0.,-pi},
|
|
|
|
{pi,pi,0.},{-pi,pi,0.},{pi,-pi,0.},{-pi,-pi,0.}
|
|
|
|
};
|
2019-12-01 15:04:25 +01:00
|
|
|
|
2019-11-27 16:23:43 +01:00
|
|
|
|
2019-12-01 15:04:25 +01:00
|
|
|
//! === vectors, reciprocals and normals for the FCC lattice ===
|
|
|
|
const int charge_fac_fcc = 4;
|
2020-03-29 14:45:43 +02:00
|
|
|
const mat3_t<real_t> mat_bravais_fcc{
|
2019-12-01 15:04:25 +01:00
|
|
|
0.0, 0.5, 0.0,
|
|
|
|
0.5, 0.0, 1.0,
|
|
|
|
0.5, 0.5, 0.0,
|
|
|
|
};
|
2020-03-29 14:45:43 +02:00
|
|
|
const mat3_t<real_t> mat_reciprocal_fcc{
|
2019-12-01 15:04:25 +01:00
|
|
|
-fourpi, fourpi, twopi,
|
|
|
|
0.0, 0.0, twopi,
|
|
|
|
fourpi, 0.0, -twopi,
|
|
|
|
};
|
2020-03-29 14:45:43 +02:00
|
|
|
const mat3_t<int> mat_invrecip_fcc{
|
2019-12-04 14:26:42 +01:00
|
|
|
0, 1, 1,
|
|
|
|
1, 0, 1,
|
|
|
|
0, 2, 0,
|
|
|
|
};
|
2020-03-29 14:45:43 +02:00
|
|
|
const std::vector<vec3_t<real_t>> normals_fcc{
|
2019-12-01 15:04:25 +01:00
|
|
|
{twopi,0.,0.},{-twopi,0.,0.},
|
|
|
|
{0.,twopi,0.},{0.,-twopi,0.},
|
|
|
|
{0.,0.,twopi},{0.,0.,-twopi},
|
|
|
|
{+pi,+pi,+pi},{+pi,+pi,-pi},
|
|
|
|
{+pi,-pi,+pi},{+pi,-pi,-pi},
|
|
|
|
{-pi,+pi,+pi},{-pi,+pi,-pi},
|
|
|
|
{-pi,-pi,+pi},{-pi,-pi,-pi},
|
|
|
|
};
|
2019-12-01 14:34:28 +01:00
|
|
|
|
2019-12-01 15:04:25 +01:00
|
|
|
//! select the properties for the chosen lattice
|
2019-12-02 01:04:03 +01:00
|
|
|
const int ilat = lattice_type; // 0 = sc, 1 = bcc, 2 = fcc
|
2019-12-01 15:04:25 +01:00
|
|
|
|
2019-12-02 01:04:03 +01:00
|
|
|
const auto mat_bravais = (ilat==2)? mat_bravais_fcc : (ilat==1)? mat_bravais_bcc : mat_bravais_sc;
|
|
|
|
const auto mat_reciprocal = (ilat==2)? mat_reciprocal_fcc : (ilat==1)? mat_reciprocal_bcc : mat_reciprocal_sc;
|
2019-12-04 14:26:42 +01:00
|
|
|
const auto mat_invrecip = (ilat==2)? mat_invrecip_fcc : (ilat==1)? mat_invrecip_bcc : mat_invrecip_sc;
|
2019-12-02 01:04:03 +01:00
|
|
|
const auto normals = (ilat==2)? normals_fcc : (ilat==1)? normals_bcc : normals_sc;
|
|
|
|
const auto charge_fac = (ilat==2)? charge_fac_fcc : (ilat==1)? charge_fac_bcc : charge_fac_sc;
|
2019-12-01 14:34:28 +01:00
|
|
|
|
2019-12-01 18:52:53 +01:00
|
|
|
const ptrdiff_t nlattice = ngrid_;
|
2019-12-01 14:34:28 +01:00
|
|
|
const real_t dx = 1.0/real_t(nlattice);
|
|
|
|
|
2019-12-02 01:04:03 +01:00
|
|
|
const real_t eta = 4.0; // Ewald cutoff shall be 4 cells
|
2019-12-01 14:34:28 +01:00
|
|
|
const real_t alpha = 1.0/std::sqrt(2)/eta;
|
|
|
|
const real_t alpha2 = alpha*alpha;
|
|
|
|
const real_t alpha3 = alpha2*alpha;
|
|
|
|
|
2019-12-01 15:04:25 +01:00
|
|
|
const real_t charge = 1.0/std::pow(real_t(nlattice),3)/charge_fac;
|
2019-12-01 14:34:28 +01:00
|
|
|
const real_t fft_norm = 1.0/std::pow(real_t(nlattice),3.0);
|
2019-12-01 11:28:17 +01:00
|
|
|
const real_t fft_norm12 = 1.0/std::pow(real_t(nlattice),1.5);
|
2019-12-01 14:34:28 +01:00
|
|
|
|
|
|
|
//! just a Kronecker \delta_ij
|
|
|
|
auto kronecker = []( int i, int j ) -> real_t { return (i==j)? 1.0 : 0.0; };
|
|
|
|
|
2019-12-01 15:40:24 +01:00
|
|
|
//! Ewald summation: short-range Green's function
|
2020-03-29 14:45:43 +02:00
|
|
|
auto add_greensftide_sr = [&]( mat3_t<real_t>& D, const vec3_t<real_t>& d ) -> void {
|
2019-12-01 14:34:28 +01:00
|
|
|
auto r = d.norm();
|
|
|
|
if( r< 1e-14 ) return; // return zero for r=0
|
|
|
|
|
|
|
|
const real_t r2(r*r), r3(r2*r), r5(r3*r2);
|
|
|
|
const real_t K1( -alpha3/pi32 * std::exp(-alpha2*r2)/r2 );
|
|
|
|
const real_t K2( (std::erfc(alpha*r) + 2.0*alpha/sqrtpi*std::exp(-alpha2*r2)*r)/fourpi );
|
|
|
|
|
|
|
|
for( int mu=0; mu<3; ++mu ){
|
|
|
|
for( int nu=mu; nu<3; ++nu ){
|
|
|
|
real_t dd( d[mu]*d[nu] * K1 + (kronecker(mu,nu)/r3 - 3.0 * (d[mu]*d[nu])/r5) * K2 );
|
|
|
|
D(mu,nu) += dd;
|
|
|
|
D(nu,mu) += (mu!=nu)? dd : 0.0;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
};
|
|
|
|
|
2019-12-01 15:40:24 +01:00
|
|
|
//! Ewald summation: long-range Green's function
|
2020-03-29 14:45:43 +02:00
|
|
|
auto add_greensftide_lr = [&]( mat3_t<real_t>& D, const vec3_t<real_t>& k, const vec3_t<real_t>& r ) -> void {
|
2019-12-01 14:34:28 +01:00
|
|
|
real_t kmod2 = k.norm_squared();
|
|
|
|
real_t term = std::exp(-kmod2/(4*alpha2))*std::cos(k.dot(r)) / kmod2 * fft_norm;
|
|
|
|
for( int mu=0; mu<3; ++mu ){
|
|
|
|
for( int nu=mu; nu<3; ++nu ){
|
|
|
|
auto dd = k[mu] * k[nu] * term;
|
|
|
|
D(mu,nu) += dd;
|
|
|
|
D(nu,mu) += (mu!=nu)? dd : 0.0;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
};
|
2019-12-01 15:40:24 +01:00
|
|
|
|
|
|
|
//! checks if 'vec' is in the FBZ with FBZ normal vectors given in 'normals'
|
|
|
|
auto check_FBZ = []( const auto& normals, const auto& vec ) -> bool {
|
|
|
|
for( const auto& n : normals ){
|
2019-12-01 18:52:53 +01:00
|
|
|
if( n.dot( vec ) > 1.0001 * n.dot(n) ){
|
2019-12-19 13:54:24 +01:00
|
|
|
return false;
|
2019-12-01 15:40:24 +01:00
|
|
|
}
|
|
|
|
}
|
2019-12-19 13:54:24 +01:00
|
|
|
return true;
|
2019-12-01 15:40:24 +01:00
|
|
|
};
|
2019-11-27 16:23:43 +01:00
|
|
|
|
2019-12-01 15:40:24 +01:00
|
|
|
constexpr ptrdiff_t lnumber = 3, knumber = 3;
|
2019-12-02 01:04:03 +01:00
|
|
|
const int numb = 1; //!< search radius when shifting vectors into FBZ
|
2019-12-01 11:28:17 +01:00
|
|
|
|
2020-03-29 14:45:43 +02:00
|
|
|
vectk_.assign(D_xx_.memsize(),vec3_t<real_t>());
|
|
|
|
ico_.assign(D_xx_.memsize(),vec3_t<int>());
|
|
|
|
vecitk_.assign(D_xx_.memsize(),vec3_t<int>());
|
2019-12-01 11:28:17 +01:00
|
|
|
|
2019-12-01 14:34:28 +01:00
|
|
|
#pragma omp parallel
|
|
|
|
{
|
|
|
|
//... temporary to hold values of the dynamical matrix
|
2020-03-29 14:45:43 +02:00
|
|
|
mat3_t<real_t> matD(0.0);
|
2019-12-01 11:28:17 +01:00
|
|
|
|
2019-12-01 14:34:28 +01:00
|
|
|
#pragma omp for
|
2019-12-02 01:04:03 +01:00
|
|
|
for( ptrdiff_t i=0; i<nlattice; ++i ){
|
|
|
|
for( ptrdiff_t j=0; j<nlattice; ++j ){
|
|
|
|
for( ptrdiff_t k=0; k<nlattice; ++k ){
|
2019-12-01 14:34:28 +01:00
|
|
|
// compute lattice site vector from (i,j,k) multiplying Bravais base matrix, and wrap back to box
|
2020-03-29 14:45:43 +02:00
|
|
|
const vec3_t<real_t> x_ijk({dx*real_t(i),dx*real_t(j),dx*real_t(k)});
|
|
|
|
const vec3_t<real_t> ar = (mat_bravais * x_ijk).wrap_abs();
|
2019-12-01 14:34:28 +01:00
|
|
|
|
|
|
|
//... zero temporary matrix
|
|
|
|
matD.zero();
|
|
|
|
|
|
|
|
// add real-space part of dynamical matrix, periodic copies
|
|
|
|
for( ptrdiff_t ix=-lnumber; ix<=lnumber; ix++ ){
|
|
|
|
for( ptrdiff_t iy=-lnumber; iy<=lnumber; iy++ ){
|
|
|
|
for( ptrdiff_t iz=-lnumber; iz<=lnumber; iz++ ){
|
2020-03-29 14:45:43 +02:00
|
|
|
const vec3_t<real_t> n_ijk({real_t(ix),real_t(iy),real_t(iz)});
|
|
|
|
const vec3_t<real_t> dr(ar - mat_bravais * n_ijk);
|
2019-12-01 14:34:28 +01:00
|
|
|
add_greensftide_sr(matD, dr);
|
2019-12-01 11:28:17 +01:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2019-12-01 14:34:28 +01:00
|
|
|
|
|
|
|
// add k-space part of dynamical matrix
|
|
|
|
for( ptrdiff_t ix=-knumber; ix<=knumber; ix++ ){
|
|
|
|
for( ptrdiff_t iy=-knumber; iy<=knumber; iy++ ){
|
|
|
|
for( ptrdiff_t iz=-knumber; iz<=knumber; iz++ ){
|
|
|
|
if(std::abs(ix)+std::abs(iy)+std::abs(iz) != 0){
|
2020-03-29 14:45:43 +02:00
|
|
|
const vec3_t<real_t> k_ijk({real_t(ix)/nlattice,real_t(iy)/nlattice,real_t(iz)/nlattice});
|
|
|
|
const vec3_t<real_t> ak( mat_reciprocal * k_ijk);
|
2019-12-01 14:34:28 +01:00
|
|
|
|
|
|
|
add_greensftide_lr(matD, ak, ar );
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
D_xx_.relem(i,j,k) = matD(0,0) * charge;
|
|
|
|
D_xy_.relem(i,j,k) = matD(0,1) * charge;
|
|
|
|
D_xz_.relem(i,j,k) = matD(0,2) * charge;
|
|
|
|
D_yy_.relem(i,j,k) = matD(1,1) * charge;
|
|
|
|
D_yz_.relem(i,j,k) = matD(1,2) * charge;
|
|
|
|
D_zz_.relem(i,j,k) = matD(2,2) * charge;
|
|
|
|
}
|
2019-11-27 16:23:43 +01:00
|
|
|
}
|
|
|
|
}
|
2019-12-01 14:34:28 +01:00
|
|
|
} // end omp parallel region
|
2019-11-27 16:23:43 +01:00
|
|
|
|
2019-12-01 11:28:17 +01:00
|
|
|
// fix r=0 with background density (added later in Fourier space)
|
|
|
|
D_xx_.relem(0,0,0) = 1.0/3.0;
|
|
|
|
D_xy_.relem(0,0,0) = 0.0;
|
|
|
|
D_xz_.relem(0,0,0) = 0.0;
|
|
|
|
D_yy_.relem(0,0,0) = 1.0/3.0;
|
|
|
|
D_yz_.relem(0,0,0) = 0.0;
|
|
|
|
D_zz_.relem(0,0,0) = 1.0/3.0;
|
2019-11-27 16:23:43 +01:00
|
|
|
|
2019-12-01 11:28:17 +01:00
|
|
|
D_xx_.FourierTransformForward();
|
|
|
|
D_xy_.FourierTransformForward();
|
|
|
|
D_xz_.FourierTransformForward();
|
|
|
|
D_yy_.FourierTransformForward();
|
|
|
|
D_yz_.FourierTransformForward();
|
|
|
|
D_zz_.FourierTransformForward();
|
|
|
|
|
2019-12-02 01:04:03 +01:00
|
|
|
#ifndef PRODUCTION
|
2019-12-01 11:28:17 +01:00
|
|
|
if (CONFIG::MPI_task_rank == 0)
|
|
|
|
unlink("debug.hdf5");
|
|
|
|
D_xx_.Write_to_HDF5("debug.hdf5","Dxx");
|
|
|
|
D_xy_.Write_to_HDF5("debug.hdf5","Dxy");
|
|
|
|
D_xz_.Write_to_HDF5("debug.hdf5","Dxz");
|
|
|
|
D_yy_.Write_to_HDF5("debug.hdf5","Dyy");
|
|
|
|
D_yz_.Write_to_HDF5("debug.hdf5","Dyz");
|
|
|
|
D_zz_.Write_to_HDF5("debug.hdf5","Dzz");
|
|
|
|
|
|
|
|
std::ofstream ofs2("test_brillouin.txt");
|
2019-12-02 01:04:03 +01:00
|
|
|
#endif
|
2020-03-29 14:45:43 +02:00
|
|
|
using map_t = std::map<vec3_t<int>,size_t>;
|
2019-12-05 05:43:49 +01:00
|
|
|
map_t iimap;
|
2019-12-02 01:04:03 +01:00
|
|
|
|
2019-12-05 05:43:49 +01:00
|
|
|
//!=== Make temporary copies before resorting to std. Fourier grid ========!//
|
|
|
|
Grid_FFT<real_t,false>
|
|
|
|
temp1({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}),
|
|
|
|
temp2({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}),
|
|
|
|
temp3({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0});
|
|
|
|
|
|
|
|
temp1.FourierTransformForward(false);
|
|
|
|
temp2.FourierTransformForward(false);
|
|
|
|
temp3.FourierTransformForward(false);
|
|
|
|
|
|
|
|
#pragma omp parallel for
|
|
|
|
for( size_t i=0; i<D_xx_.size(0); i++ )
|
|
|
|
{
|
|
|
|
for( size_t j=0; j<D_xx_.size(1); j++ )
|
2019-11-27 16:23:43 +01:00
|
|
|
{
|
2019-12-05 05:43:49 +01:00
|
|
|
for( size_t k=0; k<D_xx_.size(2); k++ )
|
2019-12-01 11:28:17 +01:00
|
|
|
{
|
2019-12-05 05:43:49 +01:00
|
|
|
temp1.kelem(i,j,k) = ccomplex_t(std::real(D_xx_.kelem(i,j,k)),std::real(D_xy_.kelem(i,j,k)));
|
|
|
|
temp2.kelem(i,j,k) = ccomplex_t(std::real(D_xz_.kelem(i,j,k)),std::real(D_yy_.kelem(i,j,k)));
|
|
|
|
temp3.kelem(i,j,k) = ccomplex_t(std::real(D_yz_.kelem(i,j,k)),std::real(D_zz_.kelem(i,j,k)));
|
2019-12-02 01:04:03 +01:00
|
|
|
}
|
|
|
|
}
|
2019-12-05 05:43:49 +01:00
|
|
|
}
|
|
|
|
D_xx_.zero(); D_xy_.zero(); D_xz_.zero();
|
|
|
|
D_yy_.zero(); D_yz_.zero(); D_zz_.zero();
|
|
|
|
|
|
|
|
|
|
|
|
//!=== Diagonalise and resort to std. Fourier grid ========!//
|
|
|
|
#pragma omp parallel
|
|
|
|
{
|
|
|
|
// thread private matrix representation
|
2020-03-29 14:45:43 +02:00
|
|
|
mat3_t<real_t> D;
|
|
|
|
vec3_t<real_t> eval, evec1, evec2, evec3_t;
|
2019-12-01 11:28:17 +01:00
|
|
|
|
2019-12-05 05:43:49 +01:00
|
|
|
#pragma omp for
|
|
|
|
for( size_t i=0; i<D_xx_.size(0); i++ )
|
|
|
|
{
|
|
|
|
for( size_t j=0; j<D_xx_.size(1); j++ )
|
2019-12-02 01:04:03 +01:00
|
|
|
{
|
2019-12-05 05:43:49 +01:00
|
|
|
for( size_t k=0; k<D_xx_.size(2); k++ )
|
2019-12-02 01:04:03 +01:00
|
|
|
{
|
2020-03-29 14:45:43 +02:00
|
|
|
vec3_t<real_t> kv = D_xx_.get_k<real_t>(i,j,k);
|
2019-12-05 05:43:49 +01:00
|
|
|
|
|
|
|
// put matrix elements into actual matrix
|
|
|
|
D(0,0) = std::real(temp1.kelem(i,j,k)) / fft_norm12;
|
|
|
|
D(0,1) = D(1,0) = std::imag(temp1.kelem(i,j,k)) / fft_norm12;
|
|
|
|
D(0,2) = D(2,0) = std::real(temp2.kelem(i,j,k)) / fft_norm12;
|
|
|
|
D(1,1) = std::imag(temp2.kelem(i,j,k)) / fft_norm12;
|
|
|
|
D(1,2) = D(2,1) = std::real(temp3.kelem(i,j,k)) / fft_norm12;
|
|
|
|
D(2,2) = std::imag(temp3.kelem(i,j,k)) / fft_norm12;
|
|
|
|
|
|
|
|
// compute eigenstructure of matrix
|
2020-03-29 14:45:43 +02:00
|
|
|
D.eigen(eval, evec1, evec2, evec3_t);
|
|
|
|
evec3_t /= (twopi*ngrid_);
|
2019-12-05 05:43:49 +01:00
|
|
|
|
|
|
|
// now determine to which modes on the regular lattice this contributes
|
2020-03-29 14:45:43 +02:00
|
|
|
vec3_t<real_t> ar = kv / (twopi*ngrid_);
|
|
|
|
vec3_t<real_t> a(mat_reciprocal * ar);
|
2019-12-05 05:43:49 +01:00
|
|
|
|
|
|
|
// translate the k-vectors into the "candidate" FBZ
|
|
|
|
for( int l1=-numb; l1<=numb; ++l1 ){
|
|
|
|
for( int l2=-numb; l2<=numb; ++l2 ){
|
|
|
|
for( int l3=-numb; l3<=numb; ++l3 ){
|
|
|
|
// need both halfs of Fourier space since we use real transforms
|
|
|
|
for( int isign=0; isign<=1; ++isign ){
|
2019-12-19 13:54:24 +01:00
|
|
|
const real_t sign = 2.0*real_t(isign)-1.0;
|
2020-03-29 14:45:43 +02:00
|
|
|
const vec3_t<real_t> vshift({real_t(l1),real_t(l2),real_t(l3)});
|
2019-12-02 01:04:03 +01:00
|
|
|
|
2020-03-29 14:45:43 +02:00
|
|
|
vec3_t<real_t> vectk = sign * a + mat_reciprocal * vshift;
|
2019-12-02 01:04:03 +01:00
|
|
|
|
|
|
|
if( check_FBZ( normals, vectk ) )
|
|
|
|
{
|
|
|
|
int ix = std::round(vectk.x*(ngrid_)/twopi);
|
|
|
|
int iy = std::round(vectk.y*(ngrid_)/twopi);
|
|
|
|
int iz = std::round(vectk.z*(ngrid_)/twopi);
|
|
|
|
|
2019-12-05 05:43:49 +01:00
|
|
|
#pragma omp critical
|
2020-03-29 14:45:43 +02:00
|
|
|
{iimap.insert( std::pair<vec3_t<int>,size_t>({ix,iy,iz}, D_xx_.get_idx(i,j,k)) );}
|
2019-12-05 05:43:49 +01:00
|
|
|
|
|
|
|
temp1.kelem(i,j,k) = ccomplex_t(eval[2],eval[1]);
|
2020-03-29 14:45:43 +02:00
|
|
|
temp2.kelem(i,j,k) = ccomplex_t(eval[0],evec3_t.x);
|
|
|
|
temp3.kelem(i,j,k) = ccomplex_t(evec3_t.y,evec3_t.z);
|
2019-12-01 15:40:24 +01:00
|
|
|
}
|
2019-12-05 05:43:49 +01:00
|
|
|
}//sign
|
|
|
|
} //l3
|
|
|
|
} //l2
|
|
|
|
} //l1
|
|
|
|
} //k
|
|
|
|
} //j
|
|
|
|
} //i
|
|
|
|
}
|
2019-12-01 18:52:53 +01:00
|
|
|
|
2019-12-05 05:43:49 +01:00
|
|
|
D_xx_.kelem(0,0,0) = 1.0;
|
|
|
|
D_xy_.kelem(0,0,0) = 0.0;
|
|
|
|
D_xz_.kelem(0,0,0) = 0.0;
|
2019-12-04 14:26:42 +01:00
|
|
|
|
2019-12-05 05:43:49 +01:00
|
|
|
D_yy_.kelem(0,0,0) = 1.0;
|
|
|
|
D_yz_.kelem(0,0,0) = 0.0;
|
|
|
|
D_zz_.kelem(0,0,0) = 0.0;
|
2019-12-01 18:52:53 +01:00
|
|
|
|
2019-12-02 01:04:03 +01:00
|
|
|
//... approximate infinite lattice by inerpolating to sites not convered by current resolution...
|
2019-12-05 05:43:49 +01:00
|
|
|
#pragma omp parallel for
|
|
|
|
for( size_t i=0; i<D_xx_.size(0); i++ ){
|
|
|
|
for( size_t j=0; j<D_xx_.size(1); j++ ){
|
|
|
|
for( size_t k=0; k<D_xx_.size(2); k++ ){
|
|
|
|
int ii = (int(i)>nlattice/2)? int(i)-nlattice : int(i);
|
|
|
|
int jj = (int(j)>nlattice/2)? int(j)-nlattice : int(j);
|
|
|
|
int kk = (int(k)>nlattice/2)? int(k)-nlattice : int(k);
|
2020-03-29 14:45:43 +02:00
|
|
|
vec3_t<real_t> kv({real_t(ii),real_t(jj),real_t(kk)});
|
2019-12-05 05:43:49 +01:00
|
|
|
|
2020-03-29 14:45:43 +02:00
|
|
|
auto align_with_k = [&]( const vec3_t<real_t>& v ) -> vec3_t<real_t>{
|
2019-12-05 05:43:49 +01:00
|
|
|
return v*((v.dot(kv)<0.0)?-1.0:1.0);
|
|
|
|
};
|
|
|
|
|
2020-03-29 14:45:43 +02:00
|
|
|
vec3_t<real_t> v, l;
|
2019-12-05 05:43:49 +01:00
|
|
|
map_t::iterator it;
|
|
|
|
|
|
|
|
if( !is_in(i,j,k,mat_invrecip) ){
|
2020-03-29 14:45:43 +02:00
|
|
|
auto average_lv = [&]( const auto& t1, const auto& t2, const auto& t3, vec3_t<real_t>& v, vec3_t<real_t>& l ) {
|
2019-12-05 05:43:49 +01:00
|
|
|
v = 0.0; l = 0.0;
|
|
|
|
int count(0);
|
|
|
|
|
|
|
|
auto add_lv = [&]( auto it ) -> void {
|
|
|
|
auto q = it->second;++count;
|
2020-03-29 14:45:43 +02:00
|
|
|
l += vec3_t<real_t>({std::real(t1.kelem(q)),std::imag(t1.kelem(q)),std::real(t2.kelem(q))});
|
|
|
|
v += align_with_k(vec3_t<real_t>({std::imag(t2.kelem(q)),std::real(t3.kelem(q)),std::imag(t3.kelem(q))}));
|
2019-12-02 01:04:03 +01:00
|
|
|
};
|
2019-12-05 05:43:49 +01:00
|
|
|
map_t::iterator it;
|
|
|
|
if( (it = iimap.find({ii-1,jj,kk}))!=iimap.end() ){ add_lv(it); }
|
|
|
|
if( (it = iimap.find({ii+1,jj,kk}))!=iimap.end() ){ add_lv(it); }
|
|
|
|
if( (it = iimap.find({ii,jj-1,kk}))!=iimap.end() ){ add_lv(it); }
|
|
|
|
if( (it = iimap.find({ii,jj+1,kk}))!=iimap.end() ){ add_lv(it); }
|
|
|
|
if( (it = iimap.find({ii,jj,kk-1}))!=iimap.end() ){ add_lv(it); }
|
|
|
|
if( (it = iimap.find({ii,jj,kk+1}))!=iimap.end() ){ add_lv(it); }
|
|
|
|
l/=real_t(count); v/=real_t(count);
|
|
|
|
};
|
2019-12-01 20:10:58 +01:00
|
|
|
|
2019-12-05 05:43:49 +01:00
|
|
|
average_lv(temp1,temp2,temp3,v,l);
|
|
|
|
|
|
|
|
}else{
|
|
|
|
if( (it = iimap.find({ii,jj,kk}))!=iimap.end() ){
|
|
|
|
auto q = it->second;
|
2020-03-29 14:45:43 +02:00
|
|
|
l = vec3_t<real_t>({std::real(temp1.kelem(q)),std::imag(temp1.kelem(q)),std::real(temp2.kelem(q))});
|
|
|
|
v = align_with_k(vec3_t<real_t>({std::imag(temp2.kelem(q)),std::real(temp3.kelem(q)),std::imag(temp3.kelem(q))}));
|
2019-12-01 20:10:58 +01:00
|
|
|
}
|
|
|
|
}
|
2019-12-05 05:43:49 +01:00
|
|
|
D_xx_.kelem(i,j,k) = l[0];
|
|
|
|
D_xy_.kelem(i,j,k) = l[1];
|
|
|
|
D_xz_.kelem(i,j,k) = l[2];
|
|
|
|
D_yy_.kelem(i,j,k) = v[0];
|
|
|
|
D_yz_.kelem(i,j,k) = v[1];
|
|
|
|
D_zz_.kelem(i,j,k) = v[2];
|
2019-12-01 20:10:58 +01:00
|
|
|
}
|
|
|
|
}
|
2019-11-27 16:23:43 +01:00
|
|
|
}
|
2019-12-01 20:10:58 +01:00
|
|
|
|
2019-12-02 01:04:03 +01:00
|
|
|
#ifdef PRODUCTION
|
|
|
|
#pragma omp parallel for
|
|
|
|
for( size_t i=0; i<D_xx_.size(0); i++ ){
|
|
|
|
for( size_t j=0; j<D_xx_.size(1); j++ ){
|
|
|
|
for( size_t k=0; k<D_xx_.size(2); k++ )
|
|
|
|
{
|
2020-03-29 14:45:43 +02:00
|
|
|
vec3_t<real_t> kv = D_xx_.get_k<real_t>(i,j,k);
|
2019-12-02 01:04:03 +01:00
|
|
|
|
|
|
|
double mu1 = std::real(D_xx_.kelem(i,j,k));
|
2019-12-04 14:26:42 +01:00
|
|
|
// double mu2 = std::real(D_xy_.kelem(i,j,k));
|
|
|
|
// double mu3 = std::real(D_xz_.kelem(i,j,k));
|
|
|
|
|
2020-03-29 14:45:43 +02:00
|
|
|
vec3_t<real_t> evec1({std::real(D_yy_.kelem(i,j,k)),std::real(D_yz_.kelem(i,j,k)),std::real(D_zz_.kelem(i,j,k))});
|
2019-12-04 14:26:42 +01:00
|
|
|
evec1 /= evec1.norm();
|
|
|
|
|
2020-03-01 17:21:17 +01:00
|
|
|
// ///////////////////////////////////
|
|
|
|
// // project onto spherical coordinate vectors
|
|
|
|
|
|
|
|
real_t kr = kv.norm(), kphi = kr>0.0? std::atan2(kv.y,kv.x) : 0.0, ktheta = kr>0.0? std::acos( kv.z / kr ): 0.0;
|
|
|
|
real_t st = std::sin(ktheta), ct = std::cos(ktheta), sp = std::sin(kphi), cp = std::cos(kphi);
|
2020-03-29 14:45:43 +02:00
|
|
|
vec3_t<real_t> e_r( st*cp, st*sp, ct), e_theta( ct*cp, ct*sp, -st), e_phi( -sp, cp, 0.0 );
|
2019-12-04 14:26:42 +01:00
|
|
|
|
2020-03-01 17:21:17 +01:00
|
|
|
// re-normalise to that longitudinal amplitude is exact
|
|
|
|
double renorm = evec1.dot( e_r ); if( renorm < 0.01 ) renorm = 1.0;
|
2020-02-28 16:15:37 +01:00
|
|
|
|
2020-03-01 17:21:17 +01:00
|
|
|
// -- store in diagonal components of D_ij
|
|
|
|
D_xx_.kelem(i,j,k) = 1.0;
|
|
|
|
D_yy_.kelem(i,j,k) = evec1.dot( e_theta ) / renorm;
|
|
|
|
D_zz_.kelem(i,j,k) = evec1.dot( e_phi ) / renorm;
|
2019-12-04 14:26:42 +01:00
|
|
|
|
2020-03-01 17:21:17 +01:00
|
|
|
// spatially dependent correction to vfact = \dot{D_+}/D_+
|
|
|
|
D_xy_.kelem(i,j,k) = 1.0/(0.25*(std::sqrt(1.+24*mu1)-1.));
|
2019-12-02 01:04:03 +01:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
D_xy_.kelem(0,0,0) = 1.0;
|
2020-03-01 17:21:17 +01:00
|
|
|
D_xx_.kelem(0,0,0) = 1.0;
|
2019-12-04 14:26:42 +01:00
|
|
|
D_yy_.kelem(0,0,0) = 0.0;
|
|
|
|
D_zz_.kelem(0,0,0) = 0.0;
|
|
|
|
|
|
|
|
// unlink("debug.hdf5");
|
|
|
|
// D_xy_.Write_to_HDF5("debug.hdf5","mu1");
|
|
|
|
// D_xx_.Write_to_HDF5("debug.hdf5","e1x");
|
|
|
|
// D_yy_.Write_to_HDF5("debug.hdf5","e1y");
|
|
|
|
// D_zz_.Write_to_HDF5("debug.hdf5","e1z");
|
|
|
|
|
2019-12-02 01:04:03 +01:00
|
|
|
#else
|
|
|
|
D_xx_.Write_to_HDF5("debug.hdf5","mu1");
|
2019-12-01 11:28:17 +01:00
|
|
|
D_xy_.Write_to_HDF5("debug.hdf5","mu2");
|
|
|
|
D_xz_.Write_to_HDF5("debug.hdf5","mu3");
|
|
|
|
D_yy_.Write_to_HDF5("debug.hdf5","e1x");
|
|
|
|
D_yz_.Write_to_HDF5("debug.hdf5","e1y");
|
|
|
|
D_zz_.Write_to_HDF5("debug.hdf5","e1z");
|
2019-12-04 14:26:42 +01:00
|
|
|
#endif
|
2019-11-27 16:23:43 +01:00
|
|
|
}
|
2019-12-01 11:28:17 +01:00
|
|
|
|
2019-11-14 14:09:24 +01:00
|
|
|
|
2019-11-05 19:14:14 +01:00
|
|
|
public:
|
2019-11-06 17:55:09 +01:00
|
|
|
// real_t boxlen, size_t ngridother
|
2020-04-04 20:55:24 +02:00
|
|
|
explicit lattice_gradient( config_file& the_config, size_t ngridself=64 )
|
|
|
|
: boxlen_( the_config.get_value<double>("setup", "BoxLength") ),
|
|
|
|
aini_ ( 1.0/(1.0+the_config.get_value<double>("setup", "zstart")) ),
|
|
|
|
ngmapto_( the_config.get_value<size_t>("setup", "GridRes") ),
|
2019-11-06 17:55:09 +01:00
|
|
|
ngrid_( ngridself ), ngrid32_( std::pow(ngrid_, 1.5) ), mapratio_(real_t(ngrid_)/real_t(ngmapto_)),
|
2020-04-04 20:55:24 +02:00
|
|
|
XmL_ ( the_config.get_value<double>("cosmology", "Omega_L") / the_config.get_value<double>("cosmology", "Omega_m") ),
|
2019-11-05 19:14:14 +01:00
|
|
|
D_xx_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_xy_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}),
|
|
|
|
D_xz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_yy_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}),
|
|
|
|
D_yz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_zz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}),
|
|
|
|
grad_x_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), grad_y_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}),
|
2019-12-02 01:04:03 +01:00
|
|
|
grad_z_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0})
|
2019-11-05 19:14:14 +01:00
|
|
|
{
|
2020-04-04 20:27:51 +02:00
|
|
|
music::ilog << "-------------------------------------------------------------------------------" << std::endl;
|
2020-04-04 20:55:24 +02:00
|
|
|
std::string lattice_str = the_config.get_value_safe<std::string>("setup","ParticleLoad","sc");
|
2019-11-06 17:55:09 +01:00
|
|
|
const lattice lattice_type =
|
|
|
|
((lattice_str=="bcc")? lattice_bcc
|
|
|
|
: ((lattice_str=="fcc")? lattice_fcc
|
|
|
|
: ((lattice_str=="rsc")? lattice_rsc
|
|
|
|
: lattice_sc)));
|
|
|
|
|
2020-04-04 20:27:51 +02:00
|
|
|
music::ilog << "PLT corrections for " << lattice_str << " lattice will be computed on " << ngrid_ << "**3 mesh" << std::endl;
|
2019-11-06 17:55:09 +01:00
|
|
|
|
2019-11-05 19:14:14 +01:00
|
|
|
double wtime = get_wtime();
|
2020-04-04 20:27:51 +02:00
|
|
|
music::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing PLT eigenmodes "<< std::flush;
|
2019-11-05 19:14:14 +01:00
|
|
|
|
2019-12-02 01:04:03 +01:00
|
|
|
init_D( lattice_type );
|
2019-12-01 18:52:53 +01:00
|
|
|
// init_D__old();
|
2019-11-05 19:14:14 +01:00
|
|
|
|
2020-04-04 20:27:51 +02:00
|
|
|
music::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
2019-11-05 19:14:14 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
inline ccomplex_t gradient( const int idim, std::array<size_t,3> ijk ) const
|
|
|
|
{
|
|
|
|
real_t ix = ijk[0]*mapratio_, iy = ijk[1]*mapratio_, iz = ijk[2]*mapratio_;
|
2020-02-28 16:15:37 +01:00
|
|
|
|
2020-03-01 17:21:17 +01:00
|
|
|
auto kv = D_xx_.get_k<real_t>( ix, iy, iz );
|
|
|
|
auto kmod = kv.norm() / mapratio_ / boxlen_;
|
2020-02-28 16:15:37 +01:00
|
|
|
|
2020-02-28 18:35:26 +01:00
|
|
|
// // project onto spherical coordinate vectors
|
2020-03-01 17:21:17 +01:00
|
|
|
auto D_r = std::real(D_xx_.get_cic_kspace({ix,iy,iz}));
|
|
|
|
auto D_theta = std::real(D_yy_.get_cic_kspace({ix,iy,iz}));
|
|
|
|
auto D_phi = std::real(D_zz_.get_cic_kspace({ix,iy,iz}));
|
2020-02-28 16:15:37 +01:00
|
|
|
|
2020-03-01 17:21:17 +01:00
|
|
|
real_t kr = kv.norm(), kphi = kr>0.0? std::atan2(kv.y,kv.x) : 0.0, ktheta = kr>0.0? std::acos( kv.z / kr ) : 0.0;
|
|
|
|
real_t st = std::sin(ktheta), ct = std::cos(ktheta), sp = std::sin(kphi), cp = std::cos(kphi);
|
2020-02-28 18:35:26 +01:00
|
|
|
|
2020-03-01 17:21:17 +01:00
|
|
|
if( idim == 0 ){
|
|
|
|
return ccomplex_t(0.0, kmod*(D_r * st * cp + D_theta * ct * cp - D_phi * sp));
|
|
|
|
}
|
|
|
|
else if( idim == 1 ){
|
|
|
|
return ccomplex_t(0.0, kmod*(D_r * st * sp + D_theta * ct * sp + D_phi * cp));
|
|
|
|
}
|
|
|
|
return ccomplex_t(0.0, kmod*(D_r * ct - D_theta * st));
|
2019-11-06 14:06:19 +01:00
|
|
|
}
|
|
|
|
|
2020-02-28 16:15:37 +01:00
|
|
|
inline real_t vfac_corr( std::array<size_t,3> ijk ) const
|
2019-11-06 14:06:19 +01:00
|
|
|
{
|
|
|
|
real_t ix = ijk[0]*mapratio_, iy = ijk[1]*mapratio_, iz = ijk[2]*mapratio_;
|
2020-02-28 16:15:37 +01:00
|
|
|
const real_t alpha = 1.0/std::real(D_xy_.get_cic_kspace({ix,iy,iz}));
|
|
|
|
return 1.0/alpha;
|
2020-03-01 17:21:17 +01:00
|
|
|
// // below is for LCDM, but it is a tiny correction for typical starting redshifts:
|
2020-02-28 16:15:37 +01:00
|
|
|
//! X = \Omega_\Lambda / \Omega_m
|
|
|
|
// return 1.0 / (alpha - (2*std::pow(aini_,3)*alpha*(2 + alpha)*XmL_*Hypergeometric2F1((3 + alpha)/3.,(5 + alpha)/3.,
|
|
|
|
// (13 + 4*alpha)/6.,-(std::pow(aini_,3)*XmL_)))/
|
|
|
|
// ((7 + 4*alpha)*Hypergeometric2F1(alpha/3.,(2 + alpha)/3.,(7 + 4*alpha)/6.,-(std::pow(aini_,3)*XmL_))));
|
2019-11-05 19:14:14 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
};
|
|
|
|
|
2019-11-03 15:54:17 +01:00
|
|
|
}
|