Skip to content
Permalink
Browse files

Merge #3144

3144: external force density in first integration step r=RudolfWeeber a=KaiSzuttor

Fixes #2875

Description of changes:
 - initialize the force density member with the external force density set by the user
 - if node velocity is set from the interface, the force density on this node is reset

Co-authored-by: Kai Szuttor <kai@icp.uni-stuttgart.de>
Co-authored-by: RudolfWeeber <weeber@icp.uni-stuttgart.de>
  • Loading branch information...
3 people committed Sep 11, 2019
2 parents 5c69842 + 583da85 commit ea237af5d31d385aaa01f4332369722529975407
@@ -186,7 +186,7 @@ def main():
mass = np.array([0.17])

num_types = len(wca_sig)

def mix_eps(eps1, eps2, rule='LB'):
return math.sqrt(eps1 * eps2)

@@ -172,6 +172,51 @@ static double fluidstep = 0.0;
#include "grid.hpp"

/********************** The Main LB Part *************************************/

/**
* @brief Initialize fluid nodes.
* @param[out] fields Vector containing the fluid nodes
* @param[in] lb_parameters Parameters for the LB
* @param[in] lb_lattice Lattice instance
*/
void lb_initialize_fields(std::vector<LB_FluidNode> &fields,
LB_Parameters const &lb_parameters,
Lattice const &lb_lattice) {
fields.resize(lb_lattice.halo_grid_volume);
for (auto &field : fields) {
field.force_density = lb_parameters.ext_force_density;
#ifdef LB_BOUNDARIES
field.boundary = false;
#endif // LB_BOUNDARIES
}
}

/** (Re-)allocate memory for the fluid and initialize pointers. */
void lb_realloc_fluid(LB_FluidData &lb_fluid_a, LB_FluidData &lb_fluid_b,
const Lattice::index_t halo_grid_volume,
LB_Fluid &lb_fluid, LB_Fluid &lb_fluid_post) {
const std::array<int, 2> size = {{D3Q19::n_vel, halo_grid_volume}};

lb_fluid_a.resize(size);
lb_fluid_b.resize(size);

using Utils::Span;
for (int i = 0; i < size[0]; i++) {
lb_fluid[i] = Span<double>(lb_fluid_a[i].origin(), size[1]);
lb_fluid_post[i] = Span<double>(lb_fluid_b[i].origin(), size[1]);
}
}

void lb_set_equilibrium_populations(const Lattice &lb_lattice,
const LB_Parameters &lb_parameters) {
for (Lattice::index_t index = 0; index < lb_lattice.halo_grid_volume;
++index) {
lb_set_population_from_density_momentum_density_stress(
index, lb_parameters.density, Utils::Vector3d{} /*momentum density*/,
Utils::Vector6d{} /*stress*/);
}
}

void lb_init(const LB_Parameters &lb_parameters) {
if (lb_parameters.agrid <= 0.0) {
runtimeErrorMsg()
@@ -181,7 +226,6 @@ void lb_init(const LB_Parameters &lb_parameters) {
return;

/* initialize the local lattice domain */

try {
lblattice = Lattice(lb_parameters.agrid, 0.5 /*offset*/, 1 /*halo size*/,
local_geo.length(), local_geo.my_right(),
@@ -193,40 +237,28 @@ void lb_init(const LB_Parameters &lb_parameters) {

/* allocate memory for data structures */
lb_realloc_fluid(lbfluid_a, lbfluid_b, lblattice.halo_grid_volume, lbfluid,
lbfluid_post, lbfields);
lbfluid_post);

lb_initialize_fields(lbfields, lbpar, lblattice);

/* prepare the halo communication */
lb_prepare_communication(update_halo_comm, lblattice);

/* initialize derived parameters */
lb_reinit_parameters(lbpar);

/* setup the initial populations */
lb_reinit_fluid(lbfields, lblattice, lbpar);
}

void lb_reinit_fluid(std::vector<LB_FluidNode> &lb_fields,
const Lattice &lb_lattice,
const LB_Parameters &lb_parameters) {
std::fill(lb_fields.begin(), lb_fields.end(), LB_FluidNode());
/* default values for fields in lattice units */
Utils::Vector3d momentum_density{};
Utils::Vector6d stress{};

for (Lattice::index_t index = 0; index < lb_lattice.halo_grid_volume;
++index) {
// sets equilibrium distribution
lb_set_population_from_density_momentum_density_stress(
index, lb_parameters.density, momentum_density, stress);

#ifdef LB_BOUNDARIES
lb_fields[index].boundary = 0;
#endif // LB_BOUNDARIES
}
lb_set_equilibrium_populations(lblattice, lbpar);

#ifdef LB_BOUNDARIES
LBBoundaries::lb_init_boundaries();
#endif // LB_BOUNDARIES
#endif
}

void lb_reinit_fluid(std::vector<LB_FluidNode> &lb_fields,
Lattice const &lb_lattice,
LB_Parameters const &lb_parameters) {
lb_set_equilibrium_populations(lb_lattice, lb_parameters);
lb_initialize_fields(lbfields, lb_parameters, lb_lattice);
}

void lb_reinit_parameters(LB_Parameters &lb_parameters) {
@@ -587,25 +619,6 @@ void lb_fluid_set_rng_state(uint64_t counter) {

/***********************************************************************/

/** (Re-)allocate memory for the fluid and initialize pointers. */
void lb_realloc_fluid(LB_FluidData &lb_fluid_a, LB_FluidData &lb_fluid_b,
const Lattice::index_t halo_grid_volume,
LB_Fluid &lb_fluid, LB_Fluid &lb_fluid_post,
std::vector<LB_FluidNode> &lb_fields) {
const std::array<int, 2> size = {{D3Q19::n_vel, halo_grid_volume}};

lb_fluid_a.resize(size);
lb_fluid_b.resize(size);

using Utils::Span;
for (int i = 0; i < size[0]; i++) {
lb_fluid[i] = Span<double>(lb_fluid_a[i].origin(), size[1]);
lb_fluid_post[i] = Span<double>(lb_fluid_b[i].origin(), size[1]);
}

lb_fields.resize(halo_grid_volume);
}

/** Set up the structures for exchange of the halo regions.
* See also \ref halo.cpp
*/
@@ -938,18 +951,6 @@ inline void lb_collide_stream() {
}
#endif // LB_BOUNDARIES

#ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
// Safeguard the node forces so that we can later use them for the IBM
// particle update
// In the following loop the lbfields[XX].force are reset to zero
// Safeguard the node forces so that we can later use them for the IBM
// particle update In the following loop the lbfields[XX].force are reset to
// zero
for (int i = 0; i < lblattice.halo_grid_volume; ++i) {
lbfields[i].force_density_buf = lbfields[i].force_density;
}
#endif

Lattice::index_t index = lblattice.halo_offset;
for (int z = 1; z <= lblattice.grid[2]; z++) {
for (int y = 1; y <= lblattice.grid[1]; y++) {
@@ -974,6 +975,12 @@ inline void lb_collide_stream() {
auto const modes_with_forces =
lb_apply_forces(index, thermalized_modes, lbpar, lbfields);

#ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
// Safeguard the node forces so that we can later use them for the IBM
// particle update
lbfields[index].force_density_buf = lbfields[index].force_density;
#endif

/* reset the force density */
lbfields[index].force_density = lbpar.ext_force_density;

@@ -143,13 +143,6 @@ extern Lattice lblattice;

extern HaloCommunicator update_halo_comm;

void lb_realloc_fluid(boost::multi_array<double, 2> &lb_fluid_a,
boost::multi_array<double, 2> &lb_fluid_b,
Lattice::index_t halo_grid_volume,
std::array<Utils::Span<double>, 19> &lb_fluid,
std::array<Utils::Span<double>, 19> &lb_fluid_post,
std::vector<LB_FluidNode> &lb_fields);

void lb_init(const LB_Parameters &lb_parameters);

void lb_reinit_fluid(std::vector<LB_FluidNode> &lb_fields,
@@ -287,6 +280,9 @@ void lb_calc_fluid_momentum(double *result, const LB_Parameters &lb_parameters,
const std::vector<LB_FluidNode> &lb_fields,
const Lattice &lb_lattice);
void lb_collect_boundary_forces(double *result);
void lb_initialize_fields(std::vector<LB_FluidNode> &fields,
LB_Parameters const &lb_parameters,
Lattice const &lb_lattice);

/*@}*/

@@ -89,6 +89,17 @@ void mpi_lb_set_population(Utils::Vector3i const &index,

REGISTER_CALLBACK(mpi_lb_set_population)

void mpi_lb_set_force_density(Utils::Vector3i const &index,
Utils::Vector3d const &force_density) {
lb_set(index, [&](auto index) {
auto const linear_index =
get_linear_index(lblattice.local_index(index), lblattice.halo_grid);
lbfields[linear_index].force_density = force_density;
});
}

REGISTER_CALLBACK(mpi_lb_set_force_density)

auto mpi_lb_get_momentum_density(Utils::Vector3i const &index) {
return lb_calc_fluid_kernel(index, [&](auto modes, auto force_density) {
return lb_calc_momentum_density(modes, force_density);
@@ -1265,6 +1276,7 @@ void lb_lbnode_set_velocity(const Utils::Vector3i &ind,
lb_get_population_from_density_momentum_density_stress(
density, momentum_density, stress);
mpi_call_all(mpi_lb_set_population, ind, population);
mpi_call_all(mpi_lb_set_force_density, ind, Utils::Vector3d{});
}
}

@@ -1308,6 +1320,7 @@ void lb_lbfluid_on_lb_params_change(LBParam field) {
break;
case LBParam::VISCOSITY:
case LBParam::EXT_FORCE_DENSITY:
lb_initialize_fields(lbfields, lbpar, lblattice);
case LBParam::BULKVISC:
case LBParam::KT:
case LBParam::GAMMA_ODD:
@@ -213,10 +213,6 @@ void lb_init_boundaries_GPU(int n_lb_boundaries, int number_of_boundnodes,
int *host_boundary_index_list,
float *lb_bounday_velocity);
#endif
void lb_init_extern_nodeforcedensities_GPU(
int n_extern_node_force_densities,
LB_extern_nodeforcedensity_gpu *host_extern_node_force_densities,
LB_parameters_gpu *lbpar_gpu);

void lb_set_agrid_gpu(double agrid);

0 comments on commit ea237af

Please sign in to comment.
You can’t perform that action at this time.