Skip to content

Commit

Permalink
Merge branch 'development' into comments-fab
Browse files Browse the repository at this point in the history
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed Nov 28, 2023
2 parents 0858ecb + 9e35dc1 commit 52a6f70
Show file tree
Hide file tree
Showing 8 changed files with 504 additions and 89 deletions.
3 changes: 2 additions & 1 deletion .github/workflows/ascent.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,8 @@ jobs:
-DCMAKE_BUILD_TYPE=Debug \
-DAMReX_ENABLE_TESTS=ON \
-DAMReX_FORTRAN=OFF \
-DAMReX_ASCENT=ON
-DAMReX_ASCENT=ON \
-DAMReX_CONDUIT=ON
- name: Build
run: |
. /ascent_docker_setup_env.sh
Expand Down
5 changes: 2 additions & 3 deletions Src/Extern/Conduit/AMReX_Conduit_Blueprint.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,9 +96,8 @@ namespace amrex
// coordset and fields used to represent the passed particle container.
// This allows you to use unique names to wrap multiple particle containers
// into a single blueprint tree.
template <int NStructReal, int NStructInt, int NArrayReal, int NArrayInt>
void ParticleContainerToBlueprint (const ParticleContainer<NStructReal,
NStructInt,
template <typename ParticleType, int NArrayReal, int NArrayInt>
void ParticleContainerToBlueprint (const ParticleContainer_impl<ParticleType,
NArrayReal,
NArrayInt> &pc,
const Vector<std::string> &real_comp_names,
Expand Down
207 changes: 136 additions & 71 deletions Src/Extern/Conduit/AMReX_Conduit_Blueprint_ParticlesI.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,26 +20,21 @@ namespace amrex
// Note:
// This is a helper function, it's not part of the AMReX Blueprint Interface.
//---------------------------------------------------------------------------//
template <int NStructReal, int NStructInt, int NArrayReal, int NArrayInt>
template <typename ParticleType, int NArrayReal, int NArrayInt>
void
ParticleTileToBlueprint(const ParticleTile<amrex::Particle<NStructReal,
NStructInt>,
ParticleTileToBlueprint(const ParticleTile<ParticleType,
NArrayReal,
NArrayInt> &ptile,
const Vector<std::string> &real_comp_names,
const Vector<std::string> &int_comp_names,
conduit::Node &res,
const std::string &topology_name)
{
int num_particles = ptile.GetArrayOfStructs().size();
int struct_size = sizeof(Particle<NStructReal, NStructInt>);
int num_particles = ptile.size();

// knowing the above, we can zero copy the x,y,z positions + id, cpu
// and any user fields in the AOS

// get the first particle's struct
const auto &pstruct = ptile.GetArrayOfStructs();

// setup a blueprint description for the particle mesh
// create a coordinate set
std::string coordset_name = topology_name + "_coords";
Expand All @@ -63,29 +58,56 @@ ParticleTileToBlueprint(const ParticleTile<amrex::Particle<NStructReal,
//----------------------------------//
// point locations from from aos
//----------------------------------//
char* pbuf = nullptr;

const char* pbuf_const = reinterpret_cast<const char*>(pstruct.data());
char* pbuf = const_cast<char*>(pbuf_const);
if constexpr(ParticleType::is_soa_particle)
{
amrex::ignore_unused(pbuf);

ParticleReal* xp = reinterpret_cast<ParticleReal*>(pbuf); pbuf += sizeof(ParticleReal);
n_coords["values/x"].set_external(xp,
num_particles,
0,
struct_size);
const auto &soa = ptile.GetStructOfArrays();

// for soa entries, we can use standard strides,
// since these are contiguous arrays

n_coords["values/x"].set_external(const_cast<ParticleReal*>(&soa.GetRealData(0)[0]),
num_particles);
#if AMREX_SPACEDIM > 1
ParticleReal* yp = reinterpret_cast<ParticleReal*>(pbuf); pbuf += sizeof(ParticleReal);
n_coords["values/y"].set_external(yp,
num_particles,
0,
struct_size);
n_coords["values/y"].set_external(const_cast<ParticleReal*>(&soa.GetRealData(1)[0]),
num_particles);
#endif
#if AMREX_SPACEDIM > 2
ParticleReal* zp = reinterpret_cast<ParticleReal*>(pbuf); pbuf += sizeof(ParticleReal);
n_coords["values/z"].set_external(zp,
num_particles,
0,
struct_size);
n_coords["values/z"].set_external(const_cast<ParticleReal*>(&soa.GetRealData(2)[0]),
num_particles);
#endif
} else
{
// get the first particle's struct
const auto &pstruct = ptile.GetArrayOfStructs();
const int struct_size = sizeof(ParticleType);

const char* pbuf_const = reinterpret_cast<const char*>(pstruct.data());
pbuf = const_cast<char*>(pbuf_const);

ParticleReal* xp = reinterpret_cast<ParticleReal*>(pbuf); pbuf += sizeof(ParticleReal);
n_coords["values/x"].set_external(xp,
num_particles,
0,
struct_size);
#if AMREX_SPACEDIM > 1
ParticleReal* yp = reinterpret_cast<ParticleReal*>(pbuf); pbuf += sizeof(ParticleReal);
n_coords["values/y"].set_external(yp,
num_particles,
0,
struct_size);
#endif
#if AMREX_SPACEDIM > 2
ParticleReal* zp = reinterpret_cast<ParticleReal*>(pbuf); pbuf += sizeof(ParticleReal);
n_coords["values/z"].set_external(zp,
num_particles,
0,
struct_size);
#endif
}

// fields
conduit::Node &n_fields = res["fields"];
Expand All @@ -95,65 +117,104 @@ ParticleTileToBlueprint(const ParticleTile<amrex::Particle<NStructReal,
// -----------------------------

int vname_real_idx = 0;
// struct real fields, the first set are always the particle positions
// which we wrap above
for (int i = 0; i < NStructReal; i++)
if constexpr(!ParticleType::is_soa_particle)
{
ParticleReal* val = reinterpret_cast<ParticleReal*>(pbuf); pbuf += sizeof(ParticleReal);
conduit::Node &n_f = n_fields[real_comp_names.at(vname_real_idx)];
n_f["topology"] = topology_name;
n_f["association"] = "element";
n_f["values"].set_external(val,
num_particles,
0,
struct_size);
constexpr int struct_size = sizeof(ParticleType);
constexpr int NStructReal = ParticleType::NReal;

vname_real_idx++;
// struct real fields, the first set are always the particle positions
// which we wrap above
for (int i = 0; i < NStructReal; i++)
{
ParticleReal* val = reinterpret_cast<ParticleReal*>(pbuf); pbuf += sizeof(ParticleReal);
conduit::Node &n_f = n_fields[real_comp_names.at(vname_real_idx)];
n_f["topology"] = topology_name;
n_f["association"] = "element";
n_f["values"].set_external(val,
num_particles,
0,
struct_size);

vname_real_idx++;
}
}

//----------------------------------//
// standard integer fields from aos
// (id, cpu)
//----------------------------------//

// id is the first int entry
int* id = reinterpret_cast<int*>(pbuf); pbuf += sizeof(int);
conduit::Node &n_f_id = n_fields[topology_name + "_id"];
if constexpr(!ParticleType::is_soa_particle)
{
const int struct_size = sizeof(ParticleType);

// id is the first int entry
int* id = reinterpret_cast<int*>(pbuf); pbuf += sizeof(int);
conduit::Node &n_f_id = n_fields[topology_name + "_id"];

n_f_id["topology"] = topology_name;
n_f_id["association"] = "element";
n_f_id["values"].set_external(id,
num_particles,
0,
struct_size);

// cpu is the second int entry
int* cpu = reinterpret_cast<int*>(pbuf); pbuf += sizeof(int);
conduit::Node &n_f_cpu = n_fields[topology_name + "_cpu"];

n_f_cpu["topology"] = topology_name;
n_f_cpu["association"] = "element";
n_f_cpu["values"].set_external(cpu,
num_particles,
0,
struct_size);
} else {
const auto &soa = ptile.GetStructOfArrays();

// for soa entries, we can use standard strides,
// since these are contiguous arrays

n_f_id["topology"] = topology_name;
n_f_id["association"] = "element";
n_f_id["values"].set_external(id,
num_particles,
0,
struct_size);
// id is the first int entry
conduit::Node &n_f_id = n_fields[topology_name + "_id"];

// cpu is the second int entry
int* cpu = reinterpret_cast<int*>(pbuf); pbuf += sizeof(int);
conduit::Node &n_f_cpu = n_fields[topology_name + "_cpu"];
n_f_id["topology"] = topology_name;
n_f_id["association"] = "element";
n_f_id["values"].set_external(const_cast<int*>(&soa.GetIntData(0)[0]),
num_particles);

n_f_cpu["topology"] = topology_name;
n_f_cpu["association"] = "element";
n_f_cpu["values"].set_external(cpu,
num_particles,
0,
struct_size);
// cpu is the second int entry
conduit::Node &n_f_cpu = n_fields[topology_name + "_cpu"];

n_f_cpu["topology"] = topology_name;
n_f_cpu["association"] = "element";
n_f_cpu["values"].set_external(const_cast<int*>(&soa.GetIntData(0)[0]),
num_particles);

}

// --------------------------------
// user defined, integer aos fields
// --------------------------------

int vname_int_idx = 0;
for (int i = 0; i < NStructInt; i++)
if constexpr(!ParticleType::is_soa_particle)
{
int* val = reinterpret_cast<int*>(pbuf); pbuf += sizeof(int);
conduit::Node &n_f = n_fields[int_comp_names.at(vname_int_idx)];
n_f["topology"] = topology_name;
n_f["association"] = "element";
n_f["values"].set_external(val,
num_particles,
0,
struct_size);
vname_int_idx++;
constexpr int struct_size = sizeof(ParticleType);
constexpr int NStructInt = ParticleType::NInt;

for (int i = 0; i < NStructInt; i++)
{
int* val = reinterpret_cast<int*>(pbuf); pbuf += sizeof(int);
conduit::Node &n_f = n_fields[int_comp_names.at(vname_int_idx)];
n_f["topology"] = topology_name;
n_f["association"] = "element";
n_f["values"].set_external(val,
num_particles,
0,
struct_size);
vname_int_idx++;
}
}

// -------------------------
Expand Down Expand Up @@ -193,10 +254,9 @@ ParticleTileToBlueprint(const ParticleTile<amrex::Particle<NStructReal,
//---------------------------------------------------------------------------//
// Converts a AMReX Particle Container into a Conduit Mesh Blueprint Hierarchy.
//---------------------------------------------------------------------------//
template <int NStructReal, int NStructInt, int NArrayReal, int NArrayInt>
template <typename ParticleType, int NArrayReal, int NArrayInt>
void
ParticleContainerToBlueprint(const ParticleContainer<NStructReal,
NStructInt,
ParticleContainerToBlueprint(const ParticleContainer_impl<ParticleType,
NArrayReal,
NArrayInt> &pc,
const Vector<std::string> &real_comp_names,
Expand All @@ -209,8 +269,13 @@ ParticleContainerToBlueprint(const ParticleContainer<NStructReal,
// validate varnames, which are used to provide field names
// for user defined aos and soa values.

BL_ASSERT(real_comp_names.size() == (NStructReal + NArrayReal) );
BL_ASSERT(int_comp_names.size() == (NStructInt + NArrayInt) );
if constexpr(ParticleType::is_soa_particle) {
BL_ASSERT(real_comp_names.size() == NArrayReal);
BL_ASSERT(int_comp_names.size() == NArrayInt);
} else {
BL_ASSERT(real_comp_names.size() == (ParticleType::NReal + NArrayReal) );
BL_ASSERT(int_comp_names.size() == (ParticleType::NInt + NArrayInt) );
}

int num_levels = pc.maxLevel() + 1;
int num_domains = 0;
Expand All @@ -224,7 +289,7 @@ ParticleContainerToBlueprint(const ParticleContainer<NStructReal,
int rank = ParallelDescriptor::MyProc();
int nprocs = ParallelDescriptor::NProcs();

using MyParConstIter = ParConstIter<NStructReal, NStructInt, NArrayReal, NArrayInt>;
using MyParConstIter = ParConstIter_impl<ParticleType, NArrayReal, NArrayInt>;

//
// blueprint expects unique ids for each domain published
Expand Down
25 changes: 11 additions & 14 deletions Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,14 +90,12 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)

const int ncomp = sol.nComp();

MF ph = Lp.make(amrlev, mglev, sol.nGrowVect());
MF sh = Lp.make(amrlev, mglev, sol.nGrowVect());
ph.setVal(RT(0.0));
sh.setVal(RT(0.0));
MF p = Lp.make(amrlev, mglev, sol.nGrowVect());
MF r = Lp.make(amrlev, mglev, sol.nGrowVect());
p.setVal(RT(0.0)); // Make sure all entries are initialized to avoid errors
r.setVal(RT(0.0));

MF sorig = Lp.make(amrlev, mglev, nghost);
MF p = Lp.make(amrlev, mglev, nghost);
MF r = Lp.make(amrlev, mglev, nghost);
MF rh = Lp.make(amrlev, mglev, nghost);
MF v = Lp.make(amrlev, mglev, nghost);
MF t = Lp.make(amrlev, mglev, nghost);
Expand Down Expand Up @@ -151,8 +149,7 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
MF::Saxpy(p, -omega, v, 0, 0, ncomp, nghost); // p += -omega*v
MF::Xpay(p, beta, r, 0, 0, ncomp, nghost); // p = r + beta*p
}
ph.LocalCopy(p,0,0,ncomp,nghost);
Lp.apply(amrlev, mglev, v, ph, MLLinOpT<MF>::BCMode::Homogeneous, MLLinOpT<MF>::StateMode::Correction);
Lp.apply(amrlev, mglev, v, p, MLLinOpT<MF>::BCMode::Homogeneous, MLLinOpT<MF>::StateMode::Correction);
Lp.normalize(amrlev, mglev, v);

RT rhTv = dotxy(rh,v);
Expand All @@ -164,9 +161,10 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
{
ret = 2; break;
}
MF::Saxpy(sol, alpha, ph, 0, 0, ncomp, nghost); // sol += alpha * ph
MF::Saxpy(r, -alpha, v, 0, 0, ncomp, nghost); // r += -alpha * v
MF::Saxpy(sol, alpha, p, 0, 0, ncomp, nghost); // sol += alpha * p
MF::Saxpy(r, -alpha, v, 0, 0, ncomp, nghost); // r += -alpha * v

rnorm = norm_inf(r);
rnorm = norm_inf(r);

if ( verbose > 2 && ParallelDescriptor::IOProcessor() )
Expand All @@ -179,8 +177,7 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)

if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) { break; }

sh.LocalCopy(r,0,0,ncomp,nghost);
Lp.apply(amrlev, mglev, t, sh, MLLinOpT<MF>::BCMode::Homogeneous, MLLinOpT<MF>::StateMode::Correction);
Lp.apply(amrlev, mglev, t, r, MLLinOpT<MF>::BCMode::Homogeneous, MLLinOpT<MF>::StateMode::Correction);
Lp.normalize(amrlev, mglev, t);
//
// This is a little funky. I want to elide one of the reductions
Expand All @@ -201,8 +198,8 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
{
ret = 3; break;
}
MF::Saxpy(sol, omega, sh, 0, 0, ncomp, nghost); // sol += omega * sh
MF::Saxpy(r, -omega, t, 0, 0, ncomp, nghost); // r += -omega * t
MF::Saxpy(sol, omega, r, 0, 0, ncomp, nghost); // sol += omega * r
MF::Saxpy(r, -omega, t, 0, 0, ncomp, nghost); // r += -omega * t

rnorm = norm_inf(r);

Expand Down
13 changes: 13 additions & 0 deletions Tests/Particles/Ascent_Insitu_SOA/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
if ( NOT AMReX_ASCENT )
return ()
endif ()

foreach(D IN LISTS AMReX_SPACEDIM)
set(_sources main.cpp)
set(_input_files inputs.rt )

setup_test(${D} _sources _input_files NTASKS 2)

unset(_sources)
unset(_input_files)
endforeach()
Loading

0 comments on commit 52a6f70

Please sign in to comment.