Skip to content
Snippets Groups Projects
MpiSimulation.cpp 7.57 KiB
Newer Older
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
#include <iostream>
#include <cmath>
#include <cfloat>
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
#include <cstring>
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
#include <Body.hpp>
#include <Box.hpp>
#include <Tree.hpp>
#include <Node.hpp>
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
#include <iostream>
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
#include <BarnesHutTree.hpp>
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
#include "MpiSimulation.hpp"

namespace nbody {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
	using namespace std;

Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
	MpiSimulation::MpiSimulation(int& argc, char**& argv) {
		this->bodyType = MPI_DATATYPE_NULL;
		this->boxType = MPI_DATATYPE_NULL;

		int bodyBlocklengths[6] = {1, 3, 3, 3, 1, 1};
		MPI_Datatype bodyDatatypes[6] = {MPI_UNSIGNED_LONG, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT};
		MPI_Aint bodyOffsets[6];
	    bodyOffsets[0] = offsetof(Body, id);
	    bodyOffsets[1] = offsetof(Body, position);
	    bodyOffsets[2] = offsetof(Body, velocity);
	    bodyOffsets[3] = offsetof(Body, acceleration);
	    bodyOffsets[4] = offsetof(Body, mass);
	    bodyOffsets[5] = offsetof(Body, refinement);
	    MPI_Type_create_struct(6, bodyBlocklengths, bodyOffsets, bodyDatatypes, &this->bodyType);
	    MPI_Type_commit(&this->bodyType);

		int boxBlocklengths[2] = {3, 3};
		MPI_Datatype boxDatatypes[2] = {MPI_DOUBLE, MPI_DOUBLE};
		MPI_Aint boxOffsets[2];
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
	    boxOffsets[0] = offsetof(Box, min);
	    boxOffsets[1] = offsetof(Box, max);
	    MPI_Type_create_struct(2, boxBlocklengths, boxOffsets, boxDatatypes, &this->boxType);
	    MPI_Type_commit(&this->boxType);
		MPI_Comm_size(MPI_COMM_WORLD, &this->parallelSize);
		MPI_Comm_rank(MPI_COMM_WORLD, &this->parallelRank);
		this->domains = new Box[this->parallelSize];
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		this->correctState = true;
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		this->tree = new BarnesHutTree(this->parallelRank);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed

		if (argc == 2) {
			this->correctState = this->readInputData(string(argv[1]));
		} else {
			this->correctState = false;
		MPI_Bcast(&this->correctState, 1, MPI_INT, 0, MPI_COMM_WORLD);
		if (!this->correctState) {
			cerr << "Error occurred: terminating ..." << endl;
			MPI_Type_free(&this->bodyType);
			MPI_Type_free(&this->boxType);
			MPI_Finalize();
			exit(-1);
		}
	void MpiSimulation::cleanup() {
		MPI_Type_free(&this->bodyType);
		MPI_Type_free(&this->boxType);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		delete this->tree;
		this->tree = NULL;
	}

	bool MpiSimulation::stateCorrect() {
		return this->correctState != 0;
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
	}

	MpiSimulation::~MpiSimulation() {
		delete[] this->domains;
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		delete this->tree;
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		while (!this->sendStores.empty()) {
			delete[] this->sendStores.back().bodies;
			this->sendStores.pop_back();
		}
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
	}

	int MpiSimulation::getNumberOfProcesses() {
		return this->parallelSize;
	}

	int MpiSimulation::getProcessId() {
		return this->parallelRank;
	MPI_Datatype* MpiSimulation::getDatatype() {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		return &this->bodyType;
	}

Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
	void MpiSimulation::send(vector<Body> bodies, int target) {
		int bodySize = bodies.size();
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		SendStore* store = this->availableSendStore(bodySize);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		memcpy(store->bodies, &(bodies[0]), bodySize * sizeof(Body));
		MPI_Isend(store->bodies, bodySize, this->bodyType, target, 0, MPI_COMM_WORLD, &store->request);
		//MPI_Send(&(bodies[0]), bodySize, this->bodyType, target, 0, MPI_COMM_WORLD);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
	void MpiSimulation::recv(vector<Body>& bodies) {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		MPI_Status status;
		int count;
		Body* lb;

		MPI_Probe(0, 0, MPI_COMM_WORLD, &status);
		MPI_Get_count(&status, this->bodyType, &count);
		lb = new Body[count];
		MPI_Recv(lb, count, this->bodyType, status.MPI_SOURCE, 0, MPI_COMM_WORLD, &status);
		bodies = vector<Body>(lb, lb + count);
		delete[] lb;
	}

	void MpiSimulation::distributeBodies() {
		if (this->parallelRank == 0) {
			vector<Node> nodes;
			init(bb);
			nodes.push_back(Node(NULL));
			nodes.front().setBodies(this->bodies);
			extendForBodies(bb, this->bodies);
			nodes.front().setBB(bb);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed

			for (unsigned int i = 0; i < nodes.size(); i++) {
				printBB(i, nodes.front().getBB());
			}
			while (nodes.size() < (unsigned int) this->parallelSize) {
				int mostBodiesIndex = 0;

Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
				for (unsigned int i = 1; i < nodes.size(); i++) {
					if (nodes[i].getBodies().size() > nodes[mostBodiesIndex].getBodies().size()) {
						mostBodiesIndex = i;
					}
				}
				vector<Box> subdomains = splitLongestSide(nodes[mostBodiesIndex].getBB());
				vector<Body> buf = nodes[mostBodiesIndex].getBodies();
				Node n(NULL);

				n.setBodies(extractBodies(subdomains[0], buf));
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
				n.setBB(subdomains[0]);
				nodes.insert(nodes.begin() + mostBodiesIndex, n);
				n = Node(NULL);
				n.setBodies(extractBodies(subdomains[1], buf));
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
				n.setBB(subdomains[1]);
				nodes.insert(nodes.begin() + mostBodiesIndex, n);
				nodes.erase(nodes.begin() + mostBodiesIndex + 2);
			}
			this->bodies = nodes[0].getBodies();
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			for (unsigned int i = 1; i < nodes.size(); i++) {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
				this->send(nodes[i].getBodies(), i);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			}
		} else {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			this->recv(this->bodies);
	}

	void MpiSimulation::distributeDomains(vector<Body> localBodies) {
		Box localDomain;

Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		init(localDomain);
		extendForBodies(localDomain, localBodies);
		this->distributeDomains(localDomain);

	void MpiSimulation::distributeDomains() {
		this->distributeDomains(this->bodies);
	}

	void MpiSimulation::distributeDomains(Box localDomain) {
		this->domains[this->parallelRank] = localDomain;
		MPI_Allgather(&this->domains[this->parallelRank], 1, this->boxType, &this->domains[0], 1, this->boxType, MPI_COMM_WORLD);
		this->overallDomain = localDomain;
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		for (int i = 0; i < this->parallelSize; i++) {
			extend(this->overallDomain, this->domains[i]);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		}
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
	SendStore* MpiSimulation::availableSendStore(int numElems) {
		vector<SendStore>::iterator it = this->sendStores.begin();

		while (it != this->sendStores.end()) {
			int completed;

			MPI_Test(&it->request, &completed, MPI_STATUS_IGNORE);
			if (it->size >= numElems && completed) {
				return &(*it);
			} else if (completed) {
				delete[] it->bodies;
				it = this->sendStores.erase(it);
			} else {
				it++;
			}
		}
		SendStore store;
		store.bodies = new Body[numElems];
		store.size = numElems;
		this->sendStores.push_back(store);
		return &(this->sendStores.back());
	}

	void MpiSimulation::distributeLETs() {
		//send out locally essential trees (local bodies needed by remote simulations)
		for (int i = 0; i < this->parallelSize; i++) {
			if (i != this->parallelRank) {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
				//printBB(this->parallelRank, this->domains[i]);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
				vector<Body> refinements = this->tree->copyRefinements(this->domains[i]);

Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
				//cout << this->parallelRank << " -> " << i << ": " << refinements.size() << endl;
				this->send(refinements, i);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed

		//receive bodies and integrate them into local tree for simulation
		int received = 0;
		while (received < this->parallelSize - 1) {
			vector<Body> refinements;

Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			this->recv(refinements);
			//this->tree->mergeLET(refinements);
			/*
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			//integrate bodies in order of arrival to do communication/computation overlapping
			this->comms[0].recvBlocking(MPI_ANY_SOURCE, refinements);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			//cout << "recv: " << refinements.size() << endl;
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			this->tree->mergeLET(refinements);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			//this->tree.getRootBB().print();
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			*/
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		if (!this->tree->isCorrect()) {
			cout << "WRONG" << endl;
		}
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
	void MpiSimulation::buildTree() {
		if (this->bodies.size() > 0 && isValid(this->overallDomain)) {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			this->tree->build(this->bodies, this->overallDomain);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		} else if (this->bodies.size() > 0) {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			this->tree->build(this->bodies);
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		}
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		if (!this->tree->isCorrect()) {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			cout << "wrong tree" << endl;
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		} else {
			cout << this->parallelRank << " tree: " << this->tree->numberOfNodes() << endl;
	void MpiSimulation::rebuildTree() {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		this->tree->rebuild(this->overallDomain);
	void MpiSimulation::runStep() {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		this->distributeLETs();
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		/*
		this->tree->computeForces();
		this->distributeDomains(this->tree->advance());
		this->rebuildTree();
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		*/
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
		if (!this->tree->isCorrect()) {
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
			cout << "WRONG" << endl;
Paul Heinzlreiter's avatar
Paul Heinzlreiter committed
}