Python GC для памяти/ссылки ctypes «удержание» ⇐ Python

Программы на Python
Ответить
Anonymous
 Python GC для памяти/ссылки ctypes «удержание»

Сообщение Anonymous »

Я написал код CFD, который я POC использовал на Python. Чтобы ускорить вычислительно затратные части, я переписал эти методы на C. Я вызываю их, используя ctypes, который, кажется, работает нормально, но иногда мне кажется, что сборщик мусора Python освобождает массив, который используется библиотекой C.
Исследуя проблему, мне сказали, что мне нужно хранить ссылки на объекты, выделенные в Python, а затем передавать их в библиотеку C, чтобы предотвратить сбор этих объектов сборщиком мусора, однако мне не совсем ясно, как правильно это делать. это. В некоторых ссылках говорится:
  • хранение ссылки на массив (numpy) (в Python),
  • некоторые говорят, что также необходимо хранить ссылку на указатель ctypes
    , который передается в библиотеку C.
В настоящее время я делаю что-то вроде этого:
lib = CDLL("./libc-compute")

# "Module" level holders to persist references for after python_method() returns
GLOBAL_BLOCK_HOLDERS = []
GLOBAL_R_HOLDERS = []
GLOBAL_V_HOLDERS = []

nvars = 5

def python_method(blocks):

R_blocks = {}
V_blocks = {}

GLOBAL_BLOCK_HOLDERS.clear() # So we dont simply accrue memory in python
GLOBAL_R_HOLDERS.clear()
GLOBAL_V_HOLDERS.clear()

blocks_to_process = list(blocks.keys())

if not use_c:
#use the slow python method
else
#use the compiled C library
for blockID in blocks_to_process:

block = blocks[blockID]
IM, JM, KM = block["IM"], block["JM"], block["KM"]
ncells = (IM - 1) * (JM - 1) * (KM - 1)

r_array = np.zeros(ncells * nvars, dtype=np.float64) # residuals: flattened with nvars
v_array = np.zeros(ncells, dtype=np.float64) # volumes: one per cell

block_holders = [] # create holders to store arrays from each struct
r_holders = []
v_holders = []
r_holders.append(r_array)
v_holders.append(v_array)

R_c = build_R_block(blockID, r_array, r_holders)
V_c = build_V_block(blockID, v_array, v_holders)
block_c = build_block(block, block_holders)

# Keep block_holder alive for the duration
GLOBAL_BLOCK_HOLDERS.append(block_holders)
GLOBAL_R_HOLDERS.append(r_holders)
GLOBAL_V_HOLDERSa.append(v_holders)

# Configure C function signature
lib.compute_flux_residuals.argtypes = [
POINTER(BLOCK), # blocks
POINTER(R_BLOCK),
POINTER(V_BLOCK)
]
lib.compute_flux_residuals.restype = None
lib.compute_flux_residuals(byref(block_c), byref(R_c), byref(V_c))

# Build NumPy views into the C-backed buffers (no copies).
# R: shape (IM-1, JM-1, KM-1, nvars)
R_flat = np.ctypeslib.as_array(R_c.r, shape=(ncells * nvars,))
R_blocks[blockID] = R_flat.reshape((IM - 1, JM - 1, KM - 1, nvars))

# V: shape (IM-1, JM-1, KM-1)
V_flat = np.ctypeslib.as_array(V_c.v, shape=(ncells,))
V_blocks[blockID] = V_flat.reshape((IM - 1, JM - 1, KM - 1))

# end loop over blocks

# Do a few other things to blocks, R_block / V_block
# in python only

# Return R_blocks, V_blocks for use in the outside main loop, blocks is not returned by also updated outside in the main loop
return R_blocks, V_blocks


Массивы C сглаживаются для скорости, поэтому я сглаживаю массивы numpy перед передачей в библиотеку. Любую память, которую я использую malloc() в C, я также освобождаю() в библиотеке. Единственная другая "память", которую я делаю в C, - это memcpy(V->v, block->V, ncells * sizeof(double)); для копирования объемов ячеек из структуры блока в массив V_blocks... в противном случае он просто обращается к сглаженным массивам numpy в C и обновляет значения...
Методы "строителя" структуры python ctype (т.е. build_R_block, build_block) добавляет каждый из массивов, включенных в структуры, в массив «держатель» (только r_array/v_array для массивов остатков и томов, несколько других для структуры «блока» данных). '
Большую часть времени код работает быстро и хорошо, но на больших сетках и в других нестандартных случаях я получаю ошибку сегментации из C (я думаю) (Process finished with exit code 139 (interrupted by signal 11:SIGSEGV)) ...
Мне было интересно, сможет ли кто-нибудь заметить что-нибудь явно неправильное, что я делаю в приведенном выше (упрощенном, но репрезентативном) виде, или вообще этот код подходит для целей «удержания ссылок». Я ценю любые указатели или советы. Заранее спасибо.
РЕДАКТИРОВАТЬ: Добавление функций C
// Allocate 4D flattened array: size = IM * JM * KM * nf
double* alloc4D(int IM, int JM, int KM, int nf) {
return calloc(IM * JM * KM * nf, sizeof(double));
}

void compute_flux_residuals(BLOCK *block, double gamma, enum RIEMANN riemann, double theta, double p_floor,
enum RECONSTRUCTION reconstruction, enum LIMITER limiter, R_BLOCK* R, V_BLOCK* V) {

printf("In C %d\n", block->id);

int IM = block->IM;
int JM = block->JM;
int KM = block->KM;

int ncells = (IM - 1) * (JM - 1) * (KM - 1);

double *rho = block->rho;
double *momx = block->momx;
double *momy = block->momy;
double *momz = block->momz;
double *rhoE = block->rhoE;
double *mu_eff = block->mu_eff;
double *a = block->a;

double * face_normals = block->n_hat;
double * face_areas = block->A;

// Copy volumes from block->V into V->v
memcpy(V->v, block->V, ncells * sizeof(double));

double *cx = block->cx;
double *cy = block->cy;
double *cz = block->cz;

const CELLOFFSET *nbr_offset = get_cell_offsets();

// Get the conserved vector
double *Uconserved = alloc4D(IM, JM, KM, NVARS);
get_conserved_vector(rho, momx, momy, momz, rhoE, IM, JM, KM, NVARS, Uconserved);

double UL[5];
double UR[5];
double UI[5];
// double UO[5];

for (int i = 0; i < IM - 1; i++) {
for (int j = 0; j < JM - 1; j++) {
for (int k = 0; k < KM - 1; k++) {

// Loop over all faces of cell (i,j,k)
for (int f = 0; f < 6; f++) {

double nx = face_normals[IDX5(i, j, k, f, 0, JM-1, KM-1, NFACES, NDIMS)];
double ny = face_normals[IDX5(i, j, k, f, 1, JM-1, KM-1, NFACES, NDIMS)];
double nz = face_normals[IDX5(i, j, k, f, 2, JM-1, KM-1, NFACES, NDIMS)];
double face_area = face_areas[IDX4(i,j,k, f, JM-1, KM-1, NFACES)];

CELLOFFSET offset = nbr_offset[f]; // get f-th offset

int di = offset.di;
int dj = offset.dj;
int dk = offset.dk;

int ni_i = i + di;
int ni_j = j + dj;
int ni_k = k + dk;

// Interior cell
if (ni_i >= 0 && ni_i < (IM - 1) &&
ni_j >= 0 && ni_j < (JM - 1) &&
ni_k >= 0 && ni_k < (KM - 1)) {

if (f == 0 || f == 2 || f == 4) {
continue; // skip this iteration
}

// Left side of face
UL[0] = Uconserved[IDX4(i, j, k, 0, JM-1, KM-1, NVARS)];
UL[1] = Uconserved[IDX4(i, j, k, 1, JM-1, KM-1, NVARS)];
UL[2] = Uconserved[IDX4(i, j, k, 2, JM-1, KM-1, NVARS)];
UL[3] = Uconserved[IDX4(i, j, k, 3, JM-1, KM-1, NVARS)];
UL[4] = Uconserved[IDX4(i, j, k, 4, JM-1, KM-1, NVARS)];

// Right side of face
UR[0] = Uconserved[IDX4(ni_i, ni_j, ni_k, 0, JM-1, KM-1, NVARS)];
UR[1] = Uconserved[IDX4(ni_i, ni_j, ni_k, 1, JM-1, KM-1, NVARS)];
UR[2] = Uconserved[IDX4(ni_i, ni_j, ni_k, 2, JM-1, KM-1, NVARS)];
UR[3] = Uconserved[IDX4(ni_i, ni_j, ni_k, 3, JM-1, KM-1, NVARS)];
UR[4] = Uconserved[IDX4(ni_i, ni_j, ni_k, 4, JM-1, KM-1, NVARS)];

double uxL = UL[1] / UL[0]; // Velocities left
double uyL = UL[2] / UL[0]; // Velocities left
double uzL = UL[3] / UL[0]; // Velocities left

double uxR = UR[1] / UR[0]; // Velocities right
double uyR = UR[2] / UR[0]; // Velocities right
double uzR = UR[3] / UR[0]; // Velocities right

double rhoL = UL[0]; // Density left
double rhoR = UR[0]; // Density right

double pL = pressure_from_conserved(UL, gamma, p_floor); // Pressure left
double pR = pressure_from_conserved(UR, gamma, p_floor); // Pressure right

double aL = a[IDX(i,j,k, JM-1, KM-1)]; // Speed of sound left
double aR = a[IDX(ni_i, ni_j, ni_k, JM-1, KM-1)]; // Speed of sound right

// Calculate the fluxes
double FL[5];
double FR[5];
double unL;
double unR;

// Calculate the inviscid fluxes
calculate_inviscid_fluxes(UL, rhoL, uxL, uyL, uzL, pL, nx, ny, nz, FL, &unL);
calculate_inviscid_fluxes(UR, rhoR, uxR, uyR, uzR, pR, nx, ny, nz, FR, &unR);

double lambda_max = fmax(fabs(unL) + aL, fabs(unR) + aR);

double Finviscid[5];
get_face_flux(UL, UR, FL, FR, lambda_max, Finviscid);
double muL = mu_eff[IDX(i, j, k, JM-1, KM-1)]; // Viscosity left
double muR = mu_eff[IDX(ni_i, ni_j, ni_k, JM-1, KM-1)]; // Viscosity right

double xL[] = {cx[IDX(i,j,k,JM-1,KM-1)], cy[IDX(i,j,k,JM-1,KM-1)], cz[IDX(i,j,k,JM-1,KM-1)]};
double xR[] = {cx[IDX(ni_i, ni_j, ni_k,JM-1,KM-1)], cy[IDX(ni_i, ni_j, ni_k,JM-1,KM-1)], cz[IDX(ni_i, ni_j, ni_k,JM-1,KM-1)]};

double dx = xR[0] - xL[0];
double dy = xR[1] - xL[1];
double dz = xR[2] - xL[2];
double d_n = fabs(dx * nx + dy * ny + dz * nz);

// Calculate the viscous fluxes
double Fviscous[5];
calculate_viscous_fluxes(uxL, uyL, uzL, uxR, uyR, uzR, muL, muR, nx, ny, nz, d_n, Fviscous);

// Calculate the total face flux
double Ftotal[5];
for (int n = 0; n < 5; n++)
Ftotal[n] = Finviscid[n] + Fviscous[n];

R->r[IDX4(i, j, k, 0, JM-1, KM-1, 5)] += Ftotal[0] * face_area;
R->r[IDX4(i, j, k, 1, JM-1, KM-1, 5)] += Ftotal[1] * face_area;
R->r[IDX4(i, j, k, 2, JM-1, KM-1, 5)] += Ftotal[2] * face_area;
R->r[IDX4(i, j, k, 3, JM-1, KM-1, 5)] += Ftotal[3] * face_area;
R->r[IDX4(i, j, k, 4, JM-1, KM-1, 5)] += Ftotal[4] * face_area;

R->r[IDX4(ni_i, ni_j, ni_k, 0, JM-1, KM-1, 5)] -= Ftotal[0] * face_area;
R->r[IDX4(ni_i, ni_j, ni_k, 1, JM-1, KM-1, 5)] -= Ftotal[1] * face_area;
R->r[IDX4(ni_i, ni_j, ni_k, 2, JM-1, KM-1, 5)] -= Ftotal[2] * face_area;
R->r[IDX4(ni_i, ni_j, ni_k, 3, JM-1, KM-1, 5)] -= Ftotal[3] * face_area;
R->r[IDX4(ni_i, ni_j, ni_k, 4, JM-1, KM-1, 5)] -= Ftotal[4] * face_area;
}
// Boundary cell
else {

UI[0] = Uconserved[IDX4(i, j, k, 0, JM-1, KM-1, NVARS)];
UI[1] = Uconserved[IDX4(i, j, k, 1, JM-1, KM-1, NVARS)];
UI[2] = Uconserved[IDX4(i, j, k, 2, JM-1, KM-1, NVARS)];
UI[3] = Uconserved[IDX4(i, j, k, 3, JM-1, KM-1, NVARS)];
UI[4] = Uconserved[IDX4(i, j, k, 4, JM-1, KM-1, NVARS)];

double rhoI = UI[0];
double pI = pressure_from_conserved(UI, gamma, p_floor);

double uI = UI[1] / UI[0];
double vI = UI[2] / UI[0];
double wI = UI[3] / UI[0];

double UR_temp[5];
if (!build_ghost_state_for_face(block, i, j, k, f, UI, gamma, UR_temp))
continue;

for (int i = 0; i < NVARS; i++) {
UR = UR_temp;
}

conserved_from_primitive(rhoI, uI, vI, wI, pI, NULL, gamma, UL);
double rhoL = UL[0];
double rhoR = UR[0];

double uL = UL[1] / rhoL;
double vL = UL[2] / rhoL;
double wL = UL[3] / rhoL;

double uR = UR[1] / rhoR;
double vR = UR[2] / rhoR;
double wR = UR[3] / rhoR;

double pL = pressure_from_conserved(UL, gamma, p_floor);
double pR = pressure_from_conserved(UR, gamma, p_floor);

double aL = sqrt(gamma * pL / rhoL);
double aR = sqrt(gamma * pR / rhoR);

// Calculate the fluxes
double FL[5];
double FR[5];
double unL;
double unR;

// Start with the inviscid fluxes
calculate_inviscid_fluxes(UL, rhoL, uL, vL, wL, pL, nx, ny, nz, FL, &unL);
calculate_inviscid_fluxes(UR, rhoR, uR, vR, wR, pR, nx, ny, nz, FR, &unR);

double lambda_max = fmax(fabs(unL) + aL, fabs(unR) + aR);

double Finviscid[5];
get_rusanov_face_flux(UL, UR, FL, FR, lambda_max, Finviscid);
R->r[IDX4(i, j, k, 0, JM-1, KM-1, NVARS)] += Finviscid[0] * face_area;
R->r[IDX4(i, j, k, 1, JM-1, KM-1, NVARS)] += Finviscid[1] * face_area;
R->r[IDX4(i, j, k, 2, JM-1, KM-1, NVARS)] += Finviscid[2] * face_area;
R->r[IDX4(i, j, k, 3, JM-1, KM-1, NVARS)] += Finviscid[3] * face_area;
R->r[IDX4(i, j, k, 4, JM-1, KM-1, NVARS)] += Finviscid[4] * face_area;
}
}
}
}
}

free(Uconserved);
}


Подробнее здесь: https://stackoverflow.com/questions/798 ... ce-holding
Вернуться к началу
Ответить

Быстрый ответ

Изменение регистра текста: 
Смайлики
:) :( :oops: :roll: :wink: :muza: :clever: :sorry: :angel: :read: *x)
Ещё смайлики…
   
К этому ответу прикреплено по крайней мере одно вложение.

Если вы не хотите добавлять вложения, оставьте поля пустыми.

Максимально разрешённый размер вложения: 15 МБ.

Вернуться в «Python»