Commit fe58c661 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

remove cython

parent f26e842f
all :
python2 setup.py build_ext --inplace
petsc :
python2 setup_petsc.py build_ext --inplace
clean :
python setup.py clean
libmbfluid.so : fluid_problem.c linear_system_banded_avx.c mesh.c mesh_find.c quadtree.c
${CC} $^ -o $@ -Wno-unused-function -O3 -march=native -mtune=native -lm -shared -fPIC
from ctypes import *
import numpy
import signal
import os
dir_path = os.path.dirname(os.path.realpath(__file__))
lib = CDLL(os.path.join(dir_path,"libmbfluid.so"))
signal.signal(signal.SIGINT, signal.SIG_DFL)
class StrongBoundary(Structure):
_fields_ = [("tag",c_char_p),("field",c_int),("v",c_double)]
class fluid_problem :
def __init__(self, mesh_file_name, mu, rho, epsilon, strong_boundaries):
nsb = len(strong_boundaries)
asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],i[2]) for i in strong_boundaries])
self._fp = lib.fluid_problem_new(mesh_file_name.encode(), c_double(mu), c_double(rho), c_double(epsilon), nsb, asb)
if self._fp == None :
raise NameError("cannot create fluid problem")
def __del__(self):
if(self._fp != None) :
lib.fluid_problem_free(self._fp)
def adapt_mesh(self, gradmin, gradmax, lcmin) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(gradmin), c_double(gradmax), c_double(lcmin))
def export(self, output_dir, t, it) :
lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))
def compute_node_force(self, dt) :
forces = numpy.ndarray([self.n_particles,2],"d")
ptr = lib.fluid_problem_compute_node_particle_force(self._fp, c_double(dt),forces.ctypes.data)
return forces
def implicit_euler(self, dt) :
lib.fluid_problem_implicit_euler(self._fp, c_double(dt))
def set_particles(self, mass, volume, position, velocity):
_mass = numpy.ascontiguousarray(mass).ctypes.data
_volume = numpy.ascontiguousarray(volume).ctypes.data
_position = numpy.ascontiguousarray(position).ctypes.data
_velocity = numpy.ascontiguousarray(velocity).ctypes.data
self.n_particles = mass.shape[0]
lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_mass,_volume,_position,_velocity)
import numpy
import signal
from cython.view cimport array as cvarray
signal.signal(signal.SIGINT, signal.SIG_DFL)
cdef extern from "linear_system.h":
ctypedef struct StrongBoundary :
const char *tag;
int field;
double v;
cdef extern from "fluid_problem.h":
ctypedef struct FluidProblem :
double mu
double epsilon
double rho
int n_particles
double *particle_force
double *node_force
pass
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int it);
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_boundary, StrongBoundary *boundaries);
void fluid_problem_free(FluidProblem *problem);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity);
cdef extern from "linear_system.h":
void initialize_linear_solver(int argc, char **argv);
include "openmpihack.pxi"
import sys
from libc.stdlib cimport malloc, free
argvdef = list([v for v in sys.argv])
if not "-pc_type" in argvdef :
argvdef += ["-pc_type", "lu"]
cdef char **argv = <char**>malloc(sizeof(char*)*(len(argvdef)));
argvbyte = [c.encode() for c in argvdef]
for i,c in enumerate(argvbyte) :
argv[i] = argvbyte[i]
initialize_linear_solver(len(argvdef), argv)
free(argv)
cdef class fluid_problem :
cdef FluidProblem *_fp
def __cinit__(self, mesh_file_name, mu, rho, epsilon, strong_boundaries):
cdef StrongBoundary* bnds = <StrongBoundary*> malloc(sizeof(StrongBoundary)*len(strong_boundaries));
tags = [bnd[0].encode() for bnd in strong_boundaries]
for i, bnd in enumerate(strong_boundaries) :
bnds[i].tag = <char*>tags[i]
bnds[i].field = bnd[1]
bnds[i].v = bnd[2]
self._fp = fluid_problem_new(mesh_file_name.encode(), mu, rho, epsilon, len(strong_boundaries),bnds)
free(bnds)
if self._fp == NULL :
raise NameError("cannot create fluid problem")
def __dealloc__(self):
if(self._fp != NULL) :
fluid_problem_free(self._fp)
def __getattr__(self, name) :
if name == "mu" :
return self._fp.mu
elif name == "rho" :
return self._fp.rho
elif name == "epsilon" :
return self._fp.epsilon
def __setattr__(self, name, val) :
if name == "mu" :
self._fp.mu = val
if name == "rho" :
self._fp.rho = val
elif name == "epsilon" :
self._fp.epsilon = val
def adapt_mesh(self, double gradmin, double gradmax, double lcmin) :
fluid_problem_adapt_mesh(self._fp, gradmin, gradmax, lcmin)
def export(self, output_dir, double t, int it) :
fluid_problem_export(self._fp, output_dir.encode(), t, it)
def compute_node_force(self, double dt) :
fluid_problem_compute_node_particle_force(self._fp, dt)
forces = numpy.ndarray([self._fp.n_particles,2],"d")
for i in range(self._fp.n_particles) :
forces[i,0] = self._fp.particle_force[i*2+0]
forces[i,1] = self._fp.particle_force[i*2+1]
return forces
def implicit_euler(self, double dt) :
fluid_problem_implicit_euler(self._fp, dt)
def set_particles(self, mass, volume, position, velocity):
cdef double[:,::1] _mass = numpy.ascontiguousarray(mass)
cdef double[:,::1] _volume = numpy.ascontiguousarray(volume)
cdef double[:,::1] _position = numpy.ascontiguousarray(position)
cdef double[:,::1] _velocity = numpy.ascontiguousarray(velocity)
fluid_problem_set_particles(self._fp,mass.shape[0],&_mass[0][0],&_volume[0][0],&_position[0][0],&_velocity[0][0])
......@@ -63,7 +63,7 @@ static void particle_drag(double u[2], double mu, double rho, double d, double c
drag[1] = GoU*u[1]/(1+(rhop+c/(1-c)*rho)*dt/(c/(1-c)*rho*mass)*GoU);
}
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt) {
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
double *porosity = problem->porosity;
Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
......@@ -115,8 +115,8 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt)
double d = 2*sqrt(vol/M_PI);
double drag[2];
particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop);
problem->particle_force[i*2+0] = -drag[0]-vol*gradp[0];
problem->particle_force[i*2+1] = -drag[1]-9.81*(mass-problem->rho*vol)-vol*gradp[1];
particle_force[i*2+0] = -drag[0]-vol*gradp[0];
particle_force[i*2+1] = -drag[1]-9.81*(mass-problem->rho*vol)-vol*gradp[1];
for (int iphi = 0; iphi < 3; ++iphi) {
problem->node_force[tri[iphi]*2+0] += drag[0]*phi[iphi];
problem->node_force[tri[iphi]*2+1] += drag[1]*phi[iphi];
......@@ -366,7 +366,6 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rh
problem->particle_element_id = malloc(0);
problem->particle_position = malloc(0);
problem->particle_mass = malloc(0);
problem->particle_force = malloc(0);
problem->particle_volume = malloc(0);
problem->particle_velocity = malloc(0);
return problem;
......@@ -383,7 +382,6 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->particle_element_id);
free(problem->particle_position);
free(problem->particle_mass);
free(problem->particle_force);
free(problem->particle_volume);
free(problem->particle_velocity);
for (int i = 0; i < problem->n_strong_boundaries; ++i)
......@@ -522,14 +520,12 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
free(problem->particle_element_id);
free(problem->particle_position);
free(problem->particle_mass);
free(problem->particle_force);
free(problem->particle_volume);
free(problem->particle_velocity);
problem->particle_uvw = malloc(sizeof(double)*2*n);
problem->particle_element_id = malloc(sizeof(int)*n);
problem->particle_position = malloc(sizeof(double)*2*n);
problem->particle_mass = malloc(sizeof(double)*n);
problem->particle_force = malloc(sizeof(double)*2*n);
problem->particle_volume = malloc(sizeof(double)*n);
problem->particle_velocity = malloc(sizeof(double)*n*2);
}
......
......@@ -23,13 +23,13 @@ typedef struct {
double *particle_volume;
double *particle_velocity;
double *particle_mass;
double *particle_force; // complete force on the particle (including gravity)
double *node_force;
int *particle_element_id;
} FluidProblem;
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter);
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt);
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity);
......
import os
from distutils.core import setup
from distutils.extension import Extension
from Cython.Build import cythonize
extensions = [Extension("fluid",["fluid.pyx",
"fluid_problem.c",
"linear_system_banded_avx.c",
#"linear_system_petsc.c",
"mesh.c",
"mesh_find.c",
"quadtree.c"
],
extra_compile_args = ["-Wno-unused-function", "-march=native", "-mtune=native"]
)]
setup(
ext_modules = cythonize(extensions)
)
......@@ -140,7 +140,7 @@ void ObjectCollection::display()
glPopMatrix();
}
}else {
int nTotPart = 4;
int nTotPart = 1;
for(size_t i = 0; i < radius.size()/nTotPart; ++i) {
glEnableClientState(GL_VERTEX_ARRAY);
glVertexPointer(2, GL_FLOAT, 0, &circleGeom[0]);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment