diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc index a6f47db12..4c9d37833 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc @@ -1,313 +1,313 @@ /** * @file material_cohesive_linear_inline_impl.cc * * @author Mauro Corrado * @author Marco Vocialta * * @date creation: Wed Apr 22 2015 * @date last modification: Thu Jan 14 2016 * * @brief Inline functions of the MaterialCohesiveLinear * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template inline Real MaterialCohesiveLinear::computeEffectiveNorm(const Matrix & stress, const Vector & normal, const Vector & tangent, Vector & normal_traction) const { normal_traction.mul(stress, normal); Real normal_contrib = normal_traction.dot(normal); /// in 3D tangential components must be summed Real tangent_contrib = 0; if (dim == 2) { Real tangent_contrib_tmp = normal_traction.dot(tangent); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } else if (dim == 3) { for (UInt s = 0; s < dim - 1; ++s) { const Vector tangent_v(tangent.storage() + s * dim, dim); Real tangent_contrib_tmp = normal_traction.dot(tangent_v); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } } tangent_contrib = std::sqrt(tangent_contrib); normal_contrib = std::max(Real(0.), normal_contrib); return std::sqrt(normal_contrib * normal_contrib + tangent_contrib * tangent_contrib * beta2_inv); } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeTractionOnQuad( Vector & traction, Real & delta_max, __attribute__((unused)) const Real & delta_max_prev, const Real & delta_c, const Vector & insertion_stress, const Real & sigma_c, Vector & opening, Vector & opening_prec, const Vector & normal, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, bool & reduction_penalty, Real & current_penalty, Vector & contact_traction, Vector & contact_opening) { /// compute normal and tangential opening vectors normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; tangential_opening = opening; tangential_opening -= normal_opening; tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; penetration = normal_opening_norm < -Math::getTolerance(); if (this->contact_after_breaking == false && Math::are_float_equal(damage, 1.)) penetration = false; /** * if during the convergence loop a cohesive element continues to * jumps from penetration to opening, and convergence is not * reached, its penalty parameter will be reduced in the * recomputation of the same incremental step. recompute is set * equal to true when convergence is not reached in the * solveStepCohesive function and the execution of the program * goes back to the main file where the variable load_reduction * is set equal to true. */ Real normal_opening_prec_norm = opening_prec.dot(normal); opening_prec = opening; if (!this->model->isExplicit() && !this->recompute) if ((normal_opening_prec_norm * normal_opening_norm) < 0) { reduction_penalty = true; } if (penetration) { if (this->recompute && reduction_penalty){ /// the penalty parameter is locally reduced - current_penalty = this->penalty / 1000.; + current_penalty = this->penalty / 100.; } else current_penalty = this->penalty; /// use penalty coefficient in case of penetration contact_traction = normal_opening; contact_traction *= current_penalty; contact_opening = normal_opening; /// don't consider penetration contribution for delta opening = tangential_opening; normal_opening.clear(); } else { delta += normal_opening_norm * normal_opening_norm; contact_traction.clear(); contact_opening.clear(); } delta = std::sqrt(delta); /// update maximum displacement and damage delta_max = std::max(delta_max, delta); damage = std::min(delta_max / delta_c, Real(1.)); /** * Compute traction @f$ \mathbf{T} = \left( * \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n * \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1- * \frac{\delta}{\delta_c} \right)@f$ */ if (Math::are_float_equal(damage, 1.)) traction.clear(); else if (Math::are_float_equal(damage, 0.)) { if (penetration) traction.clear(); else traction = insertion_stress; } else { traction = tangential_opening; traction *= this->beta2_kappa; traction += normal_opening; AKANTU_DEBUG_ASSERT(delta_max != 0., "Division by zero, tolerance might be too low"); traction *= sigma_c / delta_max * (1. - damage); } } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeTangentTractionOnQuad( Matrix & tangent, Real & delta_max, const Real & delta_c, const Real & sigma_c, Vector & opening, const Vector & normal, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, bool & reduction_penalty, Real & current_penalty, Vector & contact_opening) { /// During the update of the residual the interpenetrations are /// stored in the array "contact_opening", therefore, in the case /// of penetration, in the array "opening" there are only the /// tangential components. opening += contact_opening; /// compute normal and tangential opening vectors normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; tangential_opening = opening; tangential_opening -= normal_opening; tangential_opening_norm = tangential_opening.norm(); penetration = normal_opening_norm < -Math::getTolerance(); if (contact_after_breaking == false && Math::are_float_equal(damage, 1.)) penetration = false; Real derivative = 0; // derivative = d(t/delta)/ddelta Real t = 0; Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(normal, normal); if (penetration){ if (recompute && reduction_penalty) - current_penalty = this->penalty / 1000.; + current_penalty = this->penalty / 100.; else current_penalty = this->penalty; /// stiffness in compression given by the penalty parameter tangent += n_outer_n; tangent *= current_penalty; opening = tangential_opening; normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; } else{ delta += normal_opening_norm * normal_opening_norm; } delta = std::sqrt(delta); /// Delta has to be different from 0 to have finite values of /// tangential stiffness. At the element insertion, delta = /// 0. Therefore, a fictictious value is defined, for the /// evaluation of the first value of K. if (delta < Math::getTolerance()) delta = (delta_c)/1000.; if (delta >= delta_max){ if (delta <= delta_c){ derivative = -sigma_c / (delta * delta); t = sigma_c * (1 - delta / delta_c); } else { derivative = 0.; t = 0.; } } else if (delta < delta_max){ Real tmax = sigma_c * (1 - delta_max / delta_c); t = tmax / delta_max * delta; } /// computation of the derivative of the constitutive law (dT/ddelta) Matrix I(spatial_dimension, spatial_dimension); I.eye(this->beta2_kappa); Matrix nn(n_outer_n); nn *= (1. - this->beta2_kappa); nn += I; nn *= t/delta; Vector t_tilde(normal_opening); t_tilde *= (1. - this->beta2_kappa2); Vector mm(opening); mm *= this->beta2_kappa2; t_tilde += mm; Vector t_hat(normal_opening); t_hat += this->beta2_kappa * tangential_opening; Matrix prov(spatial_dimension, spatial_dimension); prov.outerProduct(t_hat, t_tilde); prov *= derivative/delta; prov += nn; Matrix prov_t = prov.transpose(); tangent += prov_t; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ /* -------------------------------------------------------------------------- */ #endif //__AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_inline_impl.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_inline_impl.cc index 8775a3a89..69fcef28d 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_inline_impl.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_inline_impl.cc @@ -1,318 +1,318 @@ /** * @file solid_mechanics_model_cohesive_inline_impl.cc * * @author Mauro Corrado * @author Marco Vocialta * * @date creation: Fri Jan 18 2013 * @date last modification: Thu Jan 14 2016 * * @brief Implementation of inline functions for the Cohesive element model * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include #include "material_cohesive.hh" #ifndef __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_INLINE_IMPL_CC__ #define __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_INLINE_IMPL_CC__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template bool SolidMechanicsModelCohesive::solveStepCohesive(Real tolerance, Real & error, UInt max_iteration, bool load_reduction, Real tol_increase_factor, bool do_not_factorize) { EventManager::sendEvent(SolidMechanicsModelEvent::BeforeSolveStepEvent(method)); this->implicitPred(); bool insertion_new_element = true; bool converged = false; Array * displacement_tmp = NULL; Array * velocity_tmp = NULL; Array * acceleration_tmp = NULL; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); /// Loop for the insertion of new cohesive elements while (insertion_new_element) { if (is_extrinsic) { /** * If in extrinsic the solution of the previous incremental step * is saved in temporary arrays created for displacements, * velocities and accelerations. Such arrays are used to find * the solution with the Newton-Raphson scheme (this is done by * pointing the pointer "displacement" to displacement_tmp). In * this way, inside the array "displacement" is kept the * solution of the previous incremental step, and in * "displacement_tmp" is saved the current solution. */ Array * tmp_swap; if(!displacement_tmp) { displacement_tmp = new Array(*(this->displacement)); } else { displacement_tmp->resize(this->displacement->getSize()); displacement_tmp->copy(*(this->displacement)); } tmp_swap = displacement_tmp; displacement_tmp = this->displacement; this->displacement = tmp_swap; if(!velocity_tmp) { velocity_tmp = new Array(*(this->velocity)); } else { velocity_tmp->resize(this->velocity->getSize()); velocity_tmp->copy(*(this->velocity)); } tmp_swap = velocity_tmp; velocity_tmp = this->velocity; this->velocity = tmp_swap; if(!acceleration_tmp) { acceleration_tmp = new Array(*(this->acceleration)); } else { acceleration_tmp->resize(this->acceleration->getSize()); acceleration_tmp->copy(*(this->acceleration)); } tmp_swap = acceleration_tmp; acceleration_tmp = this->acceleration; this->acceleration = tmp_swap; } this->updateResidual(); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); bool need_factorize = !do_not_factorize; if (method ==_implicit_dynamic) { AKANTU_DEBUG_ASSERT(mass_matrix != NULL, "You should first initialize the implicit solver and assemble the mass matrix"); } switch (cmethod) { case _scm_newton_raphson_tangent: case _scm_newton_raphson_tangent_not_computed: break; case _scm_newton_raphson_tangent_modified: this->assembleStiffnessMatrix(); break; default: AKANTU_DEBUG_ERROR("The resolution method " << cmethod << " has not been implemented!"); } UInt iter = 0; converged = false; error = 0.; if(criteria == _scc_residual) { converged = this->testConvergence (tolerance, error); if(converged) return converged; } /// Loop to solve the nonlinear system do { if (cmethod == _scm_newton_raphson_tangent) this->assembleStiffnessMatrix(); solve (*increment, 1., need_factorize); this->implicitCorr(); this->updateResidual(); converged = this->testConvergence (tolerance, error); iter++; AKANTU_DEBUG_INFO("[" << criteria << "] Convergence iteration " << std::setw(std::log10(max_iteration)) << iter << ": error " << error << (converged ? " < " : " > ") << tolerance); switch (cmethod) { case _scm_newton_raphson_tangent: need_factorize = true; break; case _scm_newton_raphson_tangent_not_computed: case _scm_newton_raphson_tangent_modified: need_factorize = false; break; default: AKANTU_DEBUG_ERROR("The resolution method " << cmethod << " has not been implemented!"); } } while (!converged && iter < max_iteration); /** * This is to save the obtained result and proceed with the * simulation even if the error is higher than the pre-fixed * tolerance. This is done only after loading reduction * (load_reduction = true). */ if (load_reduction && (error < tolerance * tol_increase_factor)) converged = true; if (converged) { // EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method)); // !!! add sendEvent to call computeCauchyStress !!!! if (prank==0){ std::cout << "Error after convergence: " << error << std::endl; std::cout << "no. of iterations: " << iter << std::endl; } } else if(iter == max_iteration) { if (prank==0){ AKANTU_DEBUG_WARNING("[" << criteria << "] Convergence not reached after " << std::setw(std::log10(max_iteration)) << iter << " iteration" << (iter == 1 ? "" : "s") << "!" << std::endl); std::cout << "Error after NON convergence: " << error << std::endl; } } if (is_extrinsic) { /** * If is extrinsic the pointer "displacement" is moved back to * the array displacement. In this way, the array displacement is * correctly resized during the checkCohesiveStress function (in * case new cohesive elements are added). This is possible * because the procedure called by checkCohesiveStress does not * use the displacement field (the correct one is now stored in * displacement_tmp), but directly the stress field that is * already computed. */ Array * tmp_swap; tmp_swap = displacement_tmp; displacement_tmp = this->displacement; this->displacement = tmp_swap; tmp_swap = velocity_tmp; velocity_tmp = this->velocity; this->velocity = tmp_swap; tmp_swap = acceleration_tmp; acceleration_tmp = this->acceleration; this->acceleration = tmp_swap; /// If convergence is reached, call checkCohesiveStress in order /// to check if cohesive elements have to be introduced if (converged){ UInt new_cohesive_elements = checkCohesiveStress(); if(new_cohesive_elements == 0){ insertion_new_element = false; } else { insertion_new_element = true; if (prank==0) std::cout << "No. cohesive elements inserted = " << new_cohesive_elements << std::endl; } } } if (!converged && load_reduction) insertion_new_element = false; /** * If convergence is not reached, there is the possibility to * return back to the main file and reduce the load. Before doing * this, a pre-fixed value as to be defined for the parameter * delta_max of the cohesive elements introduced in the current * incremental step. This is done by calling the function * checkDeltaMax. */ if (!converged) { insertion_new_element = false; for (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat = dynamic_cast(*materials[m]); mat.checkDeltaMax(_not_ghost); } catch(std::bad_cast&) { } } } } //end loop for the insertion of new cohesive elements /** * When the solution to the current incremental step is computed (no * more cohesive elements have to be introduced), call the function * to compute the energies. */ if ((is_extrinsic && converged)) { for(UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat = dynamic_cast(*materials[m]); mat.computeEnergies(); } catch (std::bad_cast & bce) { } } EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method)); /** * The function resetVariables is necessary to correctly set a * variable that permit to decrease locally the penalty parameter * for compression. */ for (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat = dynamic_cast(*materials[m]); mat.resetVariables(_not_ghost); } catch(std::bad_cast&) { } } /// The correct solution is saved this->displacement->copy(*displacement_tmp); this->velocity ->copy(*velocity_tmp); this->acceleration->copy(*acceleration_tmp); - delete displacement_tmp; - delete velocity_tmp; - delete acceleration_tmp; - } + delete displacement_tmp; + delete velocity_tmp; + delete acceleration_tmp; + return insertion_new_element; } __END_AKANTU__ #if defined (AKANTU_PARALLEL_COHESIVE_ELEMENT) # include "solid_mechanics_model_cohesive_parallel_inline_impl.cc" #endif #endif /* __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_INLINE_IMPL_CC__ */ diff --git a/src/synchronizer/static_communicator_mpi.cc b/src/synchronizer/static_communicator_mpi.cc index 2c2d6cbaa..a79aa0525 100644 --- a/src/synchronizer/static_communicator_mpi.cc +++ b/src/synchronizer/static_communicator_mpi.cc @@ -1,502 +1,503 @@ /** * @file static_communicator_mpi.cc * * @author Nicolas Richart * * @date creation: Sun Sep 26 2010 * @date last modification: Thu Jan 21 2016 * * @brief StaticCommunicatorMPI implementation * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "static_communicator_mpi.hh" #include "mpi_type_wrapper.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ MPI_Op MPITypeWrapper::synchronizer_operation_to_mpi_op[_so_null + 1] = { MPI_SUM, MPI_MIN, MPI_MAX, MPI_PROD, MPI_LAND, MPI_BAND, MPI_LOR, MPI_BOR, MPI_LXOR, MPI_BXOR, MPI_MINLOC, MPI_MAXLOC, MPI_OP_NULL}; class CommunicationRequestMPI : public CommunicationRequest { public: CommunicationRequestMPI(UInt source, UInt dest); ~CommunicationRequestMPI(); MPI_Request * getMPIRequest() { return request; }; private: MPI_Request * request; }; /* -------------------------------------------------------------------------- */ /* Implementation */ /* -------------------------------------------------------------------------- */ CommunicationRequestMPI::CommunicationRequestMPI(UInt source, UInt dest) : CommunicationRequest(source, dest) { request = new MPI_Request; } /* -------------------------------------------------------------------------- */ CommunicationRequestMPI::~CommunicationRequestMPI() { delete request; } /* -------------------------------------------------------------------------- */ StaticCommunicatorMPI::StaticCommunicatorMPI(int & argc, char **& argv) : RealStaticCommunicator(argc, argv) { int is_initialized = false; MPI_Initialized(&is_initialized); if (!is_initialized) { MPI_Init(&argc, &argv); } this->is_externaly_initialized = is_initialized; mpi_data = new MPITypeWrapper(*this); mpi_data->setMPICommunicator(MPI_COMM_WORLD); } /* -------------------------------------------------------------------------- */ StaticCommunicatorMPI::~StaticCommunicatorMPI() { if (!this->is_externaly_initialized) { MPI_Finalize(); + delete this->mpi_data; } } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::send(T * buffer, Int size, Int receiver, Int tag) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Send(buffer, size, type, receiver, tag, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Send."); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::receive(T * buffer, Int size, Int sender, Int tag) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Status status; MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Recv(buffer, size, type, sender, tag, communicator, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Recv."); } /* -------------------------------------------------------------------------- */ template CommunicationRequest * StaticCommunicatorMPI::asyncSend(T * buffer, Int size, Int receiver, Int tag) { MPI_Comm communicator = mpi_data->getMPICommunicator(); CommunicationRequestMPI * request = new CommunicationRequestMPI(prank, receiver); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Isend(buffer, size, type, receiver, tag, communicator, request->getMPIRequest()); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Isend."); return request; } /* -------------------------------------------------------------------------- */ template CommunicationRequest * StaticCommunicatorMPI::asyncReceive(T * buffer, Int size, Int sender, Int tag) { MPI_Comm communicator = mpi_data->getMPICommunicator(); CommunicationRequestMPI * request = new CommunicationRequestMPI(sender, prank); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Irecv(buffer, size, type, sender, tag, communicator, request->getMPIRequest()); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Irecv."); return request; } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::probe(Int sender, Int tag, CommunicationStatus & status) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Status mpi_status; #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Probe(sender, tag, communicator, &mpi_status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Probe."); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); int count; MPI_Get_count(&mpi_status, type, &count); status.setSource(mpi_status.MPI_SOURCE); status.setTag(mpi_status.MPI_TAG); status.setSize(count); } /* -------------------------------------------------------------------------- */ bool StaticCommunicatorMPI::testRequest(CommunicationRequest * request) { MPI_Status status; int flag; CommunicationRequestMPI * req_mpi = static_cast(request); MPI_Request * req = req_mpi->getMPIRequest(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Test(req, &flag, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Test."); return (flag != 0); } /* -------------------------------------------------------------------------- */ void StaticCommunicatorMPI::wait(CommunicationRequest * request) { MPI_Status status; CommunicationRequestMPI * req_mpi = static_cast(request); MPI_Request * req = req_mpi->getMPIRequest(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Wait(req, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Wait."); } /* -------------------------------------------------------------------------- */ void StaticCommunicatorMPI::waitAll( std::vector & requests) { MPI_Status status; std::vector::iterator it; for (it = requests.begin(); it != requests.end(); ++it) { MPI_Request * req = static_cast(*it)->getMPIRequest(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Wait(req, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Wait."); } } /* -------------------------------------------------------------------------- */ void StaticCommunicatorMPI::barrier() { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Barrier(communicator); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::reduce(T * values, int nb_values, const SynchronizerOperation & op, int root) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Reduce(MPI_IN_PLACE, values, nb_values, type, MPITypeWrapper::getMPISynchronizerOperation(op), root, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Allreduce."); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::allReduce(T * values, int nb_values, const SynchronizerOperation & op) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Allreduce(MPI_IN_PLACE, values, nb_values, type, MPITypeWrapper::getMPISynchronizerOperation(op), communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Allreduce."); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::allGather(T * values, int nb_values) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Allgather(MPI_IN_PLACE, nb_values, type, values, nb_values, type, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Allgather."); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::allGatherV(T * values, int * nb_values) { MPI_Comm communicator = mpi_data->getMPICommunicator(); int * displs = new int[psize]; displs[0] = 0; for (int i = 1; i < psize; ++i) { displs[i] = displs[i - 1] + nb_values[i - 1]; } MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Allgatherv(MPI_IN_PLACE, *nb_values, type, values, nb_values, displs, type, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Gather."); delete[] displs; } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::gather(T * values, int nb_values, int root) { MPI_Comm communicator = mpi_data->getMPICommunicator(); T * send_buf = NULL, * recv_buf = NULL; if (prank == root) { send_buf = (T *)MPI_IN_PLACE; recv_buf = values; } else { send_buf = values; } MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Gather(send_buf, nb_values, type, recv_buf, nb_values, type, root, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Gather."); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::gatherV(T * values, int * nb_values, int root) { MPI_Comm communicator = mpi_data->getMPICommunicator(); int * displs = NULL; if (prank == root) { displs = new int[psize]; displs[0] = 0; for (int i = 1; i < psize; ++i) { displs[i] = displs[i - 1] + nb_values[i - 1]; } } T * send_buf = NULL, * recv_buf = NULL; if (prank == root) { send_buf = (T *)MPI_IN_PLACE; recv_buf = values; } else send_buf = values; MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Gatherv(send_buf, *nb_values, type, recv_buf, nb_values, displs, type, root, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Gather."); if (prank == root) { delete[] displs; } } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::broadcast(T * values, int nb_values, int root) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Bcast(values, nb_values, type, root, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Gather."); } /* -------------------------------------------------------------------------- */ int StaticCommunicatorMPI::getMaxTag() { return MPI_TAG_UB; } /* -------------------------------------------------------------------------- */ int StaticCommunicatorMPI::getMinTag() { return 0; } /* -------------------------------------------------------------------------- */ // template // MPI_Datatype StaticCommunicatorMPI::getMPIDatatype() { // return MPI_DATATYPE_NULL; // } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_CHAR; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_FLOAT; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_DOUBLE; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_LONG_DOUBLE; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_INT; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_UNSIGNED; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_LONG; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_UNSIGNED_LONG; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_LONG_LONG; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_UNSIGNED_LONG_LONG; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype >() { return MPI_DOUBLE_INT; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype >() { return MPI_FLOAT_INT; } /* -------------------------------------------------------------------------- */ /* Template instantiation */ /* -------------------------------------------------------------------------- */ #define AKANTU_MPI_COMM_INSTANTIATE(T) \ template void StaticCommunicatorMPI::send(T * buffer, Int size, \ Int receiver, Int tag); \ template void StaticCommunicatorMPI::receive(T * buffer, Int size, \ Int sender, Int tag); \ template CommunicationRequest * StaticCommunicatorMPI::asyncSend( \ T * buffer, Int size, Int receiver, Int tag); \ template CommunicationRequest * StaticCommunicatorMPI::asyncReceive( \ T * buffer, Int size, Int sender, Int tag); \ template void StaticCommunicatorMPI::probe(Int sender, Int tag, \ CommunicationStatus & status); \ template void StaticCommunicatorMPI::allGather(T * values, \ int nb_values); \ template void StaticCommunicatorMPI::allGatherV(T * values, \ int * nb_values); \ template void StaticCommunicatorMPI::gather(T * values, int nb_values, \ int root); \ template void StaticCommunicatorMPI::gatherV(T * values, int * nb_values, \ int root); \ template void StaticCommunicatorMPI::broadcast(T * values, int nb_values, \ int root); \ template void StaticCommunicatorMPI::allReduce( \ T * values, int nb_values, const SynchronizerOperation & op) AKANTU_MPI_COMM_INSTANTIATE(Real); AKANTU_MPI_COMM_INSTANTIATE(UInt); AKANTU_MPI_COMM_INSTANTIATE(Int); AKANTU_MPI_COMM_INSTANTIATE(char); template void StaticCommunicatorMPI::send >( SCMinMaxLoc * buffer, Int size, Int receiver, Int tag); template void StaticCommunicatorMPI::receive >( SCMinMaxLoc * buffer, Int size, Int sender, Int tag); template CommunicationRequest * StaticCommunicatorMPI::asyncSend >( SCMinMaxLoc * buffer, Int size, Int receiver, Int tag); template CommunicationRequest * StaticCommunicatorMPI::asyncReceive >( SCMinMaxLoc * buffer, Int size, Int sender, Int tag); template void StaticCommunicatorMPI::probe >( Int sender, Int tag, CommunicationStatus & status); template void StaticCommunicatorMPI::allGather >( SCMinMaxLoc * values, int nb_values); template void StaticCommunicatorMPI::allGatherV >( SCMinMaxLoc * values, int * nb_values); template void StaticCommunicatorMPI::gather >( SCMinMaxLoc * values, int nb_values, int root); template void StaticCommunicatorMPI::gatherV >( SCMinMaxLoc * values, int * nb_values, int root); template void StaticCommunicatorMPI::broadcast >( SCMinMaxLoc * values, int nb_values, int root); template void StaticCommunicatorMPI::allReduce >( SCMinMaxLoc * values, int nb_values, const SynchronizerOperation & op); #if AKANTU_INTEGER_SIZE > 4 AKANTU_MPI_COMM_INSTANTIATE(int); #endif __END_AKANTU__