Changes to integrate radius

AMReX-Astro · Jan 30, 2024 · 094bea9 · 094bea9
1 parent 2f382fc
commit 094bea9
Show file tree

Hide file tree

Showing 7 changed files with 995 additions and 11 deletions.
diff --git a/Exec/ParticleFilterTest/DarkMatterParticleContainer.H b/Exec/ParticleFilterTest/DarkMatterParticleContainer.H
@@ -31,7 +31,7 @@ public:
 
     virtual void moveKickDrift (amrex::MultiFab& acceleration, int level, amrex::Real time,
                                 amrex::Real timestep,
-                                amrex::Real a_old = 1.0, amrex::Real a_half = 1.0, int where_width = 0);
+                                amrex::Real a_old = 1.0, amrex::Real a_half = 1.0, int where_width = 0, amrex::Real radius_inner = -1.e34, amrex::Real radius_outer = -1.e34);
     virtual void moveKick      (amrex::MultiFab& acceleration, int level, amrex::Real time,
                                 amrex::Real timestep,
                                 amrex::Real a_new = 1.0, amrex::Real a_half = 1.0);

diff --git a/Exec/ParticleFilterTest/DarkMatterParticleContainer.cpp b/Exec/ParticleFilterTest/DarkMatterParticleContainer.cpp
@@ -145,7 +145,9 @@ DarkMatterParticleContainer::moveKickDrift (amrex::MultiFab&       acceleration,
                                             amrex::Real            dt,
                                             amrex::Real            a_old,
                                             amrex::Real            a_half,
-                                            int                    where_width)
+                                            int                    where_width,
+                                            amrex::Real            radius_inner,
+                                            amrex::Real            radius_outer)
 {
     BL_PROFILE("DarkMatterParticleContainer::moveKickDrift()");
 
@@ -191,10 +193,8 @@ DarkMatterParticleContainer::moveKickDrift (amrex::MultiFab&       acceleration,
     auto geom_test=pc->Geom(lev);
     const GpuArray<Real,AMREX_SPACEDIM> phi=geom_test.ProbHiArray();
     const GpuArray<Real,AMREX_SPACEDIM> center({AMREX_D_DECL((phi[0]-plo[0])*0.5,(phi[1]-plo[1])*0.5,(phi[2]-plo[2])*0.5)});
-    Real time_unit = 3.0856776e19 / 31557600.0; // conversion to Julian years
-    time_unit=10.0;
-    Real radius_inner=(t)*time_unit*c_light;
-    Real radius_outer=(t+dt)*time_unit*c_light;
+    //const Real a_half = 0.5 * (a_old + a_new);
+
     auto domain=geom_test.Domain();
     auto z=1/a_old-1;
     //From write_info.cpp
@@ -307,6 +307,7 @@ DarkMatterParticleContainer::moveKickDrift (amrex::MultiFab&       acceleration,
     real_comp_names_shell.push_back("yposold");
     real_comp_names_shell.push_back("zposold");
     ShellPC->WritePlotFile(dir, name, real_comp_names_shell);
+    Print()<<"After write\t"<<ShellPC->TotalNumberOfParticles()<<std::endl;
     //    ShellPC->amrex::ParticleContainer<7,0>::WritePlotFile(dir, name, real_comp_names_shell);
     if (ac_ptr != &acceleration) delete ac_ptr;
 

diff --git a/Exec/ParticleFilterTest/Nyx.H b/Exec/ParticleFilterTest/Nyx.H
@@ -251,6 +251,7 @@ public:
     //
     amrex::Real get_comoving_a(amrex::Real time);
     static void integrate_time_given_a(const amrex::Real a0, const amrex::Real a1, amrex::Real& dt);
+    static void integrate_distance_given_a(const amrex::Real a0, const amrex::Real a1, amrex::Real& dt);
     void integrate_comoving_a(const amrex::Real old_a, amrex::Real& new_a, const amrex::Real dt);
     void integrate_comoving_a(amrex::Real time, amrex::Real dt);
 

diff --git a/Exec/ParticleFilterTest/NyxParticleContainer.H b/Exec/ParticleFilterTest/NyxParticleContainer.H
@@ -14,7 +14,8 @@ public:
 
     virtual void moveKickDrift (amrex::MultiFab& acceleration, int level, amrex::Real time,
                                 amrex::Real timestep,
-                        amrex::Real a_old = 1.0, amrex::Real a_half = 1.0, int where_width = 0) = 0;
+                                amrex::Real a_old = 1.0, amrex::Real a_half = 1.0, int where_width = 0,
+                                amrex::Real radius_inner = -1.e34, amrex::Real radius_outer = -1.e34) = 0;
     virtual void moveKick      (amrex::MultiFab& acceleration, int level, amrex::Real time,
                                 amrex::Real timestep,
                                 amrex::Real a_new = 1.0, amrex::Real a_half = 1.0) = 0;

diff --git a/Exec/ParticleFilterTest/Nyx_advance.cpp b/Exec/ParticleFilterTest/Nyx_advance.cpp
@@ -248,8 +248,15 @@ Nyx::advance_hydro_plus_particles (Real time,
             MultiFab grav_vec_old(ba, dm, AMREX_SPACEDIM, grav_n_grow);
             get_level(lev).gravity->get_old_grav_vector(lev, grav_vec_old, time);
 
-            for (int i = 0; i < Nyx::theActiveParticles().size(); i++)
-                Nyx::theActiveParticles()[i]->moveKickDrift(grav_vec_old, lev, time, dt, a_old, a_half, where_width);
+            for (int i = 0; i < Nyx::theActiveParticles().size(); i++) {
+                Real radius_inner = -1.e34;
+                Real radius_outer = -1.e34;
+                integrate_distance_given_a(a_old, 1.0, radius_outer);
+                integrate_distance_given_a(a_new, 1.0, radius_inner);
+                Print()<<radius_inner<<std::endl;
+                Print()<<radius_outer<<std::endl;
+                Nyx::theActiveParticles()[i]->moveKickDrift(grav_vec_old, lev, time, dt, a_old, a_half, where_width, radius_inner, radius_outer);
+            }
 
             // Only need the coarsest virtual particles here.
                 if (lev == level && level < finest_level)

diff --git a/Exec/ParticleFilterTest/Nyx_advance_particles.cpp b/Exec/ParticleFilterTest/Nyx_advance_particles.cpp
@@ -190,8 +190,13 @@ Nyx::advance_particles_only (Real time,
             MultiFab grav_vec_old(ba, dm, AMREX_SPACEDIM, grav_n_grow);
             get_level(lev).gravity->get_old_grav_vector(lev, grav_vec_old, time);
 
-            for (int i = 0; i < Nyx::theActiveParticles().size(); i++)
-                Nyx::theActiveParticles()[i]->moveKickDrift(grav_vec_old, lev, time, dt, a_old, a_half, where_width);
+            for (int i = 0; i < Nyx::theActiveParticles().size(); i++) {
+                Real radius_inner = -1.e34;
+                Real radius_outer = -1.e34;
+                integrate_distance_given_a(1.0, a_old, radius_inner);
+                integrate_distance_given_a(1.0, a_new, radius_outer);
+                Nyx::theActiveParticles()[i]->moveKickDrift(grav_vec_old, lev, time, dt, a_old, a_half, where_width, radius_inner, radius_outer);
+            }
 
             // Only need the coarsest virtual particles here.
             if (lev == level && level < finest_level)