plasmaReactingFoam coupled with BE solver

This commit is contained in:
ignis 2016-04-15 12:46:53 +09:00
parent 33a4b4945c
commit 862314aac3
8 changed files with 7297 additions and 0 deletions

View file

@ -1,4 +1,5 @@
EXE_INC = \
-Ibolos \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
@ -12,6 +13,19 @@ EXE_INC = \
-I$(LIB_SRC)/combustionModels/lnInclude
EXE_LIBS = \
-Lbolos \
-L/usr/lib/python2.7/config-x86_64-linux-gnu \
-L/usr/lib \
-lpthread \
-ldl \
-lutil \
-lm \
-lpython2.7 \
-Xlinker \
-export-dynamic \
-Wl,-O1 \
-Wl,-Bsymbolic-functions \
-lbolos \
-lfiniteVolume \
-lfvOptions \
-lmeshTools \

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,28 @@
CC=g++
AR=ar
CFLAGS=-c -Wall -I/usr/include/python2.7 -I/usr/include/x86_64-linux-gnu/python2.7 -fno-strict-aliasing -D_FORTIFY_SOURCE=2 -g -fstack-protector --param=ssp-buffer-size=4 -Wformat -Werror=format-security -DNDEBUG -g -fwrapv -O0 -Wall -Wstrict-prototypes
LDFLAGS= -L/usr/lib/python2.7/config-x86_64-linux-gnu -L/usr/lib -lpthread -ldl -lutil -lm -lpython2.7 -Xlinker -export-dynamic -Wl,-O1 -Wl,-Bsymbolic-functions
HEADERS=bolos.h
SOURCES=pure_bolos.c
CPPSRCS=bolos.cpp
OBJECTS=$(SOURCES:.c=.o)
OBJECTS+=$(CPPSRCS:.cpp=.o)
EXECUTABLE=call_bolos
ARCHIVE=libbolos.a
all: $(HEADERS) $(SOURCES) $(CPPSRCS) $(EXECUTABLE) $(ARCHIVE)
$(ARCHIVE): $(CPPSRCS)
$(AR) -cvq $@ $(CPPSRCS:.cpp=.o)
$(EXECUTABLE): $(OBJECTS)
$(CC) $(OBJECTS) -o $@ $(LDFLAGS)
.c.o:
$(CC) $(CFLAGS) $< -o $@
.cpp.o:
$(CC) $(CFLAGS) $< -o $@
clean:
rm -rf *.o *.a

View file

@ -0,0 +1,410 @@
#include <Python.h>
#include <iostream>
#include "bolos.h"
using std::string;
int BE_IX::Bolos::solverCount = 0;
BE_IX::Bolos::PythonCAPI BE_IX::Bolos::PythonCAPI_;
BE_IX::Bolos::PythonCAPI::PythonCAPI()
:dummy(true)
{
Py_Initialize();
PyObject *pNameParser = PyString_FromString("bolos.parser");
PyObject *pNameSolver = PyString_FromString("bolos.solver");
PyObject *pNameGrid = PyString_FromString("bolos.grid");
PyObject *pModuleParser = PyImport_Import(pNameParser);
PyObject *pModuleSolver = PyImport_Import(pNameSolver);
PyObject *pModuleGrid = PyImport_Import(pNameGrid);
pyobjs_.insert ( str_pyobj("parser", pModuleParser) );
pyobjs_.insert ( str_pyobj("solver", pModuleSolver) );
pyobjs_.insert ( str_pyobj("grid", pModuleGrid) );
Py_DECREF(pNameParser);
Py_DECREF(pNameSolver);
Py_DECREF(pNameGrid);
}
BE_IX::Bolos::PythonCAPI::~PythonCAPI()
{
Py_Finalize();
}
BE_IX::Bolos::Bolos ()
: pyobjs_(),
x0(0.0),
x1(60.0),
nLGrid(200),
y0(0.),
y1(15.),
nQGrid(200),
kB(PyFloat_AsDouble(PyObject_GetAttrString
(PythonCAPI_.pyobjs()["solver"], "KB"))),
eV(PyFloat_AsDouble(PyObject_GetAttrString
(PythonCAPI_.pyobjs()["solver"], "ELECTRONVOLT"))),
Td(PyFloat_AsDouble(PyObject_GetAttrString
(PythonCAPI_.pyobjs()["solver"], "TOWNSEND"))),
f1(NULL)
{
// Import BOLOS python module
PyObject *pModuleParser = PythonCAPI_.pyobjs()["parser"];
PyObject *pModuleSolver = PythonCAPI_.pyobjs()["solver"];
PyObject *pModuleGrid = PythonCAPI_.pyobjs()["grid"];
pyobjs_.insert ( str_pyobj("parser", pModuleParser) );
pyobjs_.insert ( str_pyobj("solver", pModuleSolver) );
pyobjs_.insert ( str_pyobj("grid", pModuleGrid) );
PyObject *pClassLinGrid =
PyObject_GetAttrString(pModuleGrid, "LinearGrid");
pyobjs_.insert ( str_pyobj("LinearGrid", pClassLinGrid) );
PyObject *pClassBSolver =
PyObject_GetAttrString(pModuleSolver, "BoltzmannSolver");
pyobjs_.insert ( str_pyobj("BoltzmannSolver", pClassBSolver) );
PyObject *pNameLoadColls = PyString_FromString ("load_collisions");
pyobjs_.insert ( str_pyobj("load_collisions", pNameLoadColls) );
PyObject *pNameConverge = PyString_FromString("converge");
pyobjs_.insert ( str_pyobj("converge", pNameConverge) );
PyObject *pNameMeanEnergy = PyString_FromString("mean_energy");
pyobjs_.insert ( str_pyobj("mean_energy", pNameMeanEnergy) );
PyObject *pClassQuadGrid =
PyObject_GetAttrString(pModuleGrid, "QuadraticGrid");
pyobjs_.insert ( str_pyobj("QuadraticGrid", pClassQuadGrid) );
PyObject *pNameInterpolate = PyString_FromString("interpolate");
pyobjs_.insert ( str_pyobj("interpolate", pNameInterpolate) );
// Create Dummy grid to create BoltzmannSolver object
// gr = grid.LinearGrid(0, 60 ,200)
PyObject *pArgsLinGrid = Py_BuildValue ("(f, f, i)", 0.0, 60.0, 200);
PyObject *pGrid = PyObject_Call (pClassLinGrid, pArgsLinGrid, NULL);
if (!pClassLinGrid)
std::cout << "Test static initializer " << std::endl;
if (!pGrid)
std::cout << "pGrid is trouble" << std::endl;
// Create BoltzmannSolver object
// boltzmann = solver.BoltzmannSolver(gr)
PyObject *pArgsBSolver = PyTuple_New(1);
PyTuple_SetItem(pArgsBSolver, 0, pGrid);
bSolver =
PyObject_Call (pClassBSolver, pArgsBSolver, NULL);
pyobjs_.insert ( str_pyobj("boltzmann", bSolver) );
// Decrease reference count for garbage collection
Py_DECREF(pArgsLinGrid);
Py_DECREF(pArgsBSolver);
/*
KB = PyFloat_AsDouble (PyObject_GetAttrString(pModuleSolver, "KB"));
eV = PyFloat_AsDouble (PyObject_GetAttrString(pModuleSolver, "eV"));
Td = PyFloat_AsDouble (PyObject_GetAttrString(pModuleSolver, "Td"));
*/
}
bool BE_IX::Bolos::readCollisions (const string filename)
{
// with open ('./LXCat-June2013.txt') as fp:
PyObject *pModuleParser = pyobjs()["parser"];
PyObject *pFuncParse = PyObject_GetAttrString(pModuleParser, "parse");
PyObject *pFP = PyFile_FromString(const_cast<char*>(filename.c_str()), "r");
PyObject *pArgsParse = PyTuple_New(1);
PyTuple_SetItem(pArgsParse, 0, pFP);
// processes = parser.parse(fp)
PyObject *pProcesses = PyObject_Call (pFuncParse, pArgsParse, NULL);
pyobjs_.insert ( str_pyobj("processes", pProcesses) );
Py_DECREF(pFP);
Py_DECREF(pFuncParse);
Py_DECREF(pArgsParse);
return true;
}
void BE_IX::Bolos::setDensity (std::string specie, double n)
{
// boltzmann.target['N2'].density = 0.8
PyObject_CallMethod
(bSolver, "set_density", "(s, f)", specie.c_str(), n);
}
void BE_IX::Bolos::setTK (const double T)
{
// boltzmann.kT = 300 * solver.KB / solver.ELECTRONVOLT
PyObject_SetAttrString(bSolver, "kT", PyFloat_FromDouble(kB*T*eV));
}
void BE_IX::Bolos::setTeV (const double T)
{
// boltzmann.kT = 300 * solver.KB / solver.ELECTRONVOLT
PyObject_SetAttrString(bSolver, "kT", PyFloat_FromDouble(T*eV));
}
void BE_IX::Bolos::setEnTd (const double EN)
{
// boltzmann.EN = 120 * solver.TOWNSEND
PyObject_SetAttrString(bSolver, "EN", PyFloat_FromDouble(EN*Td));
}
void BE_IX::Bolos::setEnSI (const double EN)
{
// boltzmann.EN = 120 * solver.TOWNSEND
PyObject_SetAttrString(bSolver, "EN", PyFloat_FromDouble(EN));
}
void BE_IX::Bolos::presolve ()
{
PyObject *pBSolver = pyobjs_["boltzmann"];
PyObject *pProcesses = pyobjs_["processes"];
// boltzmann.load_collisions(processes)
PyObject *pNameLoadColls = pyobjs_["load_collisions"];
PyObject *pLoadReturn = PyObject_CallMethodObjArgs
(pBSolver, pNameLoadColls, pProcesses, NULL);
if (pLoadReturn) Py_DECREF(pLoadReturn);
PyObject *pModuleGrid = pyobjs_["grid"];
PyObject *pClassQuadGrid = pyobjs_["QuadraticGrid"];
PyObject *pArgsQuadGrid = Py_BuildValue ("(f, f, i)", 0., 60., 200);
PyObject_Print(pArgsQuadGrid, stdout, Py_PRINT_RAW); std::cout << std::endl;
PyObject *pNewGrid = PyObject_Call (pClassQuadGrid, pArgsQuadGrid, NULL);
PyObject_Print(pNewGrid, stdout, Py_PRINT_RAW); std::cout << std::endl;
PyObject_SetAttrString(pBSolver, "grid", pNewGrid);
// boltzmann.target['N2'].density = 0.8
setDensity("N2", 0.8);
// boltzmann.target['O2'].density = 0.2
setDensity("O2", 0.2);
// boltzmann.kT = 300 * solver.KB / solver.ELECTRONVOLT
setTK (300);
// boltzmann.EN = 120 * solver.TOWNSEND
setEnSI (120.0e-21);
// boltzmann.init()
PyObject *pInitReturn = PyObject_CallMethod(pBSolver, "init", NULL);
// fMaxwell = boltzmann.maxwell(2.0)
PyObject *fMax = PyObject_CallMethod(pBSolver, "maxwell", "(f)", 2.0);
PyObject_Print(fMax, stdout, Py_PRINT_RAW);
f1 = fMax;
PyObject_Print(f1, stdout, Py_PRINT_RAW);
}
void BE_IX::Bolos::solve (double T, double En, bool initMaxwell)
{
PyObject *pModuleSolver = pyobjs_["solver"];
PyObject *pModuleGrid = pyobjs_["grid"];
// boltzmann = solver.BoltzmannSolver(gr)
PyObject *pBSolver = pyobjs_["boltzmann"];
// boltzmann.target['N2'].density = 0.8
setDensity("N2", 0.8);
// boltzmann.target['O2'].density = 0.2
setDensity("O2", 0.2);
// boltzmann.kT = 300 * solver.KB / solver.ELECTRONVOLT
setTK (T);
// boltzmann.EN = 120 * solver.TOWNSEND
setEnSI (En);
// boltzmann.init()
PyObject *pInitReturn = PyObject_CallMethod(pBSolver, "init", NULL);
// f0 = boltzmann.converge(fMaxwell, maxn=100, rtol=1e-5)
PyObject *pNameConverge = pyobjs_["converge"];
PyObject *pIntMaxN = PyInt_FromLong(100);
PyObject *pFloatRTol = PyFloat_FromDouble(1e-5);
PyObject *pSolF1 =
PyObject_CallMethodObjArgs(pBSolver, pNameConverge, f1, pIntMaxN, pFloatRTol, NULL);
Py_DECREF(f1);
Py_DECREF(pIntMaxN);
Py_DECREF(pFloatRTol);
f1 = pSolF1;
return ;
}
void BE_IX::Bolos::maxwellian (double Te)
{
PyObject *pBSolver = pyobjs_["boltzmann"];
PyObject *fMax = PyObject_CallMethod(pBSolver, "maxwell", "(f)", Te);
Py_DECREF(f1);
f1 = fMax;
}
PyObject* BE_IX::Bolos::solve (double T, double En)
{
PyObject *pModuleSolver = pyobjs_["solver"];
PyObject *pModuleGrid = pyobjs_["grid"];
// gr = grid.LinearGrid(0, 60 ,200)
PyObject *pClassLinGrid = pyobjs_["LinearGrid"];
PyObject *pArgsLinGrid = Py_BuildValue ("(f, f, i)", x0, x1, nLGrid);
PyObject *pGrid = PyObject_Call (pClassLinGrid, pArgsLinGrid, NULL);
Py_DECREF(pArgsLinGrid);
// boltzmann = solver.BoltzmannSolver(gr)
PyObject *pBSolver = pyobjs_["boltzmann"];
PyObject_SetAttrString(pBSolver, "grid", pGrid);
// boltzmann.target['N2'].density = 0.8
setDensity("N2", 0.8);
// boltzmann.target['O2'].density = 0.2
setDensity("O2", 0.2);
// boltzmann.kT = 300 * solver.KB / solver.ELECTRONVOLT
setTK (T);
// boltzmann.EN = 120 * solver.TOWNSEND
setEnSI (En);
// boltzmann.init()
PyObject *pInitReturn = PyObject_CallMethod(pBSolver, "init", NULL);
// fMaxwell = boltzmann.maxwell(2.0)
PyObject *pSolFMaxwell = PyObject_CallMethod(pBSolver, "maxwell", "(f)", 2.0);
// f0 = boltzmann.converge(fMaxwell, maxn=100, rtol=1e-5)
PyObject *pNameConverge = pyobjs_["converge"];
PyObject *pIntMaxN = PyInt_FromLong(100);
PyObject *pFloatRTol = PyFloat_FromDouble(1e-5);
PyObject *pSolF0 =
PyObject_CallMethodObjArgs(pBSolver, pNameConverge, pSolFMaxwell, pIntMaxN, pFloatRTol, NULL);
// # Calculate the mean energy according to the first EEDF
// mean_energy = boltzmann.mean_energy(f0)
PyObject *pNameMeanEnergy = pyobjs_["mean_energy"];
PyObject *pMeanEnergy = PyObject_CallMethodObjArgs(pBSolver, pNameMeanEnergy, pSolF0, NULL);
double meanEnergy = PyFloat_AsDouble(pMeanEnergy);
// # Set a new grid extending up to 15 times the mean energy.
// # Now we use a quadritic grid instead of a linear one.
// newgrid = grid.QuadraticGrid(0, 15 * mean_energy, 200)
PyObject *pClassQuadGrid = pyobjs_["QuadraticGrid"];
PyObject *pArgsQuadGrid = Py_BuildValue ("(f, f, i)", y0, y1*meanEnergy, nQGrid);
PyObject *pNewGrid = PyObject_Call (pClassQuadGrid, pArgsQuadGrid, NULL);
// # Set the new grid and update the internal
// boltzmann.grid = newgrid
PyObject_SetAttrString(pBSolver, "grid", pNewGrid);
// boltzmann.init()
pInitReturn = PyObject_CallMethod(pBSolver, "init", NULL);
// # Calculate an EEDF in the new grid by interpolating the old one
// finterp = boltzmann.grid.interpolate(f0, gr)
PyObject *pBSolverGrid = PyObject_GetAttrString(pBSolver, "grid");
PyObject *pNameInterpolate = pyobjs_["interpolate"];
PyObject *pSolFInterp =
PyObject_CallMethodObjArgs(pBSolverGrid, pNameInterpolate, pSolF0, pGrid, NULL);
// # Iterate until we have a new solution
// f1 = boltzmann.converge(finterp, maxn=200, rtol=1e-5)
PyObject *pSolF1 = PyObject_CallMethodObjArgs
(pBSolver, pNameConverge, pSolFInterp, pIntMaxN, pFloatRTol, NULL);
Py_DECREF(f1);
Py_DECREF(pGrid);
Py_DECREF(pNewGrid);
Py_DECREF(pBSolverGrid);
Py_DECREF(pArgsLinGrid);
Py_DECREF(pArgsQuadGrid);
Py_DECREF(pIntMaxN);
Py_DECREF(pFloatRTol);
f1 = pSolF1;
return pSolF1;
}
double BE_IX::Bolos::diffusivity ()
const
{
/*
*/
// diffn = boltzmann.diffusion(f1)
PyObject *pNameDiffusivity = PyString_FromString("diffusion");
PyObject *pDiff_n = PyObject_CallMethodObjArgs
( bSolver, pNameDiffusivity, f1, NULL);
double diff_e = PyFloat_AsDouble (pDiff_n);
Py_DECREF(pNameDiffusivity);
Py_DECREF(pDiff_n);
return diff_e;
}
double BE_IX::Bolos::mobility ()
const
{
/*
*/
// mun = boltzmann.mobility(f1)
PyObject *pNameMobility = PyString_FromString("mobility");
PyObject *pMu_n =
PyObject_CallMethodObjArgs(bSolver, pNameMobility, f1, NULL);
double mu_e = PyFloat_AsDouble (pMu_n);
Py_DECREF(pNameMobility);
Py_DECREF(pMu_n);
return mu_e;
}
double BE_IX::Bolos::rateCoef (std::string process)
const
{
/*
*/
// k = boltzmann.rate(f1, "N2 -> N2^+")
PyObject *pNameRateCoef = PyString_FromString("rate");
PyObject *pNameProcess = PyString_FromString("N2 -> N2^+");
PyObject *pRateCoef = PyObject_CallMethodObjArgs
(bSolver, pNameRateCoef, f1, pNameProcess, NULL);
double k = PyFloat_AsDouble (pRateCoef);
Py_DECREF(pNameRateCoef);
Py_DECREF(pRateCoef);
Py_DECREF(pNameProcess);
return k;
}
BE_IX::Bolos::~Bolos ()
{
/*
if (solverCount > 1)
{
solverCount++;
}
else
{
Py_Finalize();
solverCount = 0;
}
*/
}

View file

@ -0,0 +1,94 @@
#ifndef BOLOS_H
#define BOLOS_H
#ifndef Py_PYTHON_H
struct PyObject;
#endif
#include <string>
#include <map>
namespace BE_IX
{
// typedef for std::map
typedef std::pair<std::string, PyObject*> str_pyobj;
class Bolos
{
protected:
class PythonCAPI
{
const bool dummy;
std::map<std::string, PyObject*> pyobjs_;
public:
PythonCAPI();
~PythonCAPI();
std::map<std::string, PyObject*>& pyobjs ()
{return pyobjs_;}
};
// Bolos object count for python api init and final
static int solverCount;
static PythonCAPI PythonCAPI_;
// Constants
//- 1 Td (TOWNSEND)
const double Td;
//- 1 eV (electron volt) in Joules
const double eV;
//- Boltzmann constant
const double kB;
// References to python objects
std::map<std::string, PyObject*> pyobjs_;
PyObject* lGrid;
PyObject* qGrid;
PyObject* bSolver;
// EEDF
PyObject* f1;
// Linear Grid Parameters
double x0;
double x1;
int nLGrid;
// Quadratic Grid Parameters
double y0;
double y1;
int nQGrid;
public:
Bolos ();
~Bolos ();
// member access
std::map<std::string, PyObject*>& pyobjs () {return pyobjs_;}
PyObject* EEDF () {return f1;}
// set processes to consider
bool readCollisions (std::string filename);
void setDensity (std::string specie, double n) ;
void setTK (const double T) ;
void setTeV (const double T) ;
void setEnTd (const double En) ;
void setEnSI (const double En) ;
// calculate EEDF
void presolve ();
PyObject* solve (double T, double En);
void solve (double T, double En, bool initMaxwell);
void maxwellian (double Te);
// calculate properties
double diffusivity () const ;
double mobility () const ;
double rateCoef (std::string process) const ;
};
};
#endif // Py_PYTHON_H

View file

@ -0,0 +1,67 @@
#include <Python.h>
#include "bolos.h"
#include <iostream>
int main(int argc, char *argv[])
{
BE_IX::Bolos bolos;
std::cout << "Success Bolos constructor " << std::endl;
// with open ('./LXCat-June2013.txt') as fp:
// processes = parser.parse(fp)
bolos.readCollisions("./LXCat-June2013.txt");
bolos.presolve();
std::cout << "Success Bolos presolve " << std::endl;
for (int i = 0; i < 100; i++)
{
bolos.solve(300.0, 120.0e-21, true);
}
PyObject *pBSolver = bolos.pyobjs()["boltzmann"];
std::cout << "Test new interface" << std::endl;
// mun = boltzmann.mobility(f1)
std::cout << "diffusivity " << bolos.diffusivity() << std::endl;
// diffn = boltzmann.diffusion(f1)
std::cout << "mobility " << bolos.mobility() << std::endl;
// k = boltzmann.rate(f1, "N2 -> N2^+")
std::cout << "rate coefficient "
<< bolos.rateCoef("N2 -> N2^+") << std::endl;
// for t, p in boltzmann.iter_all():
// print t, p
{
PyObject *iterator =
PyObject_CallMethod(pBSolver, "iter_all", "()");
PyObject *item;
if (iterator == NULL) {
/* propagate error */
std::cout << "iterator error" << std::endl;
}
while (item = PyIter_Next(iterator)) {
/* do something with item */
PyObject *pTarget = PyTuple_GetItem(item, 0);
PyObject *pProcess = PyTuple_GetItem(item, 1);
PyObject_Print(pTarget, stdout, Py_PRINT_RAW);
std::cout << std::endl;
PyObject_Print(pProcess, stdout, Py_PRINT_RAW);
std::cout << std::endl;
// release reference when done
Py_DECREF(pTarget);
Py_DECREF(pProcess);
Py_DECREF(item);
}
Py_DECREF(iterator);
}
return 0;
}

View file

@ -216,6 +216,22 @@ volScalarField Te
mesh
);
BE_IX::Bolos bolos;
bolos.readCollisions("./LXCat-June2013.txt");
bolos.presolve();
{
bolos.maxwellian(2);
forAll(rho, celli)
{
De[celli] = bolos.diffusivity();
mue[celli] = bolos.mobility();
}
De.correctBoundaryConditions();
mue.correctBoundaryConditions();
}
Info<< "Calculating face flux field ve\n" << endl;
surfaceScalarField ve
(

View file

@ -36,6 +36,8 @@ Description
#include "pimpleControl.H"
#include "fvIOoptionList.H"
#include "bolos.h"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])