From 7a4b65c66a0ee391603ca7aa747748ad0b85a57f Mon Sep 17 00:00:00 2001
From: Oliver Hahn <oliverjhahn@gmail.com>
Date: Thu, 4 Jul 2019 10:17:46 +0200
Subject: [PATCH] added test problem (hard coded)

---
 src/ic_generator.cc | 36 +++++++++++++++++++++++++++++-------
 1 file changed, 29 insertions(+), 7 deletions(-)

diff --git a/src/ic_generator.cc b/src/ic_generator.cc
index 7530efa..de6c02a 100644
--- a/src/ic_generator.cc
+++ b/src/ic_generator.cc
@@ -53,11 +53,11 @@ int Run( ConfigFile& the_config )
     const real_t Dplus0 = the_cosmo_calc->CalcGrowthFactor(astart) / the_cosmo_calc->CalcGrowthFactor(1.0);
     const real_t vfac   = the_cosmo_calc->CalcVFact(astart);
 
-    const double g1  = Dplus0;
+    const double g1  = -Dplus0;
     const double g2  = (LPTorder>1)? -3.0/7.0*Dplus0*Dplus0 : 0.0;
     const double g3a = (LPTorder>2)? -1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0;
     const double g3b = (LPTorder>2)? 10.0/21.*Dplus0*Dplus0*Dplus0 : 0.0;
-    const double g3c = (LPTorder>2)? -1.0/7.0*Dplus0*Dplus0*Dplus0 : (bSymplecticPT)? g1*g2 : 0.0;
+    const double g3c = (LPTorder>2)? -1.0/7.0*Dplus0*Dplus0*Dplus0 : 0.0;
 
     const double vfac1 =  vfac;
     const double vfac2 =  2*vfac1;
@@ -92,17 +92,20 @@ int Run( ConfigFile& the_config )
 
     //--------------------------------------------------------------------
     
-    //--------------------------------------------------------------------
-    // Fill the grid with a Gaussian white noise field
-    //--------------------------------------------------------------------
-    the_random_number_generator->Fill_Grid( phi );
-
+    
     //======================================================================
     //... compute 1LPT displacement potential ....
     //======================================================================
     // phi = - delta / k^2
     double wtime = get_wtime();    
     csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(1) term" << std::flush;
+
+    #if 0 //  random ICs
+    //--------------------------------------------------------------------
+    // Fill the grid with a Gaussian white noise field
+    //--------------------------------------------------------------------
+    the_random_number_generator->Fill_Grid( phi );
+
     phi.FourierTransformForward();
     
     phi.apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {
@@ -113,6 +116,25 @@ int Run( ConfigFile& the_config )
     });
 
     phi.zero_DC_mode();
+    #else // ICs with a given phi(1) potential function
+    constexpr real_t twopi{2.0*M_PI};
+    constexpr real_t epsilon_q1d{0.25};
+    
+    phi.FourierTransformBackward(false);
+
+    phi.apply_function_r_dep([&](auto v, auto r) -> real_t {
+        real_t q2 = r[0]/boxlen * twopi;
+        real_t q1 = r[1]/boxlen * twopi;
+
+        // std::cerr << q1  << " " << q2 << std::endl;
+        
+        return (std::cos(q1) + epsilon_q1d * std::cos(q2))/twopi/twopi; ////( std::cos(q1) + epsilon_q1d * std::cos(q2) );// / volfac;
+        //return -delta / (kmod * kmod) / volfac;
+    });
+    phi.FourierTransformForward();
+
+
+    #endif
     csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
     
     //======================================================================