Tagging based on vorticity made simulation slower (time taken per timestep increases) - AMReX incflo

[brahmandas]

Hi,

I modified the incflo tagging code for cell tagging based on the vorticity magnitude (code attached below with this message). However, every time I simulate a case with this function enabled, the time taken per timestep increases exponentially (as shown in the attached figure). This happens for every simulation and is not unique to a specific timestep or time reached. Do you know what could be the cause in this case? Since the number of cells due to the adaptive mesh refinement (based on vorticity) did not change significantly (the timestep that took 50s has similar number of cells for the timestep which took only 5s).

Any help/suggestion would be greatly appreciated, thank you!

Here is the code I made:

#include <incflo.H>
#include <AMReX_buildInfo.H>
#include <AMReX_MultiFab.H>
#ifdef AMREX_USE_EB
#include <AMReX_EBAmrUtil.H>
#endif

using namespace amrex;

// tag cells for refinement
// overrides the pure virtual function in AmrCore
void incflo::ErrorEst (int lev, TagBoxArray& tags, Real time, int /*ngrow*/)
{
    BL_PROFILE("incflo::ErrorEst()");

    static bool first = true;
    static Vector<Real> rhoerr_v, gradrhoerr_v, velerr_v, gradvelerr_v, vorterr_v;

    static bool tag_region;

    if (first) {
        first = false;
        ParmParse pp("incflo");

        pp.queryarr("rhoerr", rhoerr_v);
        if (!rhoerr_v.empty()) {
            Real last = rhoerr_v.back();
            rhoerr_v.resize(max_level+1, last);
        }
 
 	pp.queryarr("velerr", velerr_v);
        if (!velerr_v.empty()) {
            Real last = velerr_v.back();
            velerr_v.resize(max_level+1, last);
        }

        pp.queryarr("gradrhoerr", gradrhoerr_v);
        if (!gradrhoerr_v.empty()) {
            Real last = gradrhoerr_v.back();
            gradrhoerr_v.resize(max_level+1, last);
        }
	
	pp.queryarr("gradvelerr", gradvelerr_v);
        if (!gradvelerr_v.empty()) {
            Real last = gradvelerr_v.back();
            gradvelerr_v.resize(max_level+1, last);
        }
	
	pp.queryarr("vorterr", vorterr_v);
        if (!vorterr_v.empty()) {
            Real last = vorterr_v.back();
            vorterr_v.resize(max_level+1, last);
        }

        tag_region_lo.resize(3);
        tag_region_hi.resize(3);

        tag_region = false;
        pp.query("tag_region", tag_region);

        pp.queryarr("tag_region_lo", tag_region_lo);
        pp.queryarr("tag_region_hi", tag_region_hi);
    }

    const auto   tagval = TagBox::SET;

    bool tag_rho = lev < rhoerr_v.size();
    bool tag_gradrho = lev < gradrhoerr_v.size();
    bool tag_gradvel = lev < gradvelerr_v.size();
    bool tag_vel     = lev < velerr_v.size();
    bool tag_vort    = lev < vorterr_v.size();

    if (tag_gradrho) {
        fillpatch_density(lev, time, m_leveldata[lev]->density, 1);
    }
    if (tag_gradvel || tag_vort) {
	fillpatch_velocity(lev, time, m_leveldata[lev]->velocity, 1);
    }

    AMREX_D_TERM(const Real l_dx = geom[lev].CellSize(0);,
                 const Real l_dy = geom[lev].CellSize(1);,
                 const Real l_dz = geom[lev].CellSize(2););

#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for (MFIter mfi(m_leveldata[lev]->velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
    	Box const& bx = mfi.tilebox();
	auto const& tag = tags.array(mfi);
	Vector<MultiFab> mf(finest_level + 1);

	if (tag_gradvel || tag_vel || tag_vort)
	{
	   Array4<Real const> const& vel = m_leveldata[lev]->velocity.const_array(mfi);
	   Real gradvelerr = tag_gradvel ? gradvelerr_v[lev]: std::numeric_limits<Real>::max();
	   Real velerr = tag_vel ? velerr_v[lev]: std::numeric_limits<Real>::max();
	   Real vorterr = tag_vort ? vorterr_v[lev]: std::numeric_limits<Real>::max();

	   amrex::ParallelFor (bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
	   {
	   	if (tag_vort) {
		Real vx = (Real(0.5) * (vel(i+1,j,k,1) - vel(i-1,j,k,1))) / l_dx;
		Real wx = (Real(0.5) * (vel(i+1,j,k,2) - vel(i-1,j,k,2))) / l_dx;

		Real uy = (Real(0.5) * (vel(i,j+1,k,0) - vel(i,j-1,k,0))) / l_dy;
		Real wy = (Real(0.5) * (vel(i,j+1,k,2) - vel(i,j-1,k,2))) / l_dy;

#if (AMREX_SPACEDIM == 2)
		if (tag_vort && (vx-uy)  > vorterr) {
		   tag(i,j,k) = tagval;
		   }
#elif (AMREX_SPACEDIM == 3)
		Real uz = (Real(0.5) * (vel(i,j,k+1,0) - vel(i,j,k-1,0))) / l_dz;
		Real vz = (Real(0.5) * (vel(i,j,k+1,1) - vel(i,j,k-1,1))) / l_dz;

		if (tag_vort && std::sqrt((wy-vz)*(wy-vz) + (uz-wx)*(uz-wx) + (vx-uy)*(vx-uy)) > vorterr) {
		   tag(i,j,k) = tagval;
		}
#endif
		}

	   	if (tag_vel && vel(i,j,k) > velerr) {
		   tag(i,j,k) = tagval;
		}

	   	if (tag_gradvel) {
	   	// gradvel in x-dir
	 	Real ux = amrex::Math::abs(vel(i+1,j,k,0) - vel(i,j,k,0));
	   	Real uy = amrex::Math::abs(vel(i,j+1,k,0) - vel(i,j,k,0));
	   	ux = amrex::max(ux, amrex::Math::abs(vel(i,j,k,0) - vel(i-1,j,k,0)));
	   	uy = amrex::max(uy, amrex::Math::abs(vel(i,j,k,0) - vel(i,j-1,k,0)));
	   	// gradvel in y-dir
	   	Real vx = amrex::Math::abs(vel(i+1,j,k,1) - vel(i,j,k,1));
		Real vy = amrex::Math::abs(vel(i,j+1,k,1) - vel(i,j,k,1));
		vx = amrex::max(vx, amrex::Math::abs(vel(i,j,k,1) - vel(i-1,j,k,1)));
		vy = amrex::max(vy, amrex::Math::abs(vel(i,j,k,1) - vel(i,j-1,k,1)));
		// gradvel in z-dir
		Real wx = amrex::Math::abs(vel(i+1,j,k,2) - vel(i,j,k,2));
		Real wy = amrex::Math::abs(vel(i,j+1,k,2) - vel(i,j,k,2));
		wx = amrex::max(wx, amrex::Math::abs(vel(i,j,k,2) - vel(i-1,j,k,2)));
		wy = amrex::max(wy, amrex::Math::abs(vel(i,j,k,2) - vel(i,j-1,k,2)));
#if (AMREX_SPACEDIM == 2)
	   	if (amrex::max(ux, uy, vx, vy, wx, wy) >= gradvelerr) {
	   	   tag(i,j,k) = tagval;
	   	}
#elif (AMREX_SPACEDIM == 3)
	   	//gradvel in x-dir
	   	Real uz = amrex::Math::abs(vel(i,j,k+1,0) - vel(i,j,k,0));
	   	uz = amrex::max(uz, amrex::Math::abs(vel(i,j,k,0) - vel(i,j,k-1,0)));
	   	//gradvel in y-dir	   
	   	Real vz = amrex::Math::abs(vel(i,j,k+1,1) - vel(i,j,k,1));
	   	vz = amrex::max(vz, amrex::Math::abs(vel(i,j,k,1) - vel(i,j,k-1,1)));
	   	//gradvel in z-dir
	   	Real wz = amrex::Math::abs(vel(i,j,k+1,2) - vel(i,j,k,2));
	   	wz = amrex::max(wz, amrex::Math::abs(vel(i,j,k,2) - vel(i,j,k-1,2)));
	   
	   	if (amrex::max(ux, uy, uz, vx, vy, vz, wx, wy, wz) >= gradvelerr) {
	   	   tag(i,j,k) = tagval;
	   	}
#endif
		}
	   });
	}
    }

    for (MFIter mfi(m_leveldata[lev]->density,TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        Box const& bx = mfi.tilebox();
        auto const& tag = tags.array(mfi);
	Vector<MultiFab> mf(finest_level + 1);

        if (tag_rho || tag_gradrho)
        {
            Array4<Real const> const& rho = m_leveldata[lev]->density.const_array(mfi);
            Real rhoerr = tag_rho ? rhoerr_v[lev]: std::numeric_limits<Real>::max();
            Real gradrhoerr = tag_gradrho ? gradrhoerr_v[lev] : std::numeric_limits<Real>::max();
            amrex::ParallelFor(bx,
            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                if (tag_rho && rho(i,j,k) > rhoerr) {
                    tag(i,j,k) = tagval;
                }
                if (tag_gradrho) {
                    Real ax = amrex::Math::abs(rho(i+1,j,k) - rho(i,j,k));
                    Real ay = amrex::Math::abs(rho(i,j+1,k) - rho(i,j,k));
                    ax = amrex::max(ax,amrex::Math::abs(rho(i,j,k) - rho(i-1,j,k)));
                    ay = amrex::max(ay,amrex::Math::abs(rho(i,j,k) - rho(i,j-1,k)));
#if (AMREX_SPACEDIM == 2)
                    if (amrex::max(ax,ay) >= gradrhoerr) {
                        tag(i,j,k) = tagval;
                    }
#elif (AMREX_SPACEDIM == 3)
                    Real az = amrex::Math::abs(rho(i,j,k+1) - rho(i,j,k));
                    az = amrex::max(az,amrex::Math::abs(rho(i,j,k) - rho(i,j,k-1)));
                    if (amrex::max(ax,ay,az) >= gradrhoerr) {
                        tag(i,j,k) = tagval;
                    }
#endif
                }
            });
        }

        if (tag_region) {

            Real xlo = tag_region_lo[0];
            Real ylo = tag_region_lo[1];
            Real xhi = tag_region_hi[0];
            Real yhi = tag_region_hi[1];
            auto const& problo = geom[lev].ProbLoArray();

#if (AMREX_SPACEDIM == 2)

            amrex::ParallelFor(bx,
            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                 Real x = problo[0] + (i+0.5)*l_dx;
                 Real y = problo[1] + (j+0.5)*l_dy;

                 // Tag if we are inside the specified box
                 if (x >= xlo && x <= xhi && y >= ylo && y <= yhi)
                 {
                    tag(i,j,k) = tagval;
                 }
            });

#else
            Real zlo = tag_region_lo[2];
            Real zhi = tag_region_hi[2];

            amrex::ParallelFor(bx,
            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                 Real x = problo[0] + Real(i+0.5)*l_dx;
                 Real y = problo[1] + Real(j+0.5)*l_dy;
                 Real z = problo[2] + Real(k+0.5)*l_dz;

                 // Tag if we are inside the specified box
                 if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi)
                 {
                    tag(i,j,k) = tagval;
                 }
            });
#endif
        }
    }

#ifdef AMREX_USE_EB
    m_refine_cutcells = true;
    // Refine on cut cells
    if (m_refine_cutcells)
    {
        amrex::TagCutCells(tags, m_leveldata[lev]->velocity);
    }
#endif
}

[MaybeBin]

I encounter a similar issue! I endeavor to utilize my program, which is based on AMReX, to simulate the evolution of a flame using over 1000 CPU cores in parallel computation. I have observed that the duration of each time step progressively increases as the calculation progresses. However, when I halt the calculation and resume from the checkpoint file, the time required for each time step initially decreases and then gradually increases again.
These images illustrate the phenomenon: it takes approximately 10 hours to complete 400 steps, but after the restart from chk12100, the calculation time reduces to only 1 hour. This perplexes me greatly, as I am uncertain whether I am encountering the same issues as you or if it pertains to the efficiency of parallelization. I am so confused in which direction should I strive to resolve this issue?

before restart from the checkpoint

after restart from the checkpoint

i compiled my code with:

openmpi/4.1.2
gcc/11.2.0
AMReX/23.11

I genuinely appreciate for any response and assistance. Thank you for your prompt reply.

[WeiqunZhang]

Could you do some profiling? Is it slow everywhere or just certain parts of the code? What does your code use? Linear solves? Particles? EB? ...

[MaybeBin]

Thank you for your response : ). I have utilized the CNS example code as the foundation for my own implementation, augmenting it with an additional source term. After conducting code profiling, I have identified the prominent areas of concern, which predominantly involve several MPI communication functions within the Runge-Kutta process, as illustrated below:

TinyProfiler total time across processes [min...avg...max]: 1.324e+04 ... 1.324e+04 ... 1.324e+04

----------------------------------------------------------------------------------------------
Name                                           NCalls  Excl. Min  Excl. Avg  Excl. Max   Max %
----------------------------------------------------------------------------------------------
FabArray::ParallelCopy_finish()                 34785       4181       4224       4269  32.26%
FAST::computeRHS_box                         14937194       3054       3088       3114  23.53%
FillBoundary_finish()                            9000       1974       1994       2042  15.43%
FAST::computeRHS()                              14000       1028       1042       1058   7.99%
FabArray::ParallelCopy_nowait()                 34785      888.6      900.7      926.9   7.00%
FillBoundary_nowait()                            9000      349.1      361.2      375.1   2.83%
FillPatchInterp(MF)                             13760      315.9      349.8        374   2.83%
MultiFab::max()                                 18001      187.5      236.3      282.3   2.13%
FabArrayBase::CPC::define()                     16282      266.7      270.8        274   2.07%
VisMF::Read()                                      12      158.4      159.6      160.2   1.21%
FillPatcher::fillCFB                            12000      128.6      142.9      153.6   1.16%
Amr::regrid()                                     500      65.52      67.28      67.98   0.51%
FAST::reactionRHS()                              7000      51.44       53.2      55.13   0.42%
.......
.......
----------------------------------------------------------------------------------------------
Name                                           NCalls  Incl. Min  Incl. Avg  Incl. Max   Max %
----------------------------------------------------------------------------------------------
main()                                              1  1.324e+04  1.324e+04  1.324e+04 100.00%
Amr::coarseTimeStep()                             500  1.306e+04  1.306e+04  1.306e+04  98.67%
Amr::timeStep()                                  3500  1.302e+04  1.302e+04  1.302e+04  98.41%
FAST::advance()                                  3500  1.213e+04  1.213e+04  1.213e+04  91.65%
AmrLevel::RK()                                   3500  1.178e+04  1.182e+04  1.187e+04  89.67%
RungeKutta3                                      3500  1.178e+04  1.182e+04  1.187e+04  89.67%
AmrLevel::FillPatcherFill()                     14000       7596       7642       7669  57.94%
FillPatcher::fill()                             12000       7444       7490       7517  56.79%
FabArray::ParallelCopy()                        34785       5333       5385       5426  41.00%
FillPatcher::fillCFB                            12000       4214       4261       4315  32.60%
FabArray::ParallelCopy_finish()                 34785       4181       4224       4269  32.26%
FAST::computeRHS()                              14000       4091       4137       4166  31.48%
FillPatchSingleLevel                            21775       3860       3891       3954  29.88%
FAST::computeRHS_box                         14937194       3054       3088       3114  23.53%
FabArray::FillBoundary()                         9000       2349       2372       2421  18.29%
FillBoundary_finish()                            9000       1974       1994       2042  15.43%
FabArray::ParallelCopy_nowait()                 34785       1145       1161       1190   8.99%
Amr::regrid()                                     500      788.7      788.9      789.1   5.96%
AmrLevel::FillPatch()                            6515      758.2      764.8      770.7   5.82%
FillPatchIterator::Initialize                    6515      748.9      755.5      761.6   5.75%
FillPatchIterator::FillFromLevel0()              3255      443.4      447.8      450.3   3.40%
FillBoundary_nowait()                            9000      366.5      378.6      392.5   2.97%
FillPatchInterp(MF)                             13760      315.9      349.8        374   2.83%
FAST::reactionRHS()                              7000      237.9      285.5      330.1   2.49%
FillPatchIterator::FillFromTwoLevels()           3260      299.9      302.1      306.6   2.32%
FillPatchTwoLevels                               3260      299.8      302.1      306.6   2.32%
MultiFab::max()                                 18001      187.5      236.3      282.3   2.13%
FabArrayBase::getCPC()                          36039      266.8      270.9      274.2   2.07%
FabArrayBase::CPC::define()                     16282      266.7      270.8        274   2.07%
Amr::init()                                         1      160.5      160.8      161.1   1.22%
Amr::restart()                                      1      160.5      160.8      161.1   1.22%
AmrLevel::restart()                                 3      158.7      159.8      160.3   1.21%
StateData::restartDoit()                            9      158.4      159.6      160.2   1.21%
VisMF::Read()                                      12      158.4      159.6      160.2   1.21%
.......
.......

I realized that the issues I encountered may have resulted from my improper input. The contents of my input file are displayed below. The total number of grids reaches several hundred million across both levels, and I employ 1024 CPU cores to perform parallel calculations. I am uncertain whether the small value of max_grid_size=32 contributes to the excessive subdivision of the calculation domain into numerous smaller boxes relative to the number of processes. This, in turn, leads to MPI communication between the boxes dominating the calculation time. However, I attempted increasing the value of max_grid_size, and it did result in faster computations.
I am interested in identifying the parameters that can enhance parallel efficiency on CPU or GPU, as well as any considerations I should bear in mind when configuring parallel computing settings. Unfortunately, I did not come across some relevant information in the source documentation or related issues. Thank you for your assistance!

geometry.is_periodic = 0 0 0
geometry.coord_sys   = 0 
geometry.prob_lo     =   0.0     0.0      0.0
geometry.prob_hi     =   0.53    0.053    0.053
amr.n_cell           =   1680      168    168
.......
.......
amr.max_level       = 2   
amr.ref_ratio       = 2 2 2 2 2
amr.regrid_int      = 2 2 2 2 2 
amr.blocking_factor = 8
amr.max_grid_size   = 32    
amr.n_error_buf     = 2 2 2 2 2
amr.grid_eff        = 0.8    

[WeiqunZhang]

Could you print out the boxarrays on each level? There might be a lot of small boxes because the number of cells amr.n_cell = 1680 168 168 is not a multiple of max_grid_size=32 (not even 16).

[MaybeBin]

Thank you for your valuable advice. I will make an effort to investigate the interrelationship between the number of boxes and the parallel process.

[WeiqunZhang]

@brahmandas Could you tell how to reproduce this with incflo? The code you posted has some dead code like Vector<MultiFab> mf(finest_level + 1); inside MFIter loops. Not sure why it was there in the first place. So if possible, please provide the git fork of inclfo.

[WeiqunZhang]

I cannot seem to reproduce what you saw.

[brahmandas]

Hi @WeiqunZhang
I looked into this more and found that the issue is system dependent. Out of four systems I tested this, I found the issue on two systems. May I ask if there is an amrex flag or compiler setting that can lead to such performance issues?
My logs from two machines are attached. The system where I had no issues and the simulation run faster has AMD 7713 Milan 2GHz processors whereas the system where I encountered issue has 3.0 GHz Intel Xeon Gold 6154 processors.

log_INTEL.txt
log_AMD.txt

[WeiqunZhang]

I cannot think of anything that could explain this. Could you compile with TINY_PROFILE=TRUE and run for a few hundred steps (by changing either stop_time or max_step)? Then we can look at the profiling result.

[brahmandas]

@WeiqunZhang This is the case with TINY_PROFILE=TRUE for the same case as before (same restarting point) with 100 timesteps of run. The 'log_good' shows the simulation with no increasing time per timestep and 'log_bad' shows the increasing time. The profiling result is at the end of the log file. I also could not find specifically how the different systems could affect this. Any insight/help would be appreciated.

log_bad.txt
log_good.txt

[WeiqunZhang]

I think too many MPI processes were used for this small problem. When I was testing it, I was able to run it on my desktop computer using only 4 MPI processes.