Skip to content

Commit

Permalink
better organization. egg carton partially working
Browse files Browse the repository at this point in the history
  • Loading branch information
Brig Bagley committed Feb 24, 2016
1 parent df6aae1 commit fa8ff04
Show file tree
Hide file tree
Showing 20 changed files with 19,970 additions and 461 deletions.
7 changes: 7 additions & 0 deletions src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,13 @@ include_directories("${CUSP_INSTALL_DIR}/src")
# Core / Examples
add_subdirectory(core)

include_directories(core/include)
include_directories(${CMAKE_CURRENT_DIRECTORY})

CUDA_ADD_LIBRARY(FEM_SOLVER FEMSolver.cu FEMSolver.h)
TARGET_LINK_LIBRARIES(FEM_SOLVER FEM_CORE)
ADD_DEPENDENCIES(FEM_SOLVER FEM_CORE)

if (${BUILD_EXAMPLES})
add_subdirectory(examples)
endif()
Expand Down
49 changes: 29 additions & 20 deletions src/core/cuda/FEMSolver.cu → src/FEMSolver.cu
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,9 @@
#include <FEMSolver.h>
#include "FEMSolver.h"
#include "FEM/FEM2D.h"
#include "FEM/FEM3D.h"
#include "cutil.h"
#include "amg.h"
#include <time.h>

FEMSolver::FEMSolver(
std::string fname, bool verbose) :
Expand All @@ -15,10 +20,10 @@ FEMSolver::FEMSolver(
metisSize_(90102), // the Max size of coarsest level
partitionMaxSize_(512), // the largest partition size (use getMaxThreads() to determine for your device)
aggregatorType_(1), // the Aggregator METIS (0) or MIS (1)
convergeType_(ABSOLUTE_CONVERGENCE),// Convergence tolerance algo [ABSOLUTE_CONVERGENCE, RELATIVE_CONVERGENCE]
convergeType_(0), // Convergence tolerance algo [ABSOLUTE_CONVERGENCE (0), RELATIVE_CONVERGENCE (1)]
tolerance_(1e-6), // the convergence tolerance
cycleType_(V_CYCLE), //set the cycle algorithm
solverType_(AMG_SOLVER), // the solving algorithm [AMG_SOLVER,PCG_SOLVER]
cycleType_(0), // set the cycle algorithm
solverType_(0), // the solving algorithm [AMG_SOLVER (0),PCG_SOLVER (1)]
smootherWeight_(1.0), // the weight parameter used in a smoother
proOmega_(0.67), // the weight parameter used in a prolongator smoother
device_(0), // the device number to run on
Expand Down Expand Up @@ -46,23 +51,28 @@ void FEMSolver::solveFEM(Matrix_ell_h* A_h,
//assembly step
if (this->triMesh_ != NULL) {
// 2D fem solving object
FEM2D* fem2d = new FEM2D;
FEM2D fem2d;
Vector_d_CG RHS(this->triMesh_->vertices.size(), 0.0);
fem2d->initializeWithTriMesh(this->triMesh_);
fem2d->assemble(this->triMesh_, A_device, RHS);
delete fem2d;
fem2d.initializeWithTriMesh(this->triMesh_);
fem2d.assemble(this->triMesh_, A_device, RHS);
} else {
// 3D fem solving object
FEM3D* fem3d = new FEM3D;
FEM3D fem3d;
Vector_d_CG RHS(this->tetMesh_->vertices.size(), 0.0);
fem3d->initializeWithTetMesh(this->tetMesh_);
fem3d->assemble(this->tetMesh_, A_device, RHS, true);
delete fem3d;
fem3d.initializeWithTetMesh(this->tetMesh_);
fem3d.assemble(this->tetMesh_, A_device, RHS, true);
}
//copy back to the host
cudaThreadSynchronize();
//register configuration parameters
AMG<Matrix_h, Vector_h> amg(this);
AMG<Matrix_h, Vector_h> amg(this->verbose_, this->convergeType_,
this->cycleType_, this->solverType_, this->tolerance_,
this->cycleIters_, this->maxIters_, this->maxLevels_,
this->topSize_, this->smootherWeight_, this->preInnerIters_,
this->postInnerIters_, this->postRelaxes_, this->dsType_,
this->metisSize_, this->partitionMaxSize_, this->proOmega_,
this->aggregatorType_,
this->triMesh_, this->tetMesh_);
//setup device
Matrix_ell_d A_d(A_device);
amg.setup(A_d);
Expand Down Expand Up @@ -127,7 +137,6 @@ void FEMSolver::checkMatrixForValidContents(Matrix_ell_h* A_h) {
}

void FEMSolver::getMatrixFromMesh(Matrix_ell_h* A_h) {
srand48(0);
if (this->triMesh_ == NULL && this->tetMesh_ == NULL)
exit(0);
Matrix_ell_d_CG Aell_d;
Expand Down Expand Up @@ -323,7 +332,7 @@ int FEMSolver::readMatlabSparseMatrix(const std::string &filename, Matrix_ell_h
return 0;
}

int FEMSolver::readMatlabNormalMatrix(const std::string &filename, vector<double> *A_h) {
int FEMSolver::readMatlabNormalMatrix(const std::string &filename, std::vector<double> *A_h) {
//read in the description header
std::ifstream in(filename.c_str());
if (!in.is_open()) {
Expand Down Expand Up @@ -553,12 +562,12 @@ void FEMSolver::writeVTK(std::vector <float> values, std::string fname)
}
}

void FEMSolver::printElementWithHeader(vector<double>& test, unsigned int index)
void FEMSolver::printElementWithHeader(std::vector<double>& test, unsigned int index)
{
std::cout << "element #" << index << " = " << test[index] << std::endl;
}

void FEMSolver::printMatlabReadContents(vector<double>& test)
void FEMSolver::printMatlabReadContents(std::vector<double>& test)
{
std::cout << "test result vector is size = " << test.size() << std::endl;
int incr = test.size() / 5;
Expand All @@ -573,7 +582,7 @@ int FEMSolver::importRhsVectorFromFile(std::string filename,
Vector_h_CG& targetVector, bool verbose)
{
if (filename.empty()) {
string errMsg = "No matlab file provided for RHS (b) vector.";
std::string errMsg = "No matlab file provided for RHS (b) vector.";
errMsg += " Specify the file using argument at commandline using -b switch.";
std::cerr << errMsg << std::endl;
return -1;
Expand All @@ -597,7 +606,7 @@ int FEMSolver::importStiffnessMatrixFromFile(std::string filename,
Matrix_ell_h* targetMatrix, bool verbose)
{
if (filename.empty()) {
string errMsg = "No matlab file provided for stiffness matrix (A).";
std::string errMsg = "No matlab file provided for stiffness matrix (A).";
errMsg += " Specify the file using argument at commandline using -A switch.";
std::cerr << errMsg << std::endl;
return -1;
Expand All @@ -607,7 +616,7 @@ int FEMSolver::importStiffnessMatrixFromFile(std::string filename,
return -1;
}
if (verbose) {
string msg = "Finished reading stiffness matrix.";
std::string msg = "Finished reading stiffness matrix.";
std::cout << msg << std::endl;
}
return 0;
Expand Down
42 changes: 13 additions & 29 deletions src/core/include/FEMSolver.h → src/FEMSolver.h
View file
Original file line number Diff line number Diff line change
@@ -1,30 +1,14 @@
#ifndef __FEMSOLVER_H__
#define __FEMSOLVER_H__

#include <stdio.h>
#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <signal.h>
#include <exception>
#include <fstream>
#include <types.h>
#include <TriMesh.h>
#include <tetmesh.h>
#include <cutil.h>
#include <FEM/FEM2D.h>
#include <FEM/FEM3D.h>
#include <my_timer.h>
#include <amg.h>
#include <typeinfo>
#include <time.h>
#include <cmath>
#include "smoothers/smoother.h"
#include "cycles/cycle.h"
#include "amg_level.h"

#ifdef WIN32
#include <cstdlib>
#define srand48 srand
#endif
#include <vector>
#include "TriMesh.h"
#include "tetmesh.h"
#include "types.h"

/** The class that represents all of the available options for FEM */
class FEMSolver {
Expand All @@ -39,19 +23,19 @@ class FEMSolver {
};

public:
FEMSolver(std::string fname = "../src/test/test_data/sphere334",
FEMSolver(std::string fname = "../src/test/test_data/simple",
bool verbose = false);
virtual ~FEMSolver();
public:
void solveFEM(Matrix_ell_h* A_h, Vector_h_CG* x_h, Vector_h_CG* b_h);
void getMatrixFromMesh(Matrix_ell_h* A_h);
static int readMatlabSparseMatrix(const std::string &filename, Matrix_ell_h *A_h);
static int readMatlabNormalMatrix(const std::string &filename, vector<double> *A_h);
static int readMatlabNormalMatrix(const std::string &filename, std::vector<double> *A_h);
int writeMatlabArray(const std::string &filename, const Vector_h_CG &array);
void checkMatrixForValidContents(Matrix_ell_h* A_h);
void writeVTK(std::vector <float> values, std::string fname);
void printElementWithHeader(vector<double>& test, unsigned int index);
void printMatlabReadContents(vector<double>& test);
void printElementWithHeader(std::vector<double>& test, unsigned int index);
void printMatlabReadContents(std::vector<double>& test);
int importRhsVectorFromFile(std::string filename,
Vector_h_CG& targetVector, bool verbose);
int importStiffnessMatrixFromFile(std::string filename,
Expand All @@ -75,10 +59,10 @@ class FEMSolver {
int metisSize_; // max size of coarsest level
int partitionMaxSize_; // max size of of the partition
int aggregatorType_; // aggregator metis (0) mis (1)
ConvergenceType convergeType_; // the convergence tolerance algorithm <absolute|relative>
int convergeType_; // the convergence tolerance algorithm <absolute (0)|relative (1)>
double tolerance_; // the convergence tolerance
CycleType cycleType_; // the cycle algorithm <V|W|F|KCG|PCGV|PCGW|PCGF>
SolverType solverType_; // the solving algorithm <AMG|PCG>
int cycleType_; // the cycle algorithm <V (0) | W (1) | F (2) | K (3)>
int solverType_; // the solving algorithm <AMG (0) | PCG (1)>
double smootherWeight_; // the weight parameter used in a smoother
double proOmega_; // the weight parameter used in prolongator smoother
int device_; // the GPU device number to specify
Expand Down
4 changes: 2 additions & 2 deletions src/core/cuda/aggregator.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,9 @@
* Allocates selector based on passed in type
*********************************************/
template <class Matrix, class Vector>
Aggregator<Matrix, Vector>* Aggregator<Matrix, Vector>::allocate(FEMSolver *cfg)
Aggregator<Matrix, Vector>* Aggregator<Matrix, Vector>::allocate(int type)
{
if (cfg->aggregatorType_ == 0)
if (type == 0)
return new MIS_Aggregator < Matrix, Vector > ;
else
return new RandMIS_Aggregator < Matrix, Vector > ;
Expand Down
79 changes: 54 additions & 25 deletions src/core/cuda/amg.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,36 @@
#include <my_timer.h>

template<class Matrix, class Vector>
AMG<Matrix, Vector>::AMG(FEMSolver * cfg) : fine(0), cfg(cfg) {
AMG<Matrix, Vector>::AMG(bool verbose, int convergeType, int cycleType,
int solverType, double tolerance, int cycleIters, int maxIters,
int maxLevels, int topSize, double smootherWeight,
int preInnerIters, int postInnerIters, int postRelaxes,
int dsType, int metisSize, int partitionMaxSize, double proOmega,
int aggregatorType, TriMesh* triMesh, TetMesh* tetMesh) :
fine(0), verbose_(verbose),
convergeType_(convergeType == 0 ? ABSOLUTE_CONVERGENCE : RELATIVE_CONVERGENCE),
solverType_(solverType == 0 ? AMG_SOLVER : PCG_SOLVER),
tolerance_(tolerance), cycleIters_(cycleIters),
maxIters_(maxIters), maxLevels_(maxLevels), topSize_(topSize),
smootherWeight_(smootherWeight), preInnerIters_(preInnerIters),
postInnerIters_(postInnerIters), postRelaxes_(postRelaxes),
dsType_(dsType), metisSize_(metisSize), partitionMaxSize_(partitionMaxSize),
proOmega_(proOmega), aggregatorType_(aggregatorType),
triMesh_(triMesh), tetMesh_(tetMesh) {
switch (cycleType) {
case 0:
this->cycleType_ = V_CYCLE;
break;
case 1:
this->cycleType_ = W_CYCLE;
break;
case 2:
this->cycleType_ = F_CYCLE;
break;
case 3:
this->cycleType_ = K_CYCLE;
break;
}
}

template<class Matrix, class Vector>
Expand All @@ -26,9 +55,9 @@ bool AMG<Matrix, Vector>::converged(const Vector &r, ValueType &nrm)
// nrm = get_norm(r, norm);

nrm = cusp::blas::nrm2(r);
if (this->cfg->convergeType_ == ABSOLUTE_CONVERGENCE)
if (this->convergeType_ == ABSOLUTE_CONVERGENCE)
{
return nrm <= this->cfg->tolerance_;
return nrm <= this->tolerance_;
} else //if (convergence==RELATIVE)
{
if (initial_nrm == -1)
Expand All @@ -37,7 +66,7 @@ bool AMG<Matrix, Vector>::converged(const Vector &r, ValueType &nrm)
return false;
}
//if the norm has been reduced by the tolerance then return true
if (nrm / initial_nrm <= this->cfg->tolerance_)
if (nrm / initial_nrm <= this->tolerance_)
return true;
else
return false;
Expand All @@ -61,55 +90,55 @@ void AMG<Matrix, Vector>::setup(const Matrix_d &Acsr_d) {
level->level_id = 0;
level->nn = Acsr_d.num_rows;

level->m_meshPtr = this->cfg->triMesh_;
level->m_tetmeshPtr = this->cfg->tetMesh_;
level->m_meshPtr = this->triMesh_;
level->m_tetmeshPtr = this->tetMesh_;

if (this->cfg->verbose_) std::cout << "Entering AMG setup loop." << std::endl;
if (this->verbose_) std::cout << "Entering AMG setup loop." << std::endl;
while (true)
{
int N = level->A_d.num_rows;
if (N < this->cfg->topSize_ || num_levels >= this->cfg->maxLevels_)
if (N < this->topSize_ || num_levels >= this->maxLevels_)
{
coarsestlevel = num_levels - 1;
Matrix_h Atmp = level->A_d;
cusp::array2d<ValueType, cusp::host_memory> coarse_dense(Atmp);
LU = cusp::detail::lu_solver<ValueType, cusp::host_memory >(coarse_dense);
break;
}
if (this->cfg->verbose_) std::cout << "Finished with lu_solver." << std::endl;
if (this->verbose_) std::cout << "Finished with lu_solver." << std::endl;

level->next = AMG_Level<Matrix, Vector>::allocate(this);
if (this->cfg->verbose_) std::cout << "Finished with AMG_Level_allocate." << std::endl;
level->createNextLevel(this->cfg->verbose_);
if (this->cfg->verbose_) std::cout << "Finished with createNextLevel call." << std::endl;
if (this->verbose_) std::cout << "Finished with AMG_Level_allocate." << std::endl;
level->createNextLevel(this->verbose_);
if (this->verbose_) std::cout << "Finished with createNextLevel call." << std::endl;

if (level->level_id == 0)
{
Ahyb_d_CG = level->A_d;
}
if (this->cfg->verbose_) std::cout << "Copied A_d." << std::endl;
if (this->verbose_) std::cout << "Copied A_d." << std::endl;

level->setup(); //allocate smoother !! must be after createNextLevel since A_d is used
if (this->cfg->verbose_) std::cout << "level->setup." << std::endl;
if (this->verbose_) std::cout << "level->setup." << std::endl;

level->next->level_id = num_levels;
level->next->nn = level->nnout;
level->next->m_xadj_d = level->m_xadjout_d;
level->next->m_adjncy_d = level->m_adjncyout_d;
int nextN = level->next->A_d.num_rows;
if (this->cfg->verbose_) std::cout << "level->next finished" << std::endl;
if (this->verbose_) std::cout << "level->next finished" << std::endl;

//resize vectors
level->xc_d = Vector_d(nextN, -1);
level->bc_d = Vector_d(nextN, -1);
if (this->cfg->verbose_) std::cout << "resize vectors finished" << std::endl;
if (this->verbose_) std::cout << "resize vectors finished" << std::endl;

//advance to the next level
level = level->next;

//increment the level counter
num_levels++;
if (this->cfg->verbose_)
if (this->verbose_)
std::cout << "Looping with num_levels=" << num_levels << std::endl;
}

Expand All @@ -123,17 +152,17 @@ void AMG<Matrix, Vector>::solve_iteration(const Vector_d_CG &b, Vector_d_CG &x)
{
Vector_d b_d(b);
Vector_d x_d(x);
switch (this->cfg->solverType_)
switch (this->solverType_)
{
case AMG_SOLVER:
//perform a single cycle on the amg hierarchy
fine->cycle(this->cfg->cycleType_, b_d, x_d, this->cfg->verbose_);
fine->cycle(this->cycleType_, b_d, x_d, this->verbose_);
x = Vector_d_CG(x_d);
break;
case PCG_SOLVER:
//create a single CG cycle (this will run CG immediatly)
CG_Flex_Cycle<Matrix_h_CG, Vector_h_CG >(this->cfg->cycleType_, this->cfg->cycleIters_,
fine, Ahyb_d_CG, b, x, this->cfg->tolerance_, this->cfg->maxIters_, this->cfg->verbose_); //DHL
CG_Flex_Cycle<Matrix_h_CG, Vector_h_CG >(this->cycleType_, this->cycleIters_,
fine, Ahyb_d_CG, b, x, this->tolerance_, this->maxIters_, this->verbose_); //DHL
break;
}
}
Expand All @@ -145,12 +174,12 @@ void AMG<Matrix, Vector>::solve_iteration(const Vector_d_CG &b, Vector_d_CG &x)
template <class Matrix, class Vector>
void AMG<Matrix, Vector>::solve(const Vector_d_CG &b_d, Vector_d_CG &x_d)
{
if (this->cfg->verbose_)
if (this->verbose_)
printf("AMG Solve:\n");
iterations = 0;
initial_nrm = -1;

if (this->cfg->verbose_) {
if (this->verbose_) {
std::cout << std::setw(15) << "iter" << std::setw(15) << "time(s)" << std::setw(15)
<< "residual" << std::setw(15) << "rate" << std::setw(15) << std::endl;
std::cout << " ----------------------------------------------------\n";
Expand All @@ -163,8 +192,8 @@ void AMG<Matrix, Vector>::solve(const Vector_d_CG &b_d, Vector_d_CG &x_d)
solve_iteration(b_d, x_d);

done = true; //converged(r_d, nrm);
} while (++iterations < this->cfg->maxIters_ && !done);
if (this->cfg->verbose_)
} while (++iterations < this->maxIters_ && !done);
if (this->verbose_)
std::cout << " ----------------------------------------------------\n";

solve_stop = CLOCK();
Expand Down
Loading

0 comments on commit fa8ff04

Please sign in to comment.