Merge branch 'development'

Exec/LyA/GNUmakefile

@@ -31,6 +31,7 @@ USE_GRAV = TRUE
 USE_HEATCOOL = TRUE
 USE_AGN = FALSE
 USE_CVODE = FALSE
+USE_SUNDIALS3 = FALSE
 
 Bpack := ./Make.package
 Blocs := .

Exec/Scaling/inputs

@@ -35,8 +35,8 @@ amr.data_log = runlog
 gravity.gravity_type = PoissonGrav
 gravity.no_sync      = 1
 gravity.no_composite = 1
-gravity.solve_with_cpp = 0
-gravity.solve_with_hpgmg = 1
+gravity.solve_with_mlmg = 1
+gravity.solve_with_hpgmg = 0
 
 mg.bottom_solver = 4
 

Source/AGNParticleContainer.cpp

@@ -279,8 +279,8 @@ void AGNParticleContainer::writeAllAtLevel(int lev)
     {
       auto& particles = pti.GetArrayOfStructs();
       size_t Np = pti.numParticles();
-      pout() << "There are " << Np  << " AGN particles at level " << lev
-             << " in grid " << pti.index() << std::endl;
+      Print() << "There are " << Np  << " AGN particles at level " << lev
+              << " in grid " << pti.index() << std::endl;
       for (unsigned i = 0; i < Np; ++i)
         {
           const ParticleType& p = particles[i];
@@ -294,14 +294,14 @@ void AGNParticleContainer::writeAllAtLevel(int lev)
           Real energy = p.rdata(4);
           Real mdot = p.rdata(5);
 
-          pout() << "[" << i << "]: id " << id << " cpu " << cpu
-                 << " mass " << mass
-                 << " index " << iv
-                 << " position " << xyz
-                 << " velocity " << uvw
-                 << " energy " << energy
-                 << " mdot " << mdot
-                 << endl;
+          Print() << "[" << i << "]: id " << id << " cpu " << cpu
+                  << " mass " << mass
+                  << " index " << iv
+                  << " position " << xyz
+                  << " velocity " << uvw
+                  << " energy " << energy
+                  << " mdot " << mdot
+                  << endl;
         }
     }
 }

Source/HeatCool/Make.package

@@ -5,21 +5,32 @@ f90EXE_sources += integrate_state_vode_3d.f90
 ifeq ($(USE_SDC), TRUE)
 f90EXE_sources += integrate_state_with_source_3d.f90
   ifeq ($(USE_CVODE), TRUE)
-    f90EXE_sources += fcvode_extras_src.f90
-    f90EXE_sources += integrate_state_fcvode_with_source_3d.f90
+      f90EXE_sources += fcvode_extras_src.f90
+      f90EXE_sources += integrate_state_fcvode_with_source_3d.f90
   else
-    f90EXE_sources += integrate_state_fcvode_with_source_3d_stubs.f90
+	ifeq ($(USE_SUNDIALS3), TRUE)
+	     f90EXE_sources += fcvode_extras_src.f90
+	     f90EXE_sources += integrate_state_fcvode_with_source_3d.f90
+	 else
+	     f90EXE_sources += integrate_state_fcvode_with_source_3d_stubs.f90
+	 endif
   endif
 endif
 
 ifeq ($(USE_CVODE), TRUE)
   f90EXE_sources += fcvode_extras.f90
   f90EXE_sources += integrate_state_fcvode_3d.f90
   f90EXE_sources += integrate_state_fcvode_vec_3d.f90
+else
+ifeq ($(USE_SUNDIALS3), TRUE)
+  f90EXE_sources += fcvode_extras.f90
+  f90EXE_sources += integrate_state_fcvode_3d.f90
+  f90EXE_sources += integrate_state_fcvode_vec_3d.f90
 else
   f90EXE_sources += integrate_state_fcvode_3d_stubs.f90
   f90EXE_sources += integrate_state_fcvode_vec_3d_stubs.f90
 endif
+endif
 f90EXE_sources += vode_aux.f90
 f90EXE_sources += f_rhs.f90
 else

Source/HydroFortran/Nyx_advection_3d.f90

@@ -694,12 +694,12 @@ subroutine ctoprim(lo,hi,          uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3
                          src,  src_l1, src_l2, src_l3, src_h1, src_h2, src_h3, &
                          srcQ,srcq_l1,srcq_l2,srcq_l3,srcq_h1,srcq_h2,srcq_h3, &
                          grav,gv_l1, gv_l2, gv_l3, gv_h1, gv_h2, gv_h3, &
-                         dx,dy,dz,dt,ngp,ngf,a_old,a_new)
+                         dx,dy,dz,dt,ngq,ngf,a_old,a_new)
       !
-      !     Will give primitive variables on lo-ngp:hi+ngp, and flatn on lo-ngf:hi+ngf
+      !     Will give primitive variables on lo-ngq:hi+ngq, and flatn on lo-ngf:hi+ngf
       !     if use_flattening=1.  Declared dimensions of q,c,csml,flatn are given
-      !     by DIMS(q).  This declared region is assumed to encompass lo-ngp:hi+ngp.
-      !     Also, uflaten call assumes ngp>=ngf+3 (ie, primitve data is used by the
+      !     by DIMS(q).  This declared region is assumed to encompass lo-ngq:hi+ngq.
+      !     Also, uflaten call assumes ngq>=ngf+3 (ie, primitve data is used by the
       !     routine that computes flatn).
       !
       use amrex_error_module
@@ -739,16 +739,16 @@ subroutine ctoprim(lo,hi,          uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3
       real(rt) :: dpdr, dpde
 
       integer          :: i, j, k
-      integer          :: ngp, ngf, loq(3), hiq(3)
+      integer          :: ngq, ngf, loq(3), hiq(3)
       integer          :: n, nq
       integer          :: iadv, ispec
       real(rt) :: courx, coury, courz
       real(rt) :: a_half, a_dot, rhoInv
       real(rt) :: dtdxaold, dtdyaold, dtdzaold, small_pres_over_dens
 
       do i=1,3
-         loq(i) = lo(i)-ngp
-         hiq(i) = hi(i)+ngp
+         loq(i) = lo(i)-ngq
+         hiq(i) = hi(i)+ngq
       enddo
       !
       ! Make q (all but p), except put e in slot for rho.e, fix after eos call.
@@ -862,9 +862,9 @@ subroutine ctoprim(lo,hi,          uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3
       !        IF NOT THEN THE FORMULAE BELOW ARE INCOMPLETE.
 
       ! compute srcQ terms
-      do k = lo(3)-1, hi(3)+1
-         do j = lo(2)-1, hi(2)+1
-            do i = lo(1)-1, hi(1)+1
+      do k = loq(3), hiq(3)
+         do j = loq(2), hiq(2)
+            do i = loq(1), hiq(1)
 
                rhoInv = ONE/q(i,j,k,QRHO)
 

Source/HydroFortran/make_hydro_sources_3d.f90

@@ -75,19 +75,19 @@ subroutine fort_make_hydro_sources(time,lo,hi,&
       q_h2 = hi(2)+NHYP
       q_h3 = hi(3)+NHYP
 
-      srcq_l1 = lo(1)-1
-      srcq_l2 = lo(2)-1
-      srcq_l3 = lo(3)-1
-      srcq_h1 = hi(1)+1
-      srcq_h2 = hi(2)+1
-      srcq_h3 = hi(3)+1
+      srcq_l1 = lo(1)-NHYP
+      srcq_l2 = lo(2)-NHYP
+      srcq_l3 = lo(3)-NHYP
+      srcq_h1 = hi(1)+NHYP
+      srcq_h2 = hi(2)+NHYP
+      srcq_h3 = hi(3)+NHYP
 
       call amrex_allocate(     q, lo-NHYP, hi+NHYP, QVAR)
       call amrex_allocate( flatn, lo-NHYP, hi+NHYP      )
       call amrex_allocate(     c, lo-NHYP, hi+NHYP      )
       call amrex_allocate(  csml, lo-NHYP, hi+NHYP      )
 
-      call amrex_allocate(   srcQ, lo-1, hi+1, QVAR)
+      call amrex_allocate(   srcQ, lo-NHYP, hi+NHYP, QVAR)
       call amrex_allocate(divu_nd, lo  , hi+1)
 
       dx = delta(1)

Source/Nyx.cpp

@@ -438,11 +438,15 @@ Nyx::read_params ()
     }
 
 #ifndef AMREX_USE_CVODE
+    #ifndef AMREX_USE_SUNDIALS3
     if (heat_cool_type == 5 || heat_cool_type == 7)
-        amrex::Error("Nyx:: cannot set heat_cool_type = 5 or 7 unless USE_CVODE=TRUE");
+        amrex::Error("Nyx:: cannot set heat_cool_type = 5 or 7 unless USE_CVODE=TRUE or USE_SUNDIALS3=TRUE");
+    #endif
 #else
+    #ifdef SDC
     if (heat_cool_type == 7 && sdc_split == 1)
         amrex::Error("Nyx:: cannot set heat_cool_type = 7 with sdc_split = 1");
+    #endif
 #endif
 
 #else
@@ -1771,7 +1775,7 @@ Nyx::post_regrid (int lbase,
         if (gravity->get_gravity_type() == "PoissonGrav")
 #endif
         {
-            int ngrow_for_solve = 1;
+            int ngrow_for_solve = parent->levelCount(level) + 1;
             int use_previous_phi_as_guess = 1;
             gravity->multilevel_solve_for_new_phi(level, new_finest, ngrow_for_solve, use_previous_phi_as_guess);
         }

Source/NyxParticleContainer.H

@@ -61,10 +61,9 @@ public:
     { 
         amrex::NeighborParticleContainer<NSR,NSI>::AssignCellDensitySingleLevel(0, mf, level, ncomp, particle_lvl_offset);
     }
-    virtual void AssignDensity (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& mf, int lev_min = 0, int ncomp = 1, int finest_level = -1,
-                                int ngrow = 1) const override
+    virtual void AssignDensity (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& mf, int lev_min = 0, int ncomp = 1, int finest_level = -1, int ngrow = 1) const override
     {
-       amrex::NeighborParticleContainer<NSR,NSI>::AssignDensity(0, mf, lev_min, ncomp, finest_level);
+        amrex::NeighborParticleContainer<NSR,NSI>::AssignDensity(0, mf, lev_min, ncomp, finest_level, ngrow);
     }
 
     void MultiplyParticleMass (int lev, amrex::Real mult);

Source/NyxParticles.cpp

@@ -993,7 +993,7 @@ Nyx::particle_redistribute (int lbase, bool my_init)
                    amrex::Print() << "Calling redistribute because DistMap changed " << '\n';
             }
 
-            DMPC->Redistribute(lbase);
+            DMPC->Redistribute(lbase, flev, parent->levelCount(lbase));
             //
             // Use the new BoxArray and DistMap to define ba and dm for next time.
             //

Source/Nyx_halos.cpp

@@ -290,7 +290,7 @@ Nyx::halo_find (Real dt)
        // Convert new_state to conserved variables: rho, momentum, energy.
        primitive_to_conserved(new_state);
 
-       //pout() << "Going into ComputeParticleVelocity (no energy), number of AGN particles on this proc is "
+       //Print() << "Going into ComputeParticleVelocity (no energy), number of AGN particles on this proc is "
        //       << Nyx::theAPC()->TotalNumberOfParticles(true, true) << endl;
 
        // Re-set the particle velocity (but not energy) after accretion,
@@ -300,7 +300,7 @@ Nyx::halo_find (Real dt)
        int add_energy = 0;
        Nyx::theAPC()->ComputeParticleVelocity(level, orig_state, new_state, add_energy);
 
-       //pout() << "Going into ReleaseEnergy, number of AGN particles on this proc is "
+       //Print() << "Going into ReleaseEnergy, number of AGN particles on this proc is "
        //       << Nyx::theAPC()->TotalNumberOfParticles(true, true) << endl;
 
        // AGN particles: may zero out energy.
@@ -310,7 +310,7 @@ Nyx::halo_find (Real dt)
        Nyx::theAPC()->ReleaseEnergy(level, new_state, D_new, a);
        // Now new_state = get_new_data(State_Type) has been updated.
 
-       //       pout() << "Step " << nStep() << " at end of Nyx_halos:" << endl;
+       //       Print() << "Step " << nStep() << " at end of Nyx_halos:" << endl;
        //       Nyx::theAPC()->writeAllAtLevel(level);
 
 #ifdef REEBER
@@ -332,7 +332,7 @@ Nyx::halo_merge ()
    Nyx::theAPC()->Redistribute(level, level, nghost0);
 
    int nmerge = npart -Nyx::theAPC()->TotalNumberOfParticles(true, true);
-   pout() << nStep() << '\t' << nmerge << endl;
+   Print() << nStep() << '\t' << nmerge << endl;
 
 }
 
@@ -356,7 +356,7 @@ Nyx::halo_accrete (Real dt)
    MultiFab agn_density_lost(simBA, simDM, ncomp1, nghost1);
    agn_density_lost.setVal(0.0);
 
-   //pout() << "Going into AccreteMass, number of AGN particles on this proc is "
+   //Print() << "Going into AccreteMass, number of AGN particles on this proc is "
    //       << Nyx::theAPC()->TotalNumberOfParticles(true, true) << endl;
 
    // AGN particles: increase mass and energy.
@@ -378,7 +378,7 @@ Nyx::halo_accrete (Real dt)
    // using change of momentum density in state.
    // No change to state, other than filling ghost cells.
    int add_energy = 1;
-   //pout() << "Going into ComputeParticleVelocity (and energy), number of AGN particles on this proc is "
+   //Print() << "Going into ComputeParticleVelocity (and energy), number of AGN particles on this proc is "
    //       << Nyx::theAPC()->TotalNumberOfParticles(true, true) << endl;
    Nyx::theAPC()->ComputeParticleVelocity(level, orig_state, new_state, add_energy);
    // Now new_state = get_new_data(State_Type) has been updated.

Source/Nyx_nd.f90

@@ -383,7 +383,6 @@ subroutine fort_set_method_params( &
            enddo
            npassive = npassive + nadv
            if(QFS > -1) then
-              print *,'NSPEC NAUX ',nspec, naux
               do ispec = 1, nspec+naux
                  upass_map(npassive + ispec) = UFS + ispec - 1
                  qpass_map(npassive + ispec) = QFS + ispec - 1

Source/sdc_hydro.cpp

@@ -45,14 +45,14 @@ Nyx::sdc_hydro (Real time,
     // point
     enforce_nonnegative_species(S_old);
 
-    MultiFab ext_src_old(grids, dmap, NUM_STATE, 3);
+    MultiFab ext_src_old(grids, dmap, NUM_STATE, NUM_GROW);
     ext_src_old.setVal(0.);
 
     if (add_ext_src)
        get_old_source(prev_time, dt, ext_src_old);
 
     // Define the gravity vector
-    MultiFab grav_vector(grids, dmap, BL_SPACEDIM, 3);
+    MultiFab grav_vector(grids, dmap, BL_SPACEDIM, NUM_GROW);
     grav_vector.setVal(0.);
 
 #ifdef GRAVITY

Source/strang_hydro.cpp

@@ -60,14 +60,14 @@ Nyx::strang_hydro (Real time,
     // point
     enforce_nonnegative_species(S_old);
 
-    MultiFab ext_src_old(grids, dmap, NUM_STATE, 3);
+    MultiFab ext_src_old(grids, dmap, NUM_STATE, NUM_GROW);
     ext_src_old.setVal(0.);
 
     if (add_ext_src)
        get_old_source(prev_time, dt, ext_src_old);
 
     // Define the gravity vector 
-    MultiFab grav_vector(grids, dmap, BL_SPACEDIM, 3);
+    MultiFab grav_vector(grids, dmap, BL_SPACEDIM, NUM_GROW);
     grav_vector.setVal(0.);
 
 #ifdef GRAVITY