diff --git a/packages/84_petsc.cmake b/packages/84_petsc.cmake index 9c68a4d6f..ad0483993 100644 --- a/packages/84_petsc.cmake +++ b/packages/84_petsc.cmake @@ -1,42 +1,41 @@ #=============================================================================== # @file petsc.cmake # # @author Nicolas Richart # @author Alejandro M. Aragón # @author Aurelia Cuba Ramos # # @date Mon Nov 21 18:19:15 2011 # # @brief package description for PETSc support # # @section LICENSE # # Copyright (©) 2010-2011 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 . # #=============================================================================== -option(AKANTU_USE_PETSC "Use external library PETSc" OFF) -add_optional_external_package(PETSc "Add PETSc support in akantu" OFF DEPENDS MPI ARGS COMPONENT ARGS CXX) +add_optional_external_package(PETSc "Add PETSc support in akantu" OFF ARGS COMPONENT ARGS CXX) +add_internal_package_dependencies(petsc parallel) -mark_as_advanced(AKANTU_USE_PETSC) -set(AKANTU_PETSC_FILES - solver/solver_petsc.hh - solver/solver_petsc.cc +set(AKANTU_PETSC_FILES solver/petsc_matrix.hh solver/petsc_matrix.cc solver/petsc_matrix_inline_impl.cc + solver/solver_petsc.hh + solver/solver_petsc.cc ) diff --git a/src/common/aka_error.cc b/src/common/aka_error.cc index 835979f84..e9d1bf08a 100644 --- a/src/common/aka_error.cc +++ b/src/common/aka_error.cc @@ -1,337 +1,361 @@ /** * @file aka_error.cc * * @author Nicolas Richart * * @date Mon Sep 06 00:18:58 2010 * * @brief handling of errors * * @section LICENSE * * Copyright (©) 2010-2011 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 "aka_config.hh" #include "aka_error.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include #include #include #include #include #if defined(AKANTU_USE_OBSOLETE_GETTIMEOFDAY) # include #else # include #endif #ifdef AKANTU_USE_MPI +#if defined(AKANTU_CORE_CXX11) +#if defined(__INTEL_COMPILER) +//#pragma warning ( disable : 383 ) +#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu +#elif (defined(__GNUC__) || defined(__GNUG__)) +# define GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) +# if GCC_VERSION > 40600 +# pragma GCC diagnostic push +# endif +# pragma GCC diagnostic ignored "-Wliteral-suffix" +#endif +#endif # include +#if defined(AKANTU_CORE_CXX11) +#if defined(__INTEL_COMPILER) +//#pragma warning ( disable : 383 ) +#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu +#elif defined(__GNUG__) +# if GCC_VERSION > 40600 +# pragma GCC diagnostic pop +# else +# pragma GCC diagnostic warning "-Wliteral-suffix" +# endif +#endif +#endif #endif #define __BEGIN_AKANTU_DEBUG__ namespace akantu { namespace debug { #define __END_AKANTU_DEBUG__ } } /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU_DEBUG__ static void printBacktraceAndExit(int sig) { printBacktrace(sig); debugger.exit(-50); } /* ------------------------------------------------------------------------ */ void initSignalHandler() { struct sigaction action; action.sa_handler = &printBacktraceAndExit; sigemptyset(&(action.sa_mask)); action.sa_flags = SA_RESETHAND; sigaction(SIGSEGV, &action, NULL); sigaction(SIGABRT, &action, NULL); } /* ------------------------------------------------------------------------ */ std::string demangle(const char *symbol) { int status; std::string result; char *demangled_name; if ((demangled_name = abi::__cxa_demangle(symbol, NULL, 0, &status)) != NULL) { result = demangled_name; free(demangled_name); } else { result = symbol; } return result; } /* ------------------------------------------------------------------------ */ std::string exec(std::string cmd) { FILE *pipe = popen(cmd.c_str(), "r"); if (!pipe) return ""; char buffer[1024]; std::string result = ""; while (!feof(pipe)) { if (fgets(buffer, 128, pipe) != NULL) result += buffer; } result = result.substr(0, result.size() - 1); pclose(pipe); return result; } /* ------------------------------------------------------------------------ */ void printBacktrace(__attribute__((unused)) int sig) { AKANTU_DEBUG_INFO("Caught signal " << sig << "!"); #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) std::string me = ""; char buf[1024]; /* The manpage says it won't null terminate. Let's zero the buffer. */ memset(buf, 0, sizeof(buf)); /* Note we use sizeof(buf)-1 since we may need an extra char for NUL. */ if (readlink("/proc/self/exe", buf, sizeof(buf) - 1)) me = std::string(buf); std::ifstream inmaps; inmaps.open("/proc/self/maps"); std::map addr_map; std::string line; while (inmaps.good()) { std::getline(inmaps, line); std::stringstream sstr(line); size_t first = line.find('-'); std::stringstream sstra(line.substr(0, first)); size_t addr; sstra >> std::hex >> addr; std::string lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; if (lib != "" && addr_map.find(lib) == addr_map.end()) { addr_map[lib] = addr; } } if (me != "") addr_map[me] = 0; #endif const size_t max_depth = 100; size_t stack_depth; void *stack_addrs[max_depth]; char **stack_strings; size_t i; stack_depth = backtrace(stack_addrs, max_depth); stack_strings = backtrace_symbols(stack_addrs, stack_depth); std::cerr << "BACKTRACE : " << stack_depth << " stack frames." << std::endl; size_t w = size_t(std::floor(log(double(stack_depth)) / std::log(10.)) + 1); /// -1 to remove the call to the printBacktrace function for (i = 1; i < stack_depth; i++) { std::cerr << std::dec << " [" << std::setw(w) << i << "] "; std::string bt_line(stack_strings[i]); size_t first, second; if ((first = bt_line.find('(')) != std::string::npos && (second = bt_line.find('+')) != std::string::npos) { std::string location = bt_line.substr(0, first); #if defined(READLINK_COMMAND) location = exec(std::string(BOOST_PP_STRINGIZE(READLINK_COMMAND)) + std::string(" -f ") + location); #endif std::string call = demangle(bt_line.substr(first + 1, second - first - 1).c_str()); size_t f = bt_line.find('['); size_t s = bt_line.find(']'); std::string address = bt_line.substr(f + 1, s - f - 1); std::stringstream sstra(address); size_t addr; sstra >> std::hex >> addr; std::cerr << location << " [" << call << "]"; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) std::map ::iterator it = addr_map.find(location); if (it != addr_map.end()) { std::stringstream syscom; syscom << BOOST_PP_STRINGIZE(ADDR2LINE_COMMAND) << " 0x" << std::hex << (addr - it->second) << " -i -e " << location; std::string line = exec(syscom.str()); std::cerr << " (" << line << ")" << std::endl; } else { #endif std::cerr << " (0x" << std::hex << addr << ")" << std::endl; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) } #endif } else { std::cerr << bt_line << std::endl; } } free(stack_strings); std::cerr << "END BACKTRACE" << std::endl; } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ Debugger::Debugger() { cout = &std::cerr; level = dblWarning; parallel_context = ""; file_open = false; print_backtrace = false; } /* ------------------------------------------------------------------------ */ Debugger::~Debugger() { if (file_open) { dynamic_cast (cout)->close(); delete cout; } } /* ------------------------------------------------------------------------ */ void Debugger::exit(int status) { if (status != EXIT_SUCCESS && status != -50) { //#ifndef AKANTU_NDEBUG // int * a = NULL; // *a = 1; // #endif if(this->print_backtrace) akantu::debug::printBacktrace(15); } #ifdef AKANTU_USE_MPI if (status != EXIT_SUCCESS) MPI_Abort(MPI_COMM_WORLD, MPI_ERR_UNKNOWN); #endif std::exit(status); // not called when compiled with MPI due to MPI_Abort, but // MPI_Abort does not have the noreturn attribute } /*------------------------------------------------------------------------- */ void Debugger::throwException(const std::string & info, const std::string & file, unsigned int line, __attribute__((unused)) bool silent, __attribute__((unused)) const std::string & location) const throw(akantu::debug::Exception) { #if !defined(AKANTU_NDEBUG) if (!silent) { printMessage("###", dblWarning, info + location); } #endif debug::Exception ex(info, file, line); throw ex; } /* ------------------------------------------------------------------------ */ void Debugger::printMessage(const std::string & prefix, const DebugLevel & level, const std::string & info) const { if (this->level >= level) { #if defined(AKANTU_USE_OBSOLETE_GETTIMEOFDAY) struct timeval time; gettimeofday(&time, NULL); double timestamp = time.tv_sec * 1e6 + time.tv_usec; /*in us*/ #else struct timespec time; clock_gettime(CLOCK_REALTIME_COARSE, &time); double timestamp = time.tv_sec * 1e6 + time.tv_nsec * 1e-3; /*in us*/ #endif *(cout) << parallel_context << "{" << (unsigned int)timestamp << "} " << prefix << " " << info << std::endl; } } /* ------------------------------------------------------------------------ */ void Debugger::setDebugLevel(const DebugLevel & level) { this->level = level; } /* ------------------------------------------------------------------------ */ const DebugLevel & Debugger::getDebugLevel() const { return level; } /* ------------------------------------------------------------------------ */ void Debugger::setLogFile(const std::string & filename) { if (file_open) { dynamic_cast (cout)->close(); delete cout; } std::ofstream *fileout = new std::ofstream(filename.c_str()); file_open = true; cout = fileout; } std::ostream & Debugger::getOutputStream() { return *cout; } /* ------------------------------------------------------------------------ */ void Debugger::setParallelContext(int rank, int size) { std::stringstream sstr; UInt pad = std::ceil(std::log10(size)); sstr << "<" << getpid() << ">[R" << std::setfill(' ') << std::right << std::setw(pad) << rank << "|S" << size << "] "; parallel_context = sstr.str(); } void setDebugLevel(const DebugLevel & level) { debugger.setDebugLevel(level); } const DebugLevel & getDebugLevel() { return debugger.getDebugLevel(); } /* -------------------------------------------------------------------------- */ void exit(int status) { debugger.exit(status); } __END_AKANTU_DEBUG__ diff --git a/src/solver/petsc_matrix.cc b/src/solver/petsc_matrix.cc index be982b6bd..5ad256638 100644 --- a/src/solver/petsc_matrix.cc +++ b/src/solver/petsc_matrix.cc @@ -1,717 +1,717 @@ /** * @file petsc_matrix.cc * @author Aurelia Cuba Ramos * @date Mon Jul 21 17:40:41 2014 * * @brief Implementation of PETSc matrix class * * @section LICENSE * * Copyright (©) 2010-2011 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 "petsc_matrix.hh" /* -------------------------------------------------------------------------- */ - +#include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ PetscMatrix::PetscMatrix(UInt size, const SparseMatrixType & sparse_matrix_type, UInt nb_degree_of_freedom, const ID & id, const MemoryID & memory_id) : PetscMatrix(size, sparse_matrix_type, nb_degree_of_freedom, id, memory_id) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); nb_proc = comm.getNbProc(); UInt nb_local_nodes = model.getMesh().getNbNodes(); - MatCreate(&comm, &mat); + MatCreate(&comm, mat); MatSetSizes(mat, nb_degree_of_freedom*nb_local_nodes, nb_degree_of_freedom*nb_local_nodes, size, size); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PetscMatrix::PetscMatrix(const PetscMatrix & matrix, const ID & id, const MemoryID & memory_id) : PetscMatrix(matrix, id, memory_id) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); nb_proc = comm.getNbProc(); UInt nb_local_nodes = model.getMesh().getNbNodes(); MatCreate(&comm, &mat); MatSetSizes(mat, nb_degree_of_freedom*nb_local_nodes, nb_degree_of_freedom*nb_local_nodes, size, size); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PetscMatrix::~PetscMatrix() { AKANTU_DEBUG_IN(); MatDestroy(&mat); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PetscMatrix::createAkantuToPetscMap(Mesh & mesh) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int nb_proc = comm.getNbProc(); Int rank = comm.whoAmI(); this->dof_synchronizer = &const_cast(dof_synchronizer); //initialize vector to store the number of local and master nodes for each processor Vector master_local_ndofs_per_proc(nb_proc); //find the number of dofs corresponding to master or local nodes and store their akantu global dof number Real nb_dofs = dof_sychronizer.getNbDOFs(); Real nb_local_master_dofs = 0; Array local_master_dofs(nb_dofs); Array::scalar_iterator it_dof = local_master_dofs.begin(); Array dof_types = dof_sychronizer.getDOFTypes(); Array::scalar_iterator it_dof_type = dof_types.begin(); for (UInt i = 0; i , UInt>; clearProfile(); this->dof_synchronizer = &const_cast(dof_synchronizer); coordinate_list_map::iterator irn_jcn_k_it; Int * eq_nb_val = dof_synchronizer.getGlobalDOFEquationNumbers().storage(); Mesh::type_iterator it = mesh.firstType(mesh.getSpatialDimension(), _not_ghost, _ek_not_defined); Mesh::type_iterator end = mesh.lastType (mesh.getSpatialDimension(), _not_ghost, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; UInt * conn_val = mesh.getConnectivity(*it, _not_ghost).storage(); Int * local_eq_nb_val = new Int[nb_degree_of_freedom * nb_nodes_per_element]; for (UInt e = 0; e < nb_element; ++e) { Int * tmp_local_eq_nb_val = local_eq_nb_val; for (UInt i = 0; i < nb_nodes_per_element; ++i) { UInt n = conn_val[i]; for (UInt d = 0; d < nb_degree_of_freedom; ++d) { *tmp_local_eq_nb_val++ = eq_nb_val[n * nb_degree_of_freedom + d]; } // memcpy(tmp_local_eq_nb_val, eq_nb_val + n * nb_degree_of_freedom, nb_degree_of_freedom * sizeof(Int)); // tmp_local_eq_nb_val += nb_degree_of_freedom; } for (UInt i = 0; i < size_mat; ++i) { UInt c_irn = local_eq_nb_val[i]; if(c_irn < size) { UInt j_start = (sparse_matrix_type == _symmetric) ? i : 0; for (UInt j = j_start; j < size_mat; ++j) { UInt c_jcn = local_eq_nb_val[j]; if(c_jcn < size) { MatSetValue(this->mat,1,c_irn, 1, c_jcn, INSERT_VALUES); if (irn_jcn_k_it == irn_jcn_k.end()) { irn_jcn_k[irn_jcn] = nb_non_zero; irn.push_back(c_irn + 1); jcn.push_back(c_jcn + 1); nb_non_zero++; } } } } } conn_val += nb_nodes_per_element; } delete [] local_eq_nb_val; } /// for pbc @todo correct it for parallel if(StaticCommunicator::getStaticCommunicator().getNbProc() == 1) { for (UInt i = 0; i < size; ++i) { KeyCOO irn_jcn = key(i, i); irn_jcn_k_it = irn_jcn_k.find(irn_jcn); if(irn_jcn_k_it == irn_jcn_k.end()) { irn_jcn_k[irn_jcn] = nb_non_zero; irn.push_back(i + 1); jcn.push_back(i + 1); nb_non_zero++; } } } a.resize(nb_non_zero); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PetscMatrix::applyBoundaryNormal(Array & boundary_normal, Array & EulerAngles, Array & rhs, const Array & matrix, Array & rhs_rotated) { AKANTU_DEBUG_IN(); const UInt dim = nb_degree_of_freedom; const UInt nb_nodes = boundary_normal.getSize(); Matrix Ti(dim, dim); Matrix Tj(dim, dim); Matrix small_matrix(dim, dim); Matrix Ti_small_matrix(dim, dim); Matrix Ti_small_matrix_Tj(dim, dim); Matrix small_rhs(dim, 1); Matrix Ti_small_rhs(dim, 1); //Transformation matrix based on euler angles X_1Y_2Z_3 const DOFSynchronizer::GlobalEquationNumberMap & local_eq_num_to_global = dof_synchronizer->getGlobalEquationNumberToLocal(); Int * irn_val = irn.storage(); Int * jcn_val = jcn.storage(); Int * eq_nb_val = dof_synchronizer->getGlobalDOFEquationNumbers().storage(); Int * eq_nb_rhs_val = dof_synchronizer->getLocalDOFEquationNumbers().storage(); Array * nodes_rotated = new Array (nb_nodes, dim, "nodes_rotated"); nodes_rotated->clear(); Array * local_eq_nb_val_node_i = new Array (1, nb_degree_of_freedom, "local_eq_nb_val_node_i"); local_eq_nb_val_node_i->clear(); Array * local_eq_nb_val_node_j = new Array (1, nb_degree_of_freedom, "local_eq_nb_val_node_j"); local_eq_nb_val_node_j->clear(); for (UInt i = 0; i < nb_non_zero; ++i) { UInt ni = local_eq_num_to_global.find(*irn_val - 1)->second; UInt node_i = (ni - ni % nb_degree_of_freedom) / nb_degree_of_freedom; UInt nj = local_eq_num_to_global.find(*jcn_val - 1)->second; UInt node_j = (nj - nj % nb_degree_of_freedom) / nb_degree_of_freedom; bool constrain_ij = false; for (UInt j = 0; j < dim; j++){ if ( (boundary_normal(node_i, j) || boundary_normal(node_j, j)) ) { constrain_ij = true; break; } } if(constrain_ij){ if(dim==2){ Real Theta=EulerAngles(node_i,0); Ti(0,0)=cos(Theta); Ti(0,1)=-sin(Theta); Ti(1,1)=cos(Theta); Ti(1,0)=sin(Theta); Theta=EulerAngles(node_j,0); Tj(0,0)=cos(Theta); Tj(0,1)=-sin(Theta); Tj(1,1)=cos(Theta); Tj(1,0)=sin(Theta); } else if(dim==3){ Real Theta_x=EulerAngles(node_i,0); Real Theta_y=EulerAngles(node_i,1); Real Theta_z=EulerAngles(node_i,2); Ti(0,0)=cos(Theta_y)*cos(Theta_z); Ti(0,1)=-cos(Theta_y)*sin(Theta_z); Ti(0,2)=sin(Theta_y); Ti(1,0)=cos(Theta_x)*sin(Theta_z)+cos(Theta_z)*sin(Theta_x)*sin(Theta_y); Ti(1,1)=cos(Theta_x)*cos(Theta_z)-sin(Theta_x)*sin(Theta_y)*sin(Theta_z); Ti(1,2)=-cos(Theta_y)*sin(Theta_x); Ti(2,0)=sin(Theta_x)*sin(Theta_z)-cos(Theta_x)*cos(Theta_z)*sin(Theta_y); Ti(2,1)=cos(Theta_z)*sin(Theta_x)+cos(Theta_x)*sin(Theta_y)*sin(Theta_z); Ti(2,2)=cos(Theta_x)*cos(Theta_y); Theta_x=EulerAngles(node_j,0); Theta_y=EulerAngles(node_j,1); Theta_z=EulerAngles(node_j,2); Tj(0,0)=cos(Theta_y)*cos(Theta_z); Tj(0,1)=-cos(Theta_y)*sin(Theta_z); Tj(0,2)=sin(Theta_y); Tj(1,0)=cos(Theta_x)*sin(Theta_z)+cos(Theta_z)*sin(Theta_x)*sin(Theta_y); Tj(1,1)=cos(Theta_x)*cos(Theta_z)-sin(Theta_x)*sin(Theta_y)*sin(Theta_z); Tj(1,2)=-cos(Theta_y)*sin(Theta_x); Tj(2,0)=sin(Theta_x)*sin(Theta_z)-cos(Theta_x)*cos(Theta_z)*sin(Theta_y); Tj(2,1)=cos(Theta_z)*sin(Theta_x)+cos(Theta_x)*sin(Theta_y)*sin(Theta_z); Tj(2,2)=cos(Theta_x)*cos(Theta_y); } for (UInt j = 0; j < nb_degree_of_freedom; ++j){ local_eq_nb_val_node_i->storage()[j] = eq_nb_val[node_i * nb_degree_of_freedom + j]; local_eq_nb_val_node_j->storage()[j] = eq_nb_val[node_j * nb_degree_of_freedom + j]; } for (UInt j = 0; j < nb_degree_of_freedom; ++j) { for (UInt k = 0; k < nb_degree_of_freedom; ++k) { KeyCOO jcn_irn = key(local_eq_nb_val_node_i->storage()[j], local_eq_nb_val_node_j->storage()[k]); coordinate_list_map::iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn); small_matrix(j, k) = matrix.storage()[irn_jcn_k_it->second]; } small_rhs(j,0)=rhs.storage()[eq_nb_rhs_val[node_i*dim+j]]; } Ti_small_matrix.mul < false, false > (Ti, small_matrix); Ti_small_matrix_Tj.mul < false, true> (Ti_small_matrix, Tj); Ti_small_rhs.mul(Ti, small_rhs); /*for (UInt j = 0; j < nb_degree_of_freedom; ++j) { for (UInt k = 0; k < nb_degree_of_freedom; ++k) { KeyCOO jcn_irn = key(local_eq_nb_val_node_i->storage()[j], local_eq_nb_val_node_j->storage()[k]); coordinate_list_map::iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn); a.storage()[irn_jcn_k_it->second] = T_small_matrix_T(j,k); } if(node_constrain==node_i && !( nodes_rotated->storage()[node_i])) rhs.storage()[eq_nb_rhs_val[node_i*dim+j]]=T_small_rhs(j,0); }*/ KeyCOO jcn_irn = key(ni, nj); coordinate_list_map::iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn); a.storage()[irn_jcn_k_it->second] = Ti_small_matrix_Tj(ni % nb_degree_of_freedom, nj % nb_degree_of_freedom); //this->saveMatrix("stiffness_partial.out"); if (!(nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]])) { rhs_rotated.storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = Ti_small_rhs(ni % nb_degree_of_freedom, 0); nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = true; } } else{ if (!(nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]])) { rhs_rotated.storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = rhs.storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]]; nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = true; } } irn_val++; jcn_val++; } delete local_eq_nb_val_node_i; delete local_eq_nb_val_node_j; delete nodes_rotated; this->saveMatrix("stiffness_rotated.out"); applyBoundary(boundary_normal); /*Real norm = 0; UInt n_zeros = 0; for (UInt j = 0; j < nb_degree_of_freedom; ++j) { norm += Normal.storage()[j] * Normal.storage()[j]; if (std::abs(Normal.storage()[j]) < Math::getTolerance()) n_zeros++; } norm = sqrt(norm); AKANTU_DEBUG_ASSERT(norm > Math::getTolerance(), "The normal is not right"); for (UInt j = 0; j < nb_degree_of_freedom; ++j) Normal.storage()[j] /= norm; if (n_zeros == nb_degree_of_freedom - 1) { UInt nb_nodes = boundary_normal.getSize(); for (UInt i = 0; i < nb_nodes; ++i) { if (boundary_normal(i, 0)) for (UInt j = 0; j < nb_degree_of_freedom; ++j) boundary(i, j) += Normal(0, j); } } else { const DOFSynchronizer::GlobalEquationNumberMap & local_eq_num_to_global = dof_synchronizer->getGlobalEquationNumberToLocal(); Int * irn_val = irn.storage(); Int * eq_nb_val = dof_synchronizer->getGlobalDOFEquationNumbers().storage(); Array * local_eq_nb_val = new Array (1, nb_degree_of_freedom, "local_eq_nb_val"); local_eq_nb_val->clear(); UInt nb_nodes = boundary_normal.getSize(); Array * node_set = new Array (nb_nodes, 1, "node_set"); node_set->clear(); for (UInt i = 0; i < nb_non_zero; ++i) { UInt ni = local_eq_num_to_global.find(*irn_val - 1)->second; UInt node_i = (ni - ni % nb_degree_of_freedom) / nb_degree_of_freedom; if (boundary_normal.storage()[ni]) { if ((!number.storage()[ni]) && (!node_set->storage()[node_i])) { UInt normal_component_i = ni % nb_degree_of_freedom; //DOF of node node_i to be removed if (std::abs(Normal.storage()[normal_component_i]) > Math::getTolerance()) { for (UInt j = 0; j < nb_degree_of_freedom; ++j) local_eq_nb_val->storage()[j] = eq_nb_val[node_i * nb_degree_of_freedom + j]; UInt DOF_remove = local_eq_nb_val->storage()[normal_component_i]; KeyCOO jcn_irn = key(DOF_remove, DOF_remove); coordinate_list_map::iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn); Real aii = a.storage()[irn_jcn_k_it->second]; Real normal_constant = (1 - aii) / (Normal.storage()[normal_component_i] * Normal.storage()[normal_component_i]); for (UInt j = 0; j < nb_degree_of_freedom; ++j) { UInt line_j = local_eq_nb_val->storage()[j]; for (UInt k = 0; k < nb_degree_of_freedom; ++k) { UInt column_k = local_eq_nb_val->storage()[k]; jcn_irn = key(line_j, column_k); if ((line_j == column_k) && (line_j == DOF_remove)) a.storage()[irn_jcn_k_it->second] = 1; else a.storage()[irn_jcn_k_it->second] += Normal.storage()[j] * Normal.storage()[k] * normal_constant; } } number.storage()[ni]++; node_set->storage()[node_i] = false; } } } irn_val++; } delete local_eq_nb_val; delete node_set; }*/ AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PetscMatrix::applyBoundary(const Array & boundary, Real block_val) { AKANTU_DEBUG_IN(); const DOFSynchronizer::GlobalEquationNumberMap & local_eq_num_to_global = dof_synchronizer->getGlobalEquationNumberToLocal(); Int * irn_val = irn.storage(); Int * jcn_val = jcn.storage(); Real * a_val = a.storage(); for (UInt i = 0; i < nb_non_zero; ++i) { UInt ni = local_eq_num_to_global.find(*irn_val - 1)->second; UInt nj = local_eq_num_to_global.find(*jcn_val - 1)->second; if(boundary.storage()[ni] || boundary.storage()[nj]) { if (*irn_val != *jcn_val) { *a_val = 0; } else { if(dof_synchronizer->getDOFTypes()(ni) >= 0) { *a_val = 0; } else { *a_val = block_val; } } } irn_val++; jcn_val++; a_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PetscMatrix::removeBoundary(const Array & boundary) { AKANTU_DEBUG_IN(); if(irn_save) delete irn_save; if(jcn_save) delete jcn_save; irn_save = new Array(irn, true); jcn_save = new Array(jcn, true); UInt n = boundary.getSize()*boundary.getNbComponent(); UInt * perm = new UInt[n]; size_save = size; size = 0; for (UInt i = 0; i < n; ++i) { if(!boundary.storage()[i]) { perm[i] = size; // std::cout << "perm["<< i <<"] = " << size << std::endl; size++; } } for (UInt i = 0; i < nb_non_zero;) { if(boundary.storage()[irn(i) - 1] || boundary.storage()[jcn(i) - 1]) { irn.erase(i); jcn.erase(i); a.erase(i); nb_non_zero--; } else { irn(i) = perm[irn(i) - 1] + 1; jcn(i) = perm[jcn(i) - 1] + 1; i++; } } delete [] perm; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PetscMatrix::restoreProfile() { AKANTU_DEBUG_IN(); irn.resize(irn_save->getSize()); jcn.resize(jcn_save->getSize()); nb_non_zero = irn.getSize(); a.resize(nb_non_zero); size = size_save; memcpy(irn.storage(), irn_save->storage(), irn.getSize()*sizeof(Int)); memcpy(jcn.storage(), jcn_save->storage(), jcn.getSize()*sizeof(Int)); delete irn_save; irn_save = NULL; delete jcn_save; jcn_save = NULL; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PetscMatrix::saveProfile(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.open(filename.c_str()); outfile << "%%MatrixMarket matrix coordinate pattern"; if(sparse_matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; UInt m = size; outfile << m << " " << m << " " << nb_non_zero << std::endl; for (UInt i = 0; i < nb_non_zero; ++i) { outfile << irn.storage()[i] << " " << jcn.storage()[i] << " 1" << std::endl; } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PetscMatrix::saveMatrix(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.precision(std::numeric_limits::digits10); outfile.open(filename.c_str()); outfile << "%%MatrixMarket matrix coordinate real"; if(sparse_matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; outfile << size << " " << size << " " << nb_non_zero << std::endl; for (UInt i = 0; i < nb_non_zero; ++i) { outfile << irn(i) << " " << jcn(i) << " " << a(i) << std::endl; } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Array & operator*=(Array & vect, const PetscMatrix & mat) { AKANTU_DEBUG_IN(); Array result MatMult(this->mat, vect, ) // AKANTU_DEBUG_ASSERT((vect.getSize()*vect.getNbComponent() == mat.getSize()) && // (vect.getNbComponent() == mat.getNbDegreOfFreedom()), // "The size of the matrix and the vector do not match"); const PetscMatrixType & sparse_matrix_type = mat.getPetscMatrixType(); DOFSynchronizer * dof_synchronizer = mat.getDOFSynchronizerPointer(); UInt nb_non_zero = mat.getNbNonZero(); Real * tmp = new Real [vect.getNbComponent() * vect.getSize()]; std::fill_n(tmp, vect.getNbComponent() * vect.getSize(), 0); Int * i_val = mat.getIRN().storage(); Int * j_val = mat.getJCN().storage(); Real * a_val = mat.getA().storage(); Real * vect_val = vect.storage(); for (UInt k = 0; k < nb_non_zero; ++k) { UInt i = *(i_val++); UInt j = *(j_val++); Real a = *(a_val++); UInt local_i = i - 1; UInt local_j = j - 1; if(dof_synchronizer) { local_i = dof_synchronizer->getDOFLocalID(local_i); local_j = dof_synchronizer->getDOFLocalID(local_j); } tmp[local_i] += a * vect_val[local_j]; if((sparse_matrix_type == _symmetric) && (local_i != local_j)) tmp[local_j] += a * vect_val[local_i]; } memcpy(vect_val, tmp, vect.getNbComponent() * vect.getSize() * sizeof(Real)); delete [] tmp; if(dof_synchronizer) dof_synchronizer->reduceSynchronize(vect); AKANTU_DEBUG_OUT(); return vect; } /* -------------------------------------------------------------------------- */ void PetscMatrix::copyContent(const PetscMatrix & matrix) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(nb_non_zero == matrix.getNbNonZero(), "The two matrices don't have the same profiles"); MatCopy(this->mat, matrix->mat, SAME_NONZERO_PATTERN); AKANTU_DEBUG_OUT(); } ///* -------------------------------------------------------------------------- */ //void PetscMatrix::copyProfile(const PetscMatrix & matrix) { // AKANTU_DEBUG_IN(); // irn = matrix.irn; // jcn = matrix.jcn; // nb_non_zero = matrix.nb_non_zero; // irn_jcn_k = matrix.irn_jcn_k; // a.resize(nb_non_zero); // AKANTU_DEBUG_OUT(); //} /* -------------------------------------------------------------------------- */ void PetscMatrix::add(const SparseMatrix & matrix, Real alpha) { AKANTU_DEBUG_ASSERT(nb_non_zero == matrix.getNbNonZero(), "The two matrices don't have the same profile"); // get values and local indices of matrix that has to be added UInt i_global; UInt j_global; Array::scalar_iterator i_global_it = i_global.begin(); Array::scalar_iterator i_global_end = i_global.end(); Array::scalar_iterator j_global_it = j_global.begin(); Int * i = matrix.irn.storage(); Int *j = matrix.jcn.storage(); Real * a_val = matrix.a.storage(); for (; i != i_global_end; ++i_global_it; ++j_global_it; ++i; ++j; ++a_val) { // get the akantu global dof index i_global = dof_synchronizer->getDOFGlobalID(*i); j_global = dof_synchronizer->getDOFGlobalID(*j); // get the petsc gloobal dof index AOApplicationToPetsc(&ao,2,{i_global,j_global}) MatSetValue(*mat,i,j,*a_val,INSERT_VALUES); } } /* -------------------------------------------------------------------------- */ void PetscMatrix::lump(Array & lumped) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT((lumped.getNbComponent() == nb_degree_of_freedom), "The size of the matrix and the vector do not match"); UInt vect_size = size / nb_degree_of_freedom; if(dof_synchronizer) vect_size = dof_synchronizer->getNbDOFs() / nb_degree_of_freedom; lumped.resize(vect_size); lumped.clear(); Int * i_val = irn.storage(); Int * j_val = jcn.storage(); Real * a_val = a.storage(); Real * vect_val = lumped.storage(); for (UInt k = 0; k < nb_non_zero; ++k) { UInt i = *(i_val++); UInt j = *(j_val++); Real a = *(a_val++); UInt local_i = i - 1; UInt local_j = j - 1; if(dof_synchronizer) { local_i = dof_synchronizer->getDOFLocalID(local_i); local_j = dof_synchronizer->getDOFLocalID(local_j); } vect_val[local_i] += a; if(sparse_matrix_type == _symmetric && (i != j)) vect_val[local_j] += a; } if(dof_synchronizer) dof_synchronizer->reduceSynchronize(lumped); AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/solver/petsc_matrix.hh b/src/solver/petsc_matrix.hh index 9779aa442..75d7bc308 100644 --- a/src/solver/petsc_matrix.hh +++ b/src/solver/petsc_matrix.hh @@ -1,121 +1,126 @@ /** * @file petsc_matrix.hh * @author Aurelia Cuba Ramos * @date Mon Jul 21 14:49:49 2014 * * @brief Interface for PETSc matrices * * @section LICENSE * * Copyright (©) 2010-2011 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_PETSC_MATRIX_HH__ #define __AKANTU_PETSC_MATRIX_HH__ /* -------------------------------------------------------------------------- */ #include "sparse_matrix.hh" +struct Mat; +struct AO; + /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ + + class PetscMatrix : SparseMatrix { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ class PetscMatrixValue { public: PetscMatrixValue(){} inline PetscMatrixValue operator = (Real v){ MatSetValue(*mat,i,j,v,INSERT_VALUES); return *this; }; inline PetscMatrixValue operator += (Real v){ // std::cout << "set " << std::abs(Real(i)-Real(j)) << std::endl; MatSetValue(*mat,i,j,v,ADD_VALUES); return *this; }; inline PetscMatrixValue operator -= (Real v){ MatSetValue(*mat,i,j,-1.*v,ADD_VALUES); return *this; }; friend PetscMatrix; private: UInt i,j; Mat * mat; }; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: PetscMatrix(UInt size, const SparseMatrixType & sparse_matrix_type, UInt nb_degree_of_freedom, const ID & id = "petsc_matrix", const MemoryID & memory_id = 0); PetscMatrix(const PetscMatrix & matrix, const ID & id = "petsc_matrix", const MemoryID & memory_id = 0); virtual ~PetscMatrix(); typedef std::pair KeyCOO; typedef unordered_map::type coordinate_list_map; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: - Mat mat; - AO ao; + Mat * mat; + AO * ao; }; __END_AKANTU__ #endif /* __AKANTU_PETSCI_MATRIX_HH__ */ diff --git a/src/solver/solver_petsc.hh b/src/solver/solver_petsc.hh index 00c62a381..f4b7e1c42 100644 --- a/src/solver/solver_petsc.hh +++ b/src/solver/solver_petsc.hh @@ -1,273 +1,273 @@ /** * @file solver_petsc.hh * # @author Alejandro M. Aragón * @author Nicolas Richart * * @date Mon Dec 13 10:48:06 2010 * * @brief Solver class implementation for the petsc solver * * @section LICENSE * * Copyright (©) 2010-2011 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_SOLVER_PETSC_HH__ #define __AKANTU_SOLVER_PETSC_HH__ #include #include "solver.hh" #include "static_communicator.hh" #include "sparse_matrix.hh" +struct Mat; __BEGIN_AKANTU__ - struct SolverPETSc : public Solver, public CommunicatorEventHandler { typedef double value_type; typedef sparse_vector sparse_vector_type; typedef SparseMatrix sparse_matrix_type; Mat A_; //!< linear system matrix Vec x_; //!< Solution vector KSP ksp_; //!< linear solver context bool allocated_; SolverPETSc(int argc, char *argv[]) : allocated_(false) { PetscInitialize(&argc, &argv,NULL,NULL); PetscErrorCode ierr; // create linear solver context ierr = KSPCreate(PETSC_COMM_WORLD, &ksp_);CHKERRCONTINUE(ierr); // initial nonzero guess ierr = KSPSetInitialGuessNonzero(ksp_,PETSC_TRUE); CHKERRCONTINUE(ierr); // set runtime options ierr = KSPSetFromOptions(ksp_);CHKERRCONTINUE(ierr); /* Set linear solver defaults for this problem (optional). - By extracting the KSP and PC contexts from the KSP context, we can then directly call any KSP and PC routines to set various options. - The following four statements are optional; all of these parameters could alternatively be specified at runtime via KSPSetFromOptions(); */ // ierr = KSPGetPC(ksp_,&pc);CHKERRCONTINUE(ierr); // ierr = PCSetType(pc,PCILU);CHKERRCONTINUE(ierr); // ierr = PCSetType(pc,PCJACOBI);CHKERRCONTINUE(ierr); ierr = KSPSetTolerances(ksp_,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRCONTINUE(ierr); } //! Overload operator() to solve system of linear equations sparse_vector_type operator()(const sparse_matrix_type& AA, const sparse_vector_type& bb); //! Overload operator() to obtain reaction vector sparse_vector_type operator()(const sparse_matrix_type& Kpf, const sparse_matrix_type& Kpp, const sparse_vector_type& Up); //! Overload operator() to obtain the addition two vectors sparse_vector_type operator()(const sparse_vector_type& aa, const sparse_vector_type& bb); value_type norm(const sparse_matrix_type& aa, Element_insertion_type it = Add_t); value_type norm(const sparse_vector_type& aa, Element_insertion_type it = Add_t); // NOTE: the destructor will return an error if it is called after MPI_Finalize is // called because it uses collect communication to free-up allocated memory. ~SolverPETSc() { static bool exit = false; if (!exit) { // add finalize PETSc function at exit atexit(finalize); exit = true; } if (allocated_) { PetscErrorCode ierr = MatDestroy(&A_);CHKERRCONTINUE(ierr); ierr = VecDestroy(&x_);CHKERRCONTINUE(ierr); ierr = KSPDestroy(&ksp_);CHKERRCONTINUE(ierr); } } /* from the PETSc library, these are the options that can be passed to the command line Options Database Keys -options_table - Calls PetscOptionsView() -options_left - Prints unused options that remain in the database -objects_left - Prints list of all objects that have not been freed -mpidump - Calls PetscMPIDump() -malloc_dump - Calls PetscMallocDump() -malloc_info - Prints total memory usage -malloc_log - Prints summary of memory usage Options Database Keys for Profiling -log_summary [filename] - Prints summary of flop and timing information to screen. If the filename is specified the summary is written to the file. See PetscLogView(). -log_summary_python [filename] - Prints data on of flop and timing usage to a file or screen. -log_all [filename] - Logs extensive profiling information See PetscLogDump(). -log [filename] - Logs basic profiline information See PetscLogDump(). -log_sync - Log the synchronization in scatters, inner products and norms -log_mpe [filename] - Creates a logfile viewable by the utility Upshot/Nupshot (in MPICH distribution) */ static void finalize() { static bool finalized = false; if (!finalized) { cout<<"*** INFO *** PETSc is finalizing..."< & solution); void solve(); void solveSlave(); virtual void setRHS(Array & rhs); /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const; virtual void onCommunicatorFinalize(const StaticCommunicator & communicator); private: void destroyMumpsData(); inline Int & icntl(UInt i) { return mumps_data.icntl[i - 1]; } inline Int & info(UInt i) { return mumps_data.info[i - 1]; } void initMumpsData(SolverMumpsOptions::ParallelMethod parallel_method); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// mumps data DMUMPS_STRUC_C mumps_data; /// specify if the mumps_data are initialized or not bool is_mumps_data_initialized; UInt prank; /* ------------------------------------------------------------------------ */ /* Local types */ /* ------------------------------------------------------------------------ */ private: SolverMumpsOptions::ParallelMethod parallel_method; bool rhs_is_local; enum SolverMumpsJob { _smj_initialize = -1, _smj_analyze = 1, _smj_factorize = 2, _smj_solve = 3, _smj_analyze_factorize = 4, _smj_factorize_solve = 5, _smj_complete = 6, // analyze, factorize, solve _smj_destroy = -2 }; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "solver_mumps_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const SolverMumps & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_SOLVER_PETSC_HH__ */