Commit 1010bbf6 authored by Athena Economides's avatar Athena Economides
Browse files

Oscillating Stokes flow: run + post

parent 5f91911e
/* The Rayleigh Problem */
/* Flow over oscillating plane */
/* DOMAIN SETUP */
#define XS 32
......@@ -37,7 +37,7 @@
/* FLOW TYPE */
#define oscillating true
#define period 40
#define AA 100
#define AA 1
/* FEATURES */
#define rbcs false
......
......@@ -4,15 +4,16 @@ NX=1
NY=1
NZ=1
declare -a AA=(0.001 0.01 0.1 1 10 100)
#declare -a AA=(0.001 0.01 0.1 1 10 100)
declare -a AA=(1)
# Loop over all cases
for i in ${AA[@]}
do
echo "Setting case: AA=${i}"
BNAME=oscillating_plane_A${i}
BNAME=oscillating_long_AA${i}
sed -i "/define AA/ c\#define AA ${i}" ../conf.h
sed -i "/define AA/ c\#define AA ${i}" ../conf.h
. launch_daint.sh
done
......@@ -31,7 +31,7 @@ cd ${RUNDIR}
echo "#!/bin/bash -l
#
#SBATCH --job-name=osc_A${i}
#SBATCH --time=02:00:00
#SBATCH --time=06:00:00
#SBATCH --ntasks=${NN}
#SBATCH --ntasks-per-node=1
#SBATCH --output=slurm_out.%j.o
......
CXX = CC
profiles: main.cpp
$(CXX) main.cpp -I./include -g -O3 -Wno-deprecated-declarations -Wno-unused-result -o profiles
clean:
rm -f profiles
.PHONY = clean
/*
* ArgumentParser.h
* Cubism
*
*This argument parser assumes that all arguments are optional ie, each of the argument names is preceded by a '-'
*all arguments are however NOT optional to avoid a mess with default values and returned values when not found!
*
*More converter could be required:
*add as needed
*TypeName as{TypeName}() in Value
*
* Created by Christian Conti on 6/7/10. That is a long time ago.
* Modified by Diego Rossinelli several times after his dreadlocks hair cut.
* Copyright 2010 ETH Zurich. All rights reserved.
*
*/
#pragma once
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>
#include <map>
#include <vector>
#include <string>
#include <sstream>
#include <iterator>
using namespace std;
class Value
{
private:
string content;
public:
Value() : content("") {}
Value(string content_) : content(content_) { /*printf("%s\n",content.c_str());*/ }
double asDouble(double def=0) const
{
if (content == "") return def;
return (double) atof(content.c_str());
}
int asInt(int def=0) const
{
if (content == "") return def;
return atoi(content.c_str());
}
bool asBool(bool def=false) const
{
if (content == "") return def;
if (content == "0") return false;
if (content == "false") return false;
return true;
}
string asString(string def="") const
{
if (content == "") return def;
return content;
}
vector<float> asVecFloat(const int musthave_size = -1) const
{
//printf("mycontent is %s\n", content.c_str());
std::stringstream ss(content);
//assert(ss.good());
vector<float> retval;
double e;
while (ss >> e)
{
retval.push_back(e);
// printf("reading %f\n", e);
if (ss.peek() == ',')
ss.ignore();
}
if (musthave_size > 0)
assert(musthave_size == (int)retval.size());
return retval;
}
};
class ArgumentParser
{
private:
map<string,Value> mapArguments;
const int iArgC;
const char** vArgV;
bool bStrictMode, bVerbose;
const char delimiter;
public:
Value operator()(const string arg)
{
map<string,Value>::const_iterator it = mapArguments.find(arg);
if (bStrictMode)
{
if (it == mapArguments.end())
{
printf("Runtime option NOT SPECIFIED! ABORTING! name: %s\n",arg.data());
abort();
}
}
if (bVerbose)
printf("%s is %s\n", arg.data(), mapArguments[arg].asString().data());
if (it != mapArguments.end())
return mapArguments[arg];
else
return Value();
}
bool check(const string arg) const
{
return mapArguments.find(arg) != mapArguments.end();
}
ArgumentParser(const int argc, const char ** argv, bool bVerbose = false, const char delimiter = '=') :
mapArguments(), iArgC(argc), vArgV(argv), bStrictMode(false), bVerbose(bVerbose), delimiter(delimiter)
{
for (int i = 1; i<argc; i++)
{
int sep = 0;
while(argv[i][sep] != '\0' && argv[i][sep] != delimiter)
++sep;
string value;
if (argv[i][sep] != '\0')
value = string(argv[i] + sep + 1);
else
value = "1";
mapArguments[string(argv[i], sep)] = Value(value);
}
mute();
}
ArgumentParser(vector<string> args, bool bVerbose = false, const char delimiter = '='):
mapArguments(), iArgC(args.size()), vArgV(NULL), bStrictMode(false), bVerbose(bVerbose), delimiter(delimiter)
{
for(vector<string>::iterator it = args.begin(); it != args.end(); ++it)
{
const char * arg = it->c_str();
int sep = 0;
while(arg[sep] != '\0' && arg[sep] != delimiter)
++sep;
string value;
if (arg[sep] != '\0')
value = string(arg + sep + 1);
else
value = "1";
mapArguments[string(arg, sep)] = Value(value);
}
mute();
}
int getargc() const { return iArgC; }
const char** getargv() const { return vArgV; }
void set_strict_mode()
{
bStrictMode = true;
}
void unset_strict_mode()
{
bStrictMode = false;
}
void mute()
{
bVerbose = false;
}
void loud()
{
bVerbose = true;
}
void print_arguments(FILE * f = stdout)
{
printf("PRINTOUT OF THE RUNTIME OPTIONS\n");
for(map<string,Value>::const_iterator it=mapArguments.begin(); it!=mapArguments.end(); it++)
fprintf(f, "%s: <%s>\n", it->first.c_str(), it->second.asString().c_str());
printf("END OF THE PRINTOUT.\n");
}
void print_arguments(string path2log)
{
FILE * f = fopen(path2log.c_str(), "w");
if (f == NULL)
{
printf("could not save the log to <%s>. Exiting now\n", path2log.c_str());
exit(-1);
}
print_arguments(f);
fclose(f);
}
vector<string> find(string name)
{
map<string, Value>::iterator itb = mapArguments.lower_bound(name);
map<string, Value>::iterator ite = mapArguments.end();
vector<string> retval;
for(map<string, Value>::iterator it = itb; it != ite; ++it)
{
if (it->first.find(name) == string::npos)
break;
retval.push_back(it->first);
}
return retval;
}
};
#include <mpi.h>
#include <cstdio>
#define MPI_CHECK(ans) do { mpiAssert((ans), __FILE__, __LINE__); } while(0)
inline void mpiAssert(int code, const char *file, int line, bool abort=true)
{
if (code != MPI_SUCCESS)
{
char error_string[2048];
int length_of_error_string = sizeof(error_string);
MPI_Error_string(code, error_string, &length_of_error_string);
printf("mpiAssert: %s %d %s\n", file, line, error_string);
MPI_Abort(MPI_COMM_WORLD, code);
}
}
#include <mpi.h>
#include <unistd.h>
#include <fcntl.h>
#include <cstdio>
#include <iostream>
#include <vector>
#include <sstream>
#include <cmath>
#include <argument-parser.h>
#include <mpi-check.h>
using namespace std;
// run as: find ../../mpi-dpd/stress/* | ./stress -origin=0,0,0 -extent=48,48,48 -project=1,0,1 > profile.txt
int main(int argc, const char ** argv)
{
const int ndata=6; //depends on the outputs from io.cu:stress_dump function
MPI_CHECK( MPI_Init(&argc, (char ***)&argv) );
int nranks, rank;
MPI_CHECK( MPI_Comm_rank(MPI_COMM_WORLD, &rank));
MPI_CHECK( MPI_Comm_size(MPI_COMM_WORLD, &nranks));
ArgumentParser argp(argc, argv);
const bool verbose = argp("-verbose").asBool(false);
const bool avg = argp("-average").asBool(true);
vector<float> origin = argp("-origin").asVecFloat(3);
vector<float> extent = argp("-extent").asVecFloat(3);
vector<float> projectf = argp("-project").asVecFloat(3);
vector<float> mybinsize = argp("-binsize").asVecFloat(1);
string contributions = argp("-contributions").asString("uf");
const double ufactor = contributions.find("u") != string::npos; //=1 // string::npos=-1. It is better to compare to npos instead of -1 because the code is more legible
const double ffactor = contributions.find("f") != string::npos; //=1
bool project[3];
for(int c = 0; c < 3; ++c)
project[c] = (projectf[c] != 0); //=1 if dimension will be flattened
int nprojections = 0;
for(int c = 0; c < 3; ++c)
nprojections += project[c]; //number of flattened dimensions
const int noutputchannels = 4; //ux, uy, uz, d
const size_t chunksize = (1 << 29) / (ndata * sizeof(float));
float * const pbuf = new float[ndata * chunksize];
float binsize[3];
for(int c = 0; c < 3; ++c)
binsize[c] = project[c] ? extent[c] : mybinsize[0]; //if that dimension is flattened then binsize=extent, else binsize=mybinsize
int nbins[3];
for(int c = 0; c < 3; ++c)
nbins[c] = extent[c] / binsize[c]; //number of bins per dimension after flattening: =1 if dim will be flattened, =extent[c] if dim will NOT be flattened
const int ntotbins = nbins[0] * nbins[1] * nbins[2]; // total number of bins after flattening
int * const bincount = new int[ntotbins]; //array with size = ntotbins
memset(bincount, 0, sizeof(int) * ntotbins);
const int noutput = noutputchannels * ntotbins; //number of output data: ntotbins * number of data stored per bin
double * const bindata = new double[noutput]; //ptr to all the binned-data
memset(bindata, 0, sizeof(double) * noutput);
vector<string> paths; //each element contains a path to a data-file. elements are filled in parallel
{
string myinput;
if (rank == 0)
for (string line; getline(cin, line);)
myinput += line + "\n"; // = one stress-data-file in a new line
int inputsize = myinput.size(); //Returns the length of the string in bytes. 1 byte = 1 char
MPI_CHECK( MPI_Bcast(&inputsize, 1, MPI_INTEGER, 0, MPI_COMM_WORLD));
myinput.resize(inputsize); //Resizes the string to a length of n chars = n bytes. ..why is this thing here?? we need myinput to get the data-file paths
MPI_CHECK( MPI_Bcast(&myinput[0], inputsize, MPI_CHAR, 0, MPI_COMM_WORLD)); //now all mpi ranks to have the data-file paths
int c = 0;
istringstream iss(myinput);
for (string line; getline(iss, line); ++c) //transverse over myinput (data-file paths)
if (c % nranks == rank) //each rank pushes back a different data-file path
paths.push_back(line); //Save PATHS to stress dump files
}
int numfiles = paths.size(); //total number of data-file paths
size_t totalfootprint = 0; //will store sum of size of all files read by this rank. Only used to provide info to the user
double timeIO = 0;
// 1. Loop over data-files
for(int ipath = 0; ipath < (int)paths.size(); ++ipath)
{
const char * const path = paths[ipath].c_str(); //returns a c-style string, i.e. a null-terminated string
if (verbose)
fprintf(stderr, "working on <%s>\n", path);
int fdin = open(path, O_RDONLY); //if error:-1, else returns int>0: a file descriptor for the named file that is the lowest file descriptor not currently open for this process
if (!fdin)
{
fprintf(stderr, "can't access <%s> , exiting now.\n", path);
exit(-1);
}
if (verbose)
perror("reading...\n");
const size_t filesize = lseek(fdin, 0, SEEK_END); // size of the file described by fdin. lseek(..) return: the resulting offset location as measured in bytes from the beginning of the file
totalfootprint += filesize;
lseek(fdin, 0, SEEK_SET); //repositions the file offset at beginning of file
const size_t nparticles = filesize / ndata / sizeof(float); //number of particle per data-file
if (ipath%20==0)
fprintf(stderr, "Processing fileid %d/%d\n", ipath, (int)paths.size());
assert(filesize % (ndata * sizeof(float)) == 0);
if (verbose)
{
fprintf(stderr, "i have found %d particles\n", (int)nparticles);
fprintf(stderr, "particle chunk %d\n", (int)chunksize);
}
// 2. Loop over particles
for(size_t base = 0; base < nparticles; base += chunksize)
{
const int nhotparticles = min(nparticles - base, chunksize);
const size_t nhotbytes = nhotparticles * sizeof(float) * ndata;
size_t nreadbytes = 0;
int start = 0;
while(start < nhotparticles)
{
const double tstart = MPI_Wtime();
nreadbytes += read(fdin, pbuf, nhotbytes - nreadbytes); //read 3rd arg bytes from the file associated with the open file descriptor, 1st arg, into the buffer pointed to by 2nd arg
timeIO += MPI_Wtime() - tstart;
const int stop = nreadbytes / sizeof(float) / ndata;
#ifndef NDEBUG
if (verbose)
{
float avgs[ndata];
for(int i = 0; i < ndata; ++i)
avgs[i] = 0;
for(int i = 0; i < nhotparticles; ++i)
for(int c = 0; c < ndata; ++c)
avgs[c] += pbuf[ndata * i + c];
for(int i = 0; i < ndata; ++i)
printf("AVG %d: %.3e\n", i, avgs[i] / nhotparticles);
}
#endif
for(int i = start; i < stop; ++i)
{
const int srcbase = ndata * i;
int index[3];
for(int c = 0; c < 3; ++c)
index[c] = (int)((pbuf[srcbase + c] - origin[c]) / binsize[c]);
bool valid = true;
for(int c = 0; c < 3; ++c)
valid &= index[c] >= 0 && index[c] < nbins[c];
if (!valid)
continue;
const int binid = index[0] + nbins[0] * (index[1] + nbins[1] * index[2]);
++bincount[binid];
const int dstbase = noutputchannels * binid;
for(int c = 0; c < 3; ++c) //velocities (read +3, save +6)
bindata[dstbase + c] += pbuf[srcbase + 3 + c];
}
start = stop;
}
}
close(fdin);
}
if (rank == 0 && !numfiles)
{
perror("ooops zero files were read. Exiting now.\n");
exit(-1);
}
// sum bincount, bindata, timeIO and totalfootprint from all ranks
MPI_CHECK( MPI_Reduce(rank ? &numfiles : MPI_IN_PLACE, &numfiles, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD) );
MPI_CHECK( MPI_Reduce(rank ? bincount : MPI_IN_PLACE, bincount, ntotbins, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD) );
MPI_CHECK( MPI_Reduce(rank ? bindata : MPI_IN_PLACE, bindata, noutput, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD) );
MPI_CHECK( MPI_Reduce(rank ? &timeIO : MPI_IN_PLACE, &timeIO, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD) );
MPI_CHECK( MPI_Reduce(rank ? &totalfootprint : MPI_IN_PLACE, &totalfootprint, 1, MPI_OFFSET, MPI_SUM, 0, MPI_COMM_WORLD) );
if (rank) // all ranks finish here except for rank 0 which averages (by default) and then JUST PRINTS STUFF!!
goto finalize;
if (avg) //avg=true by default
for(int i = 0; i < ntotbins; ++i)
{
for(int c = 0; c < noutputchannels-1; ++c)
bindata[noutputchannels * i + c] /= bincount[i];
bindata[noutputchannels * i + 3] = bincount[i] == 0 ? 1./0. : (float)(bincount[i])/(binsize[0]*binsize[1]*binsize[2])/(float)numfiles;
}
if (nprojections == 3) // flatten all 3 dimensions: 0D data: keep just 1 point
{
assert(noutput == noutputchannels);
for(int c = 0; c < noutputchannels; ++c)
printf("%+.3e\t", bindata[c]);
printf("\n");
}
else if (nprojections == 2) // flatten 2 dimensions: 1D data: line
{
int ctr = 0;
for(int iz = 0; iz < nbins[2]; ++iz)
for(int iy = 0; iy < nbins[1]; ++iy)
for(int ix = 0; ix < nbins[0]; ++ix)
{
printf("%03d ", ctr);
for(int c = 0; c < noutputchannels; ++c)
printf("%+.4e ", bindata[noutputchannels * ctr + c]);
printf("\n");
++ctr;
}
}
else if (nprojections == 1) // flatten only 1 dimension: 2D data: surface
{
int nx = 0;
for(int c = 0; c < 3; ++c)
if (nbins[c] > 1)
{
nx = nbins[c];
break;
}
for(int c = 0; c < noutputchannels; ++c)
{
int ctr = 0;
for(int iz = 0; iz < nbins[2]; ++iz)
for(int iy = 0; iy < nbins[1]; ++iy)
for(int ix = 0; ix < nbins[0]; ++ix)
{
printf("%+.5e ", bindata[noutputchannels * ctr + c]);
++ctr;
if (ctr % nx == 0)
printf("\n");
}
if (c < noutputchannels - 1)
printf("\n");