From 9aa17adfb7ea800c9f3e996f963cb1a64c016610 Mon Sep 17 00:00:00 2001 From: Kai Szuttor Date: Tue, 18 Feb 2020 16:38:45 +0100 Subject: [PATCH] lbboundaries: initialization logic for gpu. --- .../grid_based_algorithms/lb_boundaries.cpp | 183 +++++++++--------- 1 file changed, 95 insertions(+), 88 deletions(-) diff --git a/src/core/grid_based_algorithms/lb_boundaries.cpp b/src/core/grid_based_algorithms/lb_boundaries.cpp index 027c5fb889c..fcfb1b59dd9 100644 --- a/src/core/grid_based_algorithms/lb_boundaries.cpp +++ b/src/core/grid_based_algorithms/lb_boundaries.cpp @@ -25,6 +25,7 @@ * */ + #include "grid_based_algorithms/lb_boundaries.hpp" #include "communication.hpp" @@ -38,7 +39,9 @@ #include "lbboundaries/LBBoundary.hpp" #include + using Utils::get_linear_index; + #include #include @@ -65,53 +68,101 @@ void remove(const std::shared_ptr &b) { on_lbboundary_change(); } -/** Initialize boundary conditions for all constraints in the system. */ -void lb_init_boundaries() { - if (lattice_switch == ActiveLB::GPU) { -#if defined(CUDA) && defined(LB_BOUNDARIES_GPU) - if (this_node != 0) { - return; - } +void ek_init_boundaries() { +#if defined(CUDA) && defined(EK_BOUNDARIES) + int number_of_boundnodes = 0; - int number_of_boundnodes = 0; - std::vector host_boundary_node_list; - std::vector host_boundary_index_list; - size_t size_of_index; + std::vector host_wallcharge_species_density; + float node_wallcharge = 0.0f; + int wallcharge_species = -1, charged_boundaries = 0; + bool node_charged = false; -#ifdef EK_BOUNDARIES - std::vector host_wallcharge_species_density; - float node_wallcharge = 0.0f; - int wallcharge_species = -1, charged_boundaries = 0; - bool node_charged = false; + for (auto &lbboundarie : lbboundaries) { + (*lbboundarie).set_net_charge(0.0); + } + if (ek_initialized) { + host_wallcharge_species_density.resize(ek_parameters.number_of_nodes); for (auto &lbboundarie : lbboundaries) { - (*lbboundarie).set_net_charge(0.0); + if ((*lbboundarie).charge_density() != 0.0) { + charged_boundaries = 1; + break; + } } - if (ek_initialized) { - host_wallcharge_species_density.resize(ek_parameters.number_of_nodes); - for (auto &lbboundarie : lbboundaries) { - if ((*lbboundarie).charge_density() != 0.0) { - charged_boundaries = 1; - break; - } + for (int n = 0; n < int(ek_parameters.number_of_species); n++) + if (ek_parameters.valency[n] != 0.0) { + wallcharge_species = n; + break; } - for (int n = 0; n < int(ek_parameters.number_of_species); n++) - if (ek_parameters.valency[n] != 0.0) { - wallcharge_species = n; - break; - } + ek_gather_wallcharge_species_density(host_wallcharge_species_density.data(), + wallcharge_species); - ek_gather_wallcharge_species_density( - host_wallcharge_species_density.data(), wallcharge_species); + if (wallcharge_species == -1 && charged_boundaries) { + runtimeErrorMsg() + << "no charged species available to create wall charge\n"; + } - if (wallcharge_species == -1 && charged_boundaries) { - runtimeErrorMsg() - << "no charged species available to create wall charge\n"; + for (int z = 0; z < int(lbpar_gpu.dim_z); z++) { + for (int y = 0; y < int(lbpar_gpu.dim_y); y++) { + for (int x = 0; x < int(lbpar_gpu.dim_x); x++) { + auto const pos = static_cast(lbpar_gpu.agrid) * + (Utils::Vector3d{1. * x, 1. * y, 1. * z} + + Utils::Vector3d::broadcast(0.5)); + node_charged = false; + node_wallcharge = 0.0f; + + std::vector> boundaries; + std::copy_if(lbboundaries.begin(), lbboundaries.end(), + std::back_inserter(boundaries), [&pos](auto const lbb) { + return not lbb->shape().is_inside(pos); + }); + for (auto lbb : boundaries) { + if ((*lbb).charge_density() != 0.0f) { + node_charged = true; + node_wallcharge += (*lbb).charge_density() * ek_parameters.agrid * + ek_parameters.agrid * ek_parameters.agrid; + (*lbb).set_net_charge( + (*lbb).net_charge() + + (*lbb).charge_density() * ek_parameters.agrid * + ek_parameters.agrid * ek_parameters.agrid); + } + } + if (not boundaries.empty()) { + number_of_boundnodes++; + } + ek_parameters.number_of_boundary_nodes = number_of_boundnodes; + + if (wallcharge_species != -1) { + if (node_charged) + host_wallcharge_species_density[ek_parameters.dim_y * + ek_parameters.dim_x * z + + ek_parameters.dim_x * y + x] = + node_wallcharge / ek_parameters.valency[wallcharge_species]; + } + } } } + ek_init_species_density_wallcharge(host_wallcharge_species_density.data(), + wallcharge_species); + } #endif +} + +/** Initialize boundary conditions for all constraints in the system. */ +void lb_init_boundaries() { + if (lattice_switch == ActiveLB::GPU) { +#if defined(CUDA) && defined(LB_BOUNDARIES_GPU) + if (this_node != 0) { + return; + } + ek_init_boundaries(); + int number_of_boundnodes = 0; + std::vector host_boundary_node_list; + std::vector host_boundary_index_list; + size_t size_of_index; + for (int z = 0; z < int(lbpar_gpu.dim_z); z++) { for (int y = 0; y < int(lbpar_gpu.dim_y); y++) { for (int x = 0; x < int(lbpar_gpu.dim_x); x++) { @@ -119,66 +170,29 @@ void lb_init_boundaries() { (Utils::Vector3d{1. * x, 1. * y, 1. * z} + Utils::Vector3d::broadcast(0.5)); -#ifdef EK_BOUNDARIES - if (ek_initialized) { - node_charged = false; - node_wallcharge = 0.0f; - } -#endif - std::vector inside_boundaries; - for (auto lbb = lbboundaries.begin(); lbb != lbboundaries.end(); - ++lbb) { - auto const is_boundary = not(**lbb).shape().is_inside(pos); - inside_boundaries.push_back(is_boundary); -#ifdef EK_BOUNDARIES - if (ek_initialized) { - if (is_boundary && (**lbb).charge_density() != 0.0f) { - node_charged = true; - node_wallcharge += (**lbb).charge_density() * - ek_parameters.agrid * ek_parameters.agrid * - ek_parameters.agrid; - (**lbb).set_net_charge( - (**lbb).net_charge() + - (**lbb).charge_density() * ek_parameters.agrid * - ek_parameters.agrid * ek_parameters.agrid); - } - } -#endif - } - auto const it = - std::find_if(inside_boundaries.begin(), inside_boundaries.end(), - [](bool b) { return b; }); - if (it != inside_boundaries.end()) { + // take last boundary containing the node + auto const boundary = + std::find_if(lbboundaries.rbegin(), lbboundaries.rend(), + [&pos](auto const lbb) { + return not lbb->shape().is_inside(pos); + }); + + if (boundary != lbboundaries.rend()) { size_of_index = (number_of_boundnodes + 1) * sizeof(int); host_boundary_node_list.resize(size_of_index); host_boundary_index_list.resize(size_of_index); host_boundary_node_list[number_of_boundnodes] = x + lbpar_gpu.dim_x * y + lbpar_gpu.dim_x * lbpar_gpu.dim_y * z; auto const boundary_number = - std::distance(inside_boundaries.begin(), it); + std::distance(boundary, lbboundaries.rend()); host_boundary_index_list[number_of_boundnodes] = boundary_number + 1; number_of_boundnodes++; } - lbpar_gpu.number_of_boundnodes = number_of_boundnodes; - -#ifdef EK_BOUNDARIES - if (ek_initialized) { - ek_parameters.number_of_boundary_nodes = number_of_boundnodes; - - if (wallcharge_species != -1) { - if (node_charged) - host_wallcharge_species_density[ek_parameters.dim_y * - ek_parameters.dim_x * z + - ek_parameters.dim_x * y + x] = - node_wallcharge / ek_parameters.valency[wallcharge_species]; - } - } -#endif } } } - + lbpar_gpu.number_of_boundnodes = number_of_boundnodes; /**call of cuda fkt*/ std::vector boundary_velocity(3 * (lbboundaries.size() + 1)); int n = 0; @@ -198,13 +212,6 @@ void lb_init_boundaries() { host_boundary_index_list.data(), boundary_velocity.data()); -#ifdef EK_BOUNDARIES - if (ek_initialized) { - ek_init_species_density_wallcharge(host_wallcharge_species_density.data(), - wallcharge_species); - } -#endif - #endif /* defined ( CUDA) && defined (LB_BOUNDARIES_GPU) */ } else if (lattice_switch == ActiveLB::CPU) { #if defined(LB_BOUNDARIES)