|
funkalicious 0.1
|
|
funkalicious 0.1
|
#include <stdlib.h>#include <stdio.h>#include <fftw3.h>#include <math.h>#include <glib-2.0/glib.h>#include <libpostproc/fftw_wavefunc.h>
Include dependency graph for fftw_wavefunc.c:Go to the source code of this file.
Defines | |
| #define | FFTW_WAVEFUNC_GERROR_DOMAIN 32 |
| #define | COS cos |
Functions | |
| void | init_karray (fftw_complex *karray, int *ranges, double re, double im) |
| void | init_karray_1d (fftw_complex *karray, int range, double re, double im) |
| struct density * | get_potential_for_density (const struct density *rho, const struct density *epsilon, struct density **rhoprime_out, const struct density *dielectric_term_k, const struct density *dielectric_term_r, GError **e) |
| struct density_1d * | get_potential_for_density_1d (const struct density_1d *rho, const struct density_1d *epsilon, struct density_1d **rhoprime_out, GError **e) |
| void | set_del_phi_term (struct density *del_phi_term, const struct density *epsilon, const struct density *phi) |
| static void | set_del_phi_term_simple (struct density *del_phi, const struct density *epsilon, const struct density *phi, int periodic) |
| struct density * | get_potential_for_density_inhomog_simple (const struct density *rho, const struct density *epsilon, double close_enough, int global, int save_iterations, GError **e) |
| struct density * | get_potential_for_density_inhomog_multigrid (const struct density *rho, const struct density *epsilon, double close_enough, GError **e) |
| #define COS cos |
Definition at line 42 of file fftw_wavefunc.c.
Referenced by get_potential_for_density(), and get_potential_for_density_1d().
| #define FFTW_WAVEFUNC_GERROR_DOMAIN 32 |
Definition at line 36 of file fftw_wavefunc.c.
Referenced by get_potential_for_density(), get_potential_for_density_1d(), get_potential_for_density_inhomog_multigrid(), and get_potential_for_density_inhomog_simple().
| struct density* get_potential_for_density | ( | const struct density * | rho, |
| const struct density * | epsilon, | ||
| struct density ** | rhoprime_out, | ||
| const struct density * | dielectric_term_k, | ||
| const struct density * | dielectric_term_r, | ||
| GError ** | e | ||
| ) | [read] |
returns the electrostatic potential for a charge distribution rho and epsilon distribution epsilon.
| rho | a struct density containing the charge density. This should be in Coulombs per meter (the meter square would otherwise be multiplied internally; this increases accuracy to just do one division |
| epsilon | a struct density containing the electrostatic constant for the materials across the structure (*) |
| rhoprime_out | pointer to a pointer to a struct density to store rho divided by epsilon at each point. |
| dielectric_term_k,: | density which contains the *pre-fourier-transformed* inhomogeneous dielectric term (del ln(epsilon) . del phi) With this argument, this can be used for a step of an iterative or recursive approach. |
| dielectric_term_r,: | density which contains the realspace inhomogeneous dielectric term (del ln(epsilon) . del phi) With this argument, this can be used for a step of an iterative or recursive approach. |
Definition at line 76 of file fftw_wavefunc.c.
References ADD_EQUALS, COS, density_mul_by_const(), DIV_EQUALS, density::dx, FFTW_WAVEFUNC_GERROR_DOMAIN, free_density(), get_density_storage_at_c_(), density::N, new_density(), new_density_copy(), op_density_to_density(), and density::storage.
Referenced by get_potential_for_density_inhomog_multigrid(), get_potential_for_density_inhomog_simple(), and main().
{
for(int i=0; i<3; i++) {
if(rho->N[i] != epsilon->N[i]) {
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 100, "fftw_wavefunc.c::get_potential_for_density: charge density and dielectric constant grid sizes for dimension %d are not the same! (%d and %d, respectively)\n", i, rho->N[i], epsilon->N[i]);
return NULL;
}
}
/*If instructed to, intialize the FFTW threading*/
#ifdef OPENMP
int errval = fftw_init_threads();
if(errval == 0) {
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 100, "fftw_wavefunc.c::get_potential_for_density: Failed to initialize FFTW threading: error %d", errval);
return NULL;
}
int N_threads = omp_get_max_threads();
fftw_plan_with_nthreads(N_threads);
#endif
fprintf(stderr, "In get_potential_for_density: \n");
/*Construct rhoprime = rho/epsilon @ each lattice site*/
struct density* rhoprime = new_density_copy(rho);
if(rhoprime == NULL) {
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density: failed to allocate storage for rhoprime.");
return NULL;
}
fprintf(stderr, "\t allocated rhoprime\n");
if(op_density_to_density(rhoprime, DIV_EQUALS, epsilon)) {
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 101, "fftw_wavefunc.c::get_potential_for_density: failed to divide rho by epsilon at each point");
free_density(rhoprime);
return NULL;
}
fprintf(stderr, "\tdivided by epsilon\n");
/*Add in the realspace dielectric term if it's been specified.*/
if(dielectric_term_r != NULL) {
op_density_to_density(rhoprime, ADD_EQUALS, dielectric_term_r);
}
fprintf(stderr, "\tadded dielectric term\n");
/*Save a copy for output, if requested.*/
if(rhoprime_out != NULL) *rhoprime_out = rhoprime;
/*Centralizes the size of k-space arrays*/
int k_ranges[3] = {rhoprime->N[0]/2+1, rhoprime->N[1], rhoprime->N[2]};
//if((rhoprime->N[0])%2 == 1) k_ranges[0]++;
fftw_complex* rhoprime_k = (fftw_complex*)malloc(k_ranges[0]*k_ranges[1]*k_ranges[2]*sizeof(fftw_complex));
if(rhoprime_k == NULL) {
if(rhoprime_out == NULL) free_density(rhoprime);
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density: failed to allocate storage for k-space rhoprime.");
return NULL;
}
//fprintf(stderr, "\tinitializing rhoprime_k array with 0\n");
//init_karray(rhoprime_k, k_ranges, 0, 0);
/*Note that indices are in opposite order from how we've been using
them; is to make sure fftw gets the mapping into memory right.*/
fftw_plan p_f = fftw_plan_dft_r2c_3d(rhoprime->N[2], rhoprime->N[1], rhoprime->N[0], rhoprime->storage, rhoprime_k, FFTW_ESTIMATE);
fprintf(stderr, "\tExecuting FFT\n");
fftw_execute(p_f);
/*Free the plan now that we're done with it.*/
fftw_destroy_plan(p_f);
/*Calculate the potential in k-space*/
/*The "cosine term" is
*cos(2*pi*i/(N[0]/2))+cos(2*pi*j/N[1])+cos(2*pi*k/N[2])-4
*/
fftw_complex* phi_k = (fftw_complex*)malloc(k_ranges[0]*k_ranges[1]*k_ranges[2]*sizeof(fftw_complex));
if(phi_k == NULL ){
if(rhoprime_out == NULL) free_density(rhoprime);
free(rhoprime_k);
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density: failed to allocate storage for k-space potential.");
return NULL;
}
//fprintf(stderr, "\tinitializing phi_k with 0\n");
//init_karray(phi_k, k_ranges, 0, 0);
/*Algo is from Numerical Recipes, 3rd ed., p 1055*/
/*NOTE THAT DELTA^2 IS NO LONGER IN HERE; DIVIDE THE WAVEFUNCTION BY
A *SINGLE* GRID SIDE LENGTH BEFORE CALLING HERE!*/
int k, j, i;
fprintf(stderr, "\t setting phi_k\n");
#ifdef OPENMP
#pragma omp parallel for private(k, j, i)
#endif
for(k=0; k<k_ranges[2]; k++) {
for(j=0; j<k_ranges[1]; j++) {
for(i=0; i<k_ranges[0]; i++) {
phi_k[i + k_ranges[0]*(j + k_ranges[1]*k)][0]
/*TODO: fix this to not depend on dx->[0] only*/
= rhoprime_k[i + k_ranges[0]*(j + k_ranges[1]*k)][0]
/ (2.0*(COS(2.0*M_PI*((double)i)/((double)(rhoprime->N[0])))
+ COS(2.0*M_PI*((double)j)/((double)(rhoprime->N[1])))
+ COS(2.0*M_PI*((double)k)/((double)(rhoprime->N[2])))
- 3.0));
phi_k[i + k_ranges[0]*(j + k_ranges[1]*k)][1]
/*TODO: fix this to not depend on dx->[0] only*/
=rhoprime_k[i + k_ranges[0]*(j + k_ranges[1]*k)][1]
/ (2.0*(COS(2.0*M_PI*((double)i)/((double)(rhoprime->N[0])))
+ COS(2.0*M_PI*((double)j)/((double)(rhoprime->N[1])))
+ COS(2.0*M_PI*((double)k)/((double)(rhoprime->N[2])))
- 3.0));
if(dielectric_term_k != NULL)
/**\note the dielectric term is real*/
/**\note the d^2/dx^2 on the left-hand side introduces *dx^2, so we end up with one dx; the two dxes are implicitly cancelled in psi.*/
phi_k[i + k_ranges[0]*(j + k_ranges[1]*k)][0] += *get_density_storage_at_c_(i, j, k, dielectric_term_k);
#ifdef DEBUG3
fprintf(stderr, "i=%d j=%d k=%d -> %g + i%g\n", i, j, k, phi_k[i + k_ranges[0]*(j + k_ranges[1]*k)][0], phi_k[i + k_ranges[0]*(j + k_ranges[1]*k)][1]);
#endif
}
}
}
/*Set the 0th-component to 0; there's an inf and a nan otherwise;
rho is disconnected from phi_k at this point.*/
phi_k[0][0] = phi_k[0][1] = 0.0;
/*reverse plan*/
/*The final potential (density structure works as well for both)
will be saved in phi.*/
struct density* phi = new_density(rho->N, rho->dx, 0.0);
if(phi == NULL) {
if(rhoprime_out == NULL) free_density(rhoprime);
free(rhoprime_k);
free(phi_k);
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density: failed to allocate storage for realspace potential.");
return NULL;
}
fftw_plan p_r = fftw_plan_dft_c2r_3d(phi->N[2], phi->N[1], phi->N[0], phi_k, phi->storage, FFTW_ESTIMATE);
fprintf(stderr, "\treverse-transforming\n");
fftw_execute(p_r);
/*Free the plan now that we're done with it.*/
fftw_destroy_plan(p_r);
/*Free our structures*/
if(rhoprime == NULL) free_density(rhoprime);
free(rhoprime_k);
free(phi_k);
if(phi->storage == NULL) {
fprintf(stderr, "After reverse FFT, phi was NULL!\n");
}
/*After reverse-transforming, we need to scale the data; the
transforms are not normalized and must be scaled by 1/#real_elements
see http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data
*/
fprintf(stderr, "\tmultiplying by size of array\n");
density_mul_by_const(phi, -1.0/((phi->N[0])*(phi->N[1])*(phi->N[2])));
return phi;
}
Here is the call graph for this function: Here is the caller graph for this function:| struct density_1d* get_potential_for_density_1d | ( | const struct density_1d * | rho, |
| const struct density_1d * | epsilon, | ||
| struct density_1d ** | rhoprime_out, | ||
| GError ** | e | ||
| ) | [read] |
Definition at line 220 of file fftw_wavefunc.c.
References COS, density_1d_mul_by_const(), DIV_EQUALS, density_1d::dx, FFTW_WAVEFUNC_GERROR_DOMAIN, free_density_1d(), init_karray_1d(), density_1d::N, new_density_1d(), new_density_1d_copy(), op_density_1d_to_density_1d(), and density_1d::storage.
{
if(rho->N != epsilon->N) {
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 100, "fftw_wavefunc.c::get_potential_for_density_1d: charge density_1d and dielectric constant grid sizes for are not the same! (%d and %d, respectively)\n", rho->N, epsilon->N);
return NULL;
}
/*Construct rhoprime = rho/epsilon @ each lattice site*/
struct density_1d* rhoprime = new_density_1d_copy(rho);
if(rhoprime == NULL) {
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density_1d: failed to allocate storage for rhoprime");
return NULL;
}
if(op_density_1d_to_density_1d(rhoprime, DIV_EQUALS, epsilon)) {
free_density_1d(rhoprime);
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 101, "fftw_wavefunc.c::get_potential_for_density_1d: failed to divide rho by epsilon at each point");
return NULL;
}
/*Save a copy for output, if requested.*/
if(rhoprime_out != NULL) *rhoprime_out = rhoprime;
/*Centralizes the size of k-space arrays*/
int k_range = rhoprime->N/2 + 1;
if((rhoprime->N)%2 == 1) k_range++;
fftw_complex* rhoprime_k = (fftw_complex*)malloc(k_range*sizeof(fftw_complex));
if(rhoprime_k == NULL) {
if(rhoprime_out == NULL) free_density_1d(rhoprime);
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density_1d: failed to allocate storage for k-space rhoprime");
return NULL;
}
init_karray_1d(rhoprime_k, k_range, 0, 0);
/*Note that indices are in opposite order from how we've been using
them; is to make sure fftw gets the mapping into memory right.*/
fprintf(stderr, "rhoprime->N=%u krange=%d\n", rhoprime->N, k_range);
fftw_plan p_f = fftw_plan_dft_r2c_1d(rhoprime->N, rhoprime->storage, rhoprime_k, FFTW_ESTIMATE);
fftw_execute(p_f);
/*Calculate the potential in k-space*/
/*The "cosine term" is
*cos(2*pi*i/(N[0]/2))+cos(2*pi*j/N[1])+cos(2*pi*k/N[2])-4
*/
fftw_complex* phi_k = (fftw_complex*)malloc(k_range*sizeof(fftw_complex));
if(phi_k == NULL ){
if(rhoprime_out == NULL) free_density_1d(rhoprime);
free(rhoprime_k);
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density_1d: failed to allocate storage for k-space potential");
return NULL;
}
init_karray_1d(phi_k, k_range, 0, 0);
/*Algo is from Numerical Recipes, 3rd ed., p 1055*/
/*NOTE THAT DELTA^2 IS NO LONGER IN HERE; DIVIDE THE WAVEFUNCTION BY
A *SINGLE* GRID SIDE LENGTH BEFORE CALLING HERE!*/
for(int i=0; i<k_range; i++) {
phi_k[i][0]
/*TODO: fix this to not depend on dx->[0] only*/
= rhoprime_k[i][0]
/ (2.0*(COS(2.0*M_PI*((double)i)/((double)(rhoprime->N)))
- 1.0))*rhoprime->dx*rhoprime->dx;
phi_k[i][1]
/*TODO: fix this to not depend on dx->[0] only*/
= rhoprime_k[i][1]
/ (2.0*(COS(2.0*M_PI*((double)i)/((double)(rhoprime->N)))
- 1.0)) * rhoprime->dx*rhoprime->dx;
#ifdef DEBUG3
fprintf(stderr, "i=%d -> %g + i%g\n", i, phi_k[i][0], phi_k[i][1]);
#endif
}
/*Set the 0th-component to 0; there's an inf and a nan otherwise;
rho is disconnected from phi_k at this point.*/
//phi_k[0][0] = phi_k[0][1] = 0.0;
/*reverse plan*/
/*The final potential (density_1d structure works as well for both)
will be saved in phi.*/
struct density_1d* phi = new_density_1d(rho->N, rho->dx, 0.0);
if(phi == NULL) {
if(rhoprime_out == NULL) free_density_1d(rhoprime);
free(rhoprime_k);
free(phi_k);
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density_1d: failed to allocate storage for realspace potential");
return NULL;
}
fftw_plan p_r = fftw_plan_dft_c2r_1d(phi->N, phi_k, phi->storage, FFTW_ESTIMATE);
fftw_execute(p_r);
/*Free our structures*/
if(rhoprime == NULL) free_density_1d(rhoprime);
free(rhoprime_k);
free(phi_k);
if(phi->storage == NULL) {
fprintf(stderr, "After reverse FFT, phi was NULL!\n");
}
/*After reverse-transforming, we need to scale the data; the
transforms are not normalized and must be scaled by 1/#real_elements
see http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data
*/
density_1d_mul_by_const(phi, -1.0/(phi->N));
return phi;
}
Here is the call graph for this function:| struct density* get_potential_for_density_inhomog_multigrid | ( | const struct density * | rho, |
| const struct density * | epsilon, | ||
| double | close_enough, | ||
| GError ** | e | ||
| ) | [read] |
Solve the Poisson equation for an inhomogeneous dielectric using a multigrid method described in January 2011
| rho | charge dnesity in C/m (see comments on get_potential_for_density |
| epsilon | dielectric constant as a function of position |
| close_enough | convergence check (max abs percent change must be less than this value) |
| e | pointer to a pointer to a GError to return error information. |
Definition at line 630 of file fftw_wavefunc.c.
References density_get_max(), density::dx, FFTW_WAVEFUNC_GERROR_DOMAIN, free_density(), get_density_storage_at_c_(), get_double_resolution_density(), get_half_resolution_density(), get_potential_for_density(), density::N, new_density(), set_del_phi_term(), write_density_sweep(), and Z.
{
/*Get the higher-res grid versions; note that the density is LAPACK-based so the grid spacing is implicitly included in the value of the density (f*dx*dy*dz).*/
/*Get our initial phi*/
GError *err = NULL;
struct density* previous_phi = get_potential_for_density(rho, epsilon, NULL, NULL, NULL, &err);
if((previous_phi == NULL) || (err != NULL)) {
if(err != NULL) {
g_propagate_prefixed_error(e, err, "fftw_wavefunc.c::get_potential_for_density_inhomog_multigrid: error getting initial potential: ");
}else{
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 102, "fftw_wavefunc.c::get_potential_for_density_inhomog_multigrid: error getting initial potential: no further data because error was NULL; this is definitely a bug in the code! ");
}
if(previous_phi != NULL) free_density(previous_phi);
return NULL;
}
int level=0;
/*Dump the output*/
GString *iter_dir = g_string_new("del_phi_multigrid_iter=");
g_string_append_printf(iter_dir, "%d", level);
GFile* base_dir=g_file_new_for_path(iter_dir->str);
g_string_free(iter_dir, TRUE);
char* filename;
write_density_sweep(Z, previous_phi, &filename, &err, base_dir, FALSE);
err = NULL;
g_free(filename);
g_object_unref(base_dir);
struct density* hires_rho = NULL;
struct density* hires_epsilon = NULL;
struct density* del_phi_term = NULL;
struct density* phi = NULL;
double final_max_pct_diff;
do {
++level;
if(phi != NULL) {
/*To save space, start by removing the old del_phi term*/
free_density(del_phi_term);
/*Also, update previous_phi now*/
free_density(previous_phi);
previous_phi = phi;
/*And also to save space, only do one density at a time*/
struct density* hires_rho_ = get_double_resolution_density(rho, 3);
free_density(hires_rho);
hires_rho = hires_rho_;
struct density* hires_epsilon_ = get_double_resolution_density(epsilon, 0);
free_density(hires_epsilon);
hires_epsilon = hires_epsilon_;
}else{
hires_rho = get_double_resolution_density(hires_rho, 3);
hires_epsilon = get_double_resolution_density(hires_epsilon, 0);
}
/*Allocate and then set the del_phi term*/
del_phi_term = new_density(hires_rho->N, hires_rho->dx, 0);
set_del_phi_term(del_phi_term, epsilon, previous_phi);
/*Get the updated phi*/
struct density* phi = get_potential_for_density(hires_rho, hires_epsilon, NULL, del_phi_term, NULL, &err);
if((phi == NULL) || (err != NULL)) {
if(err != NULL) {
g_propagate_prefixed_error(e, err, "fftw_wavefunc.c::get_potential_for_density_inhomog_multigrid: error getting potential for iteration %d: ", level);
}else{
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 102, "fftw_wavefunc.c::get_potential_for_density_inhomog_multigrid: error getting initial potential for iteration %d: no further data because error was NULL; this is definitely a bug in the code! ", level);
}
if(phi != NULL) free_density(phi);
free_density(del_phi_term);
free_density(hires_rho);
free_density(hires_epsilon);
free_density(previous_phi);
return NULL;
}
iter_dir = g_string_new("del_phi_multigrid_iter=");
g_string_append_printf(iter_dir, "%d/phi", level);
base_dir=g_file_new_for_path(iter_dir->str);
g_string_free(iter_dir, TRUE);
err = NULL;
write_density_sweep(Z, phi, &filename, &err, base_dir, FALSE);
err = NULL;
g_free(filename);
g_object_unref(base_dir);
iter_dir = g_string_new("del_phi_multigrid_iter=");
g_string_append_printf(iter_dir, "%d/del_phi", level);
base_dir=g_file_new_for_path(iter_dir->str);
g_string_free(iter_dir, TRUE);
err = NULL;
write_density_sweep(Z, del_phi_term, &filename, &err, base_dir, FALSE);
err = NULL;
g_free(filename);
g_object_unref(base_dir);
/*parallelized finding maximum deviation from previous phi*/
int N_threads = 0;
#ifdef OPENMP
N_threads = omp_get_max_threads();
#else
N_threads = 1;
#endif
double max_pct_diff[N_threads];
for(int i=0; i<N_threads; i++) max_pct_diff[i] = nan("NaN");
long x, y, z;
double diff;
int thread;
double max_previous_phi = density_get_max(previous_phi, 1);
#ifdef OPENMP
#pragma omp parallel for private(y, x, diff, thread)
#endif
for(z=0; z<previous_phi->N[2]; z++) {
thread = 0;
#ifdef OPENMP
thread = omp_get_thread_num();
#endif
for(y=0; y<previous_phi->N[1]; y++) {
for(x=0; x<previous_phi->N[0]; x++) {
diff = abs(((*get_density_storage_at_c_(2*x, 2*y, 2*z, phi))-(*get_density_storage_at_c_(x, y, z, previous_phi)))/max_previous_phi);
//we use !<= since NaN will evaluate to false on any comparison
if(!(diff <= max_pct_diff[thread])) {
max_pct_diff[thread] = diff;
}
//fprintf(stderr, "fastmax@%d %ld %ld %ld: (%g - %g)/%g = %g (/->%g)\n", level, x, y, z, *get_density_storage_at_c_(2*x, 2*y, 2*z, phi), *get_density_storage_at_c_(x, y, z, previous_phi), max_previous_phi, diff, (*get_density_storage_at_c_(2*x, 2*y, 2*z, phi))/(*get_density_storage_at_c_(x, y, z, previous_phi)));
}
}
}
final_max_pct_diff = max_pct_diff[0];
for(int i=1; i<N_threads; i++) {
if(!(max_pct_diff[i] <= final_max_pct_diff)) final_max_pct_diff = max_pct_diff[i];
}
fprintf(stderr, "@%d: Max difference is %g:\n", level, final_max_pct_diff);
} while(final_max_pct_diff >= close_enough);
/*free the hires grids*/
free_density(hires_rho);
free_density(hires_epsilon);
free_density(del_phi_term);
/*Halve the grid resolution of phi until it's the original res*/
err = NULL;
for(int i=0; i<level; i++) {
struct density* phi_ = get_half_resolution_density(phi, &err, 0);
if(err != NULL) {
free_density(phi);
g_propagate_prefixed_error(e, err, "fftw_wavefunc::get_potential_for_density_inhomog@level %d: error occurred halving the resolution of a phi we got from a recursive call:\n\t", level);
return NULL;
}
free_density(phi);
phi = phi_;
}
return phi;
}
Here is the call graph for this function:| struct density* get_potential_for_density_inhomog_simple | ( | const struct density * | rho, |
| const struct density * | epsilon, | ||
| double | close_enough, | ||
| int | global, | ||
| int | save_iterations, | ||
| GError ** | e | ||
| ) | [read] |
Solve the Poisson equation for an inhomogeneous dielectric using simple symmetric finite differencing
| rho | charge dnesity in C/m (see comments on get_potential_for_density |
| epsilon | dielectric constant as a function of position |
| close_enough | convergence check (max abs percent change must be less than this value) |
| global | use maximum abs value of d1 instead of the local value of d1 for finding percent difference between iterations |
| save_iterations | save del_phi (z slices and grid dump) and phi (z slices) data from each iteration. |
| e | pointer to a pointer to a GError to return error information. |
Definition at line 546 of file fftw_wavefunc.c.
References density_get_max_abs_delta_pct(), density::dx, FFTW_WAVEFUNC_GERROR_DOMAIN, free_density(), get_potential_for_density(), density::N, new_density(), set_del_phi_term_simple(), write_density_slice(), write_density_sweep(), and Z.
Referenced by main().
{
fprintf(stderr, "get_potential_for_density_inhomog_simple: starting iteration.\n");
/*Get our initial phi*/
GError *err = NULL;
struct density* previous_phi = get_potential_for_density(rho, epsilon, NULL, NULL, NULL, &err);
if((previous_phi == NULL) || (err != NULL)) {
if(err != NULL) {
g_propagate_prefixed_error(e, err, "fftw_wavefunc.c::get_potential_for_density_inhomog_simple: error getting initial potential: ");
}else{
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 102, "fftw_wavefunc.c::get_potential_for_density_inhomog_simple: error getting initial potential: no further data because error was NULL; this is definitely a bug in the code! ");
}
if(previous_phi != NULL) free_density(previous_phi);
return NULL;
}
/*Dump the output*/
int level=0;
GString *iter_dir = g_string_new("del_phi_iter=0");
GFile* base_dir=g_file_new_for_path(iter_dir->str);
g_string_free(iter_dir, TRUE);
char* filename;
if(save_iterations) {
write_density_sweep(Z, previous_phi, &filename, &err, base_dir, FALSE);
g_free(filename);
}
err = NULL;
g_object_unref(base_dir);
struct density* phi = NULL;
struct density* del_phi = new_density(epsilon->N, epsilon->dx, 0);
double diff;
do {
++level;
set_del_phi_term_simple(del_phi, epsilon, previous_phi, 1);
phi = get_potential_for_density(rho, epsilon, NULL, NULL, del_phi, &err);
if((phi == NULL) || (err != NULL)) {
if(err != NULL) {
g_propagate_prefixed_error(e, err, "fftw_wavefunc.c::get_potential_for_density_inhomog_simple: error getting potential for iteration %d: ", level);
}else{
*e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 102, "fftw_wavefunc.c::get_potential_for_density_inhomog_simple: error getting initial potential for iteration %d: no further data because error was NULL; this is definitely a bug in the code! ", level);
}
if(phi != NULL) free_density(phi);
free_density(del_phi);
free_density(previous_phi);
return NULL;
}
iter_dir = g_string_new("del_phi_simple_iter=");
g_string_append_printf(iter_dir, "%d/phi", level);
base_dir=g_file_new_for_path(iter_dir->str);
g_string_free(iter_dir, TRUE);
err = NULL;
if(save_iterations) {
write_density_sweep(Z, phi, &filename, &err, base_dir, FALSE);
g_free(filename);
}
err = NULL;
g_object_unref(base_dir);
iter_dir = g_string_new("del_phi_simple_iter=");
g_string_append_printf(iter_dir, "%d/del_phi", level);
base_dir=g_file_new_for_path(iter_dir->str);
g_string_free(iter_dir, TRUE);
err = NULL;
if(save_iterations) {
write_density_sweep(Z, del_phi, &filename, &err, base_dir, FALSE);
g_free(filename);
}
err = NULL;
GFile* dumpfile = g_file_get_child(base_dir, "del_phi.dump");
if(save_iterations) write_density_slice(Z, NULL, del_phi, &filename, &err, dumpfile, ' ', "nm", "V");
err = NULL;
g_object_unref(base_dir);
fprintf(stderr, "@%d:\n", level);
if(global) {
diff = density_get_max_abs_delta_pct(previous_phi, phi, 1, 0);
}else{
diff = density_get_max_abs_delta_pct(previous_phi, phi, 0, 0);
}
fprintf(stderr, "\tdiff=%g\n", diff);
free_density(previous_phi);
previous_phi=phi;
} while(diff >= close_enough);
free_density(del_phi);
/**\note if we converged, then phi is contained in previous_phi from the end of the previous loop!!*/
return previous_phi;
}
Here is the call graph for this function: Here is the caller graph for this function:| void init_karray | ( | fftw_complex * | karray, |
| int * | ranges, | ||
| double | re, | ||
| double | im | ||
| ) |
Definition at line 46 of file fftw_wavefunc.c.
{
int i, j, k;
#ifdef OPENMP
#pragma omp parallel for private(i, j, k)
#endif
for(k=0; k<ranges[2]; k++) {
for(j=0; j<ranges[1]; j++) {
for(i=0; i<ranges[0]; i++) {
karray[i + ranges[0]*(j + ranges[1]*k)][0] = re;
karray[i + ranges[0]*(j + ranges[1]*k)][1] = im;
}
}
}
}
| void init_karray_1d | ( | fftw_complex * | karray, |
| int | range, | ||
| double | re, | ||
| double | im | ||
| ) |
Definition at line 61 of file fftw_wavefunc.c.
Referenced by get_potential_for_density_1d().
{
int i;
#ifdef OPENMP
#pragma omp parallel for private(i)
#endif
for(i=0; i<range; i++) {
karray[i][0] = re;
karray[i][1] = im;
}
}
Here is the caller graph for this function:| void set_del_phi_term | ( | struct density * | del_phi_term, |
| const struct density * | epsilon, | ||
| const struct density * | phi | ||
| ) |
sets the del_phi term
| del_phi_term | density struct to store the del_phi term; must be epsilon*N[i]*2-1 (i.e. doubled-grid size) |
| epsilon | dielectric distribution; must be same dimensions as phi |
| phi | electrostatic potential from the previous iteration |
even-z loop
even-y loop
The possible terms to average
Definition at line 325 of file fftw_wavefunc.c.
References density::dx, get_density_storage_index_(), density::N, and density::storage.
Referenced by get_potential_for_density_inhomog_multigrid().
{
/*First, set everything but the final x, y, and z faces*/
long x, y, z, ex, ox, ey, oy, ez, oz;
double phi_a, phi_b, phi_c, phi_d, phi_e, phi_f, phi_g, phi_h, epsilon_a, epsilon_b, epsilon_c, epsilon_d, epsilon_e, epsilon_f, epsilon_g, epsilon_h, AB, CD, DA, BC, FE, HG, EH, GF, BF, EA, DH, GC;
#ifdef OPENMP
#pragma omp parallel for private(y, x, ex, ox, ey, oy, ez, oz, phi_a, phi_b, phi_c, phi_d, phi_e, phi_f, phi_g, phi_h, epsilon_a, epsilon_b, epsilon_c, epsilon_d, epsilon_e, epsilon_f, epsilon_g, epsilon_h, AB, CD, DA, BC, FE, HG, EH, GF, BF, EA, DH, GC)
#endif
for(z=0; z<(epsilon->N[2]-1); ++z) {
ez=2*z;
oz=ez+1;
/**even-z loop*/
for(y=0; y<(epsilon->N[1]-1); ++y){
ey=2*y;
oy=ey+1;
/**even-y loop*/
for(x=0; x<(epsilon->N[0]-1); ++x){
/*Get phi and epsilon for the surrounding sites
top level: y+1 B C
y A D z
x x+1
bottom level: y+1 F G
y E H z+1
x x+1
*/
phi_a=phi->storage[get_density_storage_index_(x, y, z, phi->N[0], phi->N[1])];
phi_b=phi->storage[get_density_storage_index_(x, y+1, z, phi->N[0], phi->N[1])];
phi_c=phi->storage[get_density_storage_index_(x+1, y+1, z, phi->N[0], phi->N[1])];
phi_d=phi->storage[get_density_storage_index_(x+1, y, z, phi->N[0], phi->N[1])];
phi_e=phi->storage[get_density_storage_index_(x, y, z+1, phi->N[0], phi->N[1])];
phi_f=phi->storage[get_density_storage_index_(x, y+1, z+1, phi->N[0], phi->N[1])];
phi_g=phi->storage[get_density_storage_index_(x+1, y+1, z+1, phi->N[0], phi->N[1])];
phi_h=phi->storage[get_density_storage_index_(x+1, y, z+1, phi->N[0], phi->N[1])];
epsilon_a=epsilon->storage[get_density_storage_index_(x, y, z, epsilon->N[0], epsilon->N[1])];
epsilon_b=epsilon->storage[get_density_storage_index_(x, y+1, z, epsilon->N[0], epsilon->N[1])];
epsilon_c=epsilon->storage[get_density_storage_index_(x+1, y+1, z, epsilon->N[0], epsilon->N[1])];
epsilon_d=epsilon->storage[get_density_storage_index_(x+1, y, z, epsilon->N[0], epsilon->N[1])];
epsilon_e=epsilon->storage[get_density_storage_index_(x, y, z+1, epsilon->N[0], epsilon->N[1])];
epsilon_f=epsilon->storage[get_density_storage_index_(x, y+1, z+1, epsilon->N[0], epsilon->N[1])];
epsilon_g=epsilon->storage[get_density_storage_index_(x+1, y+1, z+1, epsilon->N[0], epsilon->N[1])];
epsilon_h=epsilon->storage[get_density_storage_index_(x+1, y, z+1, epsilon->N[0], epsilon->N[1])];
/**The possible terms to average*/
/*Top ones: clockwise*/
AB = log(epsilon_b/epsilon_a)*(phi_b - phi_a)/epsilon->dx[1];
CD = log(epsilon_d/epsilon_c)*(phi_d - phi_c)/epsilon->dx[1];
DA = log(epsilon_a/epsilon_d)*(phi_a - phi_d)/epsilon->dx[0];
BC = log(epsilon_c/epsilon_b)*(phi_c - phi_b)/epsilon->dx[0];
/*Bottom ones: counter clockwise*/
FE = log(epsilon_e/epsilon_f)*(phi_e - phi_f)/epsilon->dx[1];
HG = log(epsilon_g/epsilon_h)*(phi_g - phi_h)/epsilon->dx[1];
EH = log(epsilon_h/epsilon_e)*(phi_h - phi_e)/epsilon->dx[0];
GF = log(epsilon_f/epsilon_g)*(phi_f - phi_g)/epsilon->dx[0];
/*The z directions*/
BF = log(epsilon_f/epsilon_b)*(phi_f - phi_b)/epsilon->dx[2];
EA = log(epsilon_a/epsilon_e)*(phi_a - phi_e)/epsilon->dx[2];
DH = log(epsilon_h/epsilon_d)*(phi_h - phi_d)/epsilon->dx[2];
GC = log(epsilon_c/epsilon_g)*(phi_c - phi_g)/epsilon->dx[2];
fprintf(stderr, "del_phi_term: %ld %ld %ld: AB=%g CD=%g DA=%g BC=%g FE=%g HG=%g EH=%g GF=%g BF=%g EA=%g DH=%g GC=%g\n", x, y, z, AB, CD, DA, BC, FE, HG, EH, GF, BF, EA, DH, GC);
ex=2*x;
ox=ex+1;
/*We're setting the points from A*/
/*use oz oy ox: central point*/
del_phi_term->storage[get_density_storage_index_(ox, oy, oz, del_phi_term->N[0], del_phi_term->N[1])]
= (AB+CD+DA+BC+FE+HG+EH+GF+BF+EA+DH+GC)/12.0;
/*use oz oy ex: y-z facial*/
del_phi_term->storage[get_density_storage_index_(ex, oy, oz, del_phi_term->N[0], del_phi_term->N[1])]
= (AB + BF + FE + EA)/4.0;
/*use oz ey ox: x-z facial*/
del_phi_term->storage[get_density_storage_index_(ox, ey, oz, del_phi_term->N[0], del_phi_term->N[1])]
= (DA + EA + EH + DH)/4.0;
/*use oz ey ex: z axial */
del_phi_term->storage[get_density_storage_index_(ex, ey, oz, del_phi_term->N[0], del_phi_term->N[1])]
= EA;
/*use ez oy ox: x-y facial*/
del_phi_term->storage[get_density_storage_index_(ox, oy, ez, del_phi_term->N[0], del_phi_term->N[1])]
= (AB + BC + CD + DA)/4.0;
/*use ez oy ex: y axial*/
del_phi_term->storage[get_density_storage_index_(ex, oy, ez, del_phi_term->N[0], del_phi_term->N[1])]
= AB;
/*use ez ey ox: x axial*/
del_phi_term->storage[get_density_storage_index_(ex, ey, oz, del_phi_term->N[0], del_phi_term->N[1])]
= DA;
/*use ez ey ex: ignored*/
}
}
}
}
Here is the call graph for this function: Here is the caller graph for this function:| static void set_del_phi_term_simple | ( | struct density * | del_phi, |
| const struct density * | epsilon, | ||
| const struct density * | phi, | ||
| int | periodic | ||
| ) | [static] |
Sets the effective potential for the simple method.
| del_phi | density to be asigned the del_phi term. |
| epsilon | dielectric constant times epsilon_0 |
| phi | phi from the previous iteration |
| periodic | true if this is to be periodic boundary conditions. |
Definition at line 424 of file fftw_wavefunc.c.
References get_density_storage_at_(), get_density_storage_at_c_(), and density::N.
Referenced by get_potential_for_density_inhomog_simple().
{
/*Set the interior points with symmetric differencing*/
int x, y, z;
/**\note the dx terms cancel out (if all dx[0]==dx[1]==dx[2]) when multiplied by the dx^2 from the del^2 phi term, as is the case with rho!
double inverse_4_delta2_x = 1.0/(4.0*del_phi->dx[0]*del_phi->dx[0]);
double inverse_4_delta2_y = 1.0/(4.0*del_phi->dx[1]*del_phi->dx[1]);
double inverse_4_delta2_z = 1.0/(4.0*del_phi->dx[2]*del_phi->dx[2]);
*/
double inverse_4_delta2_x = 1.0/4.0;
double inverse_4_delta2_y = 1.0/4.0;
double inverse_4_delta2_z = 1.0/4.0;
#ifdef OPENMP
#pragma omp parallel for private(x, y, z)
#endif
for(z=1; z<del_phi->N[2]-1; z++) {
for(y=1; y<del_phi->N[1]-1; y++) {
for(x=1; x<del_phi->N[0]-1; x++) {
*get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
//fprintf(stderr, "%d %d %d: %g", x, y, z, *get_density_storage_at_c_(x, y, z, del_phi));
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
//fprintf(stderr, " %g", *get_density_storage_at_c_(x, y, z, del_phi));
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
//fprintf(stderr, " %g\n", *get_density_storage_at_c_(x, y, z, del_phi));
}
}
}
/*Set the z faces*/
z=0;
#ifdef OPENMP
#pragma omp parallel for private(x, y)
#endif
for(y=1; y<del_phi->N[1]-1; y++) {
for(x=1; x<del_phi->N[0]-1; x++) {
*get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
if(periodic) {
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, del_phi->N[2]-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, del_phi->N[2]-1, phi));
}else{
*get_density_storage_at_(x, y, z, del_phi) += 4.0*inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z, phi));
}
}
}
z=del_phi->N[2]-1;
#ifdef OPENMP
#pragma omp parallel for private(x, y)
#endif
for(y=1; y<del_phi->N[1]-1; y++) {
for(x=1; x<del_phi->N[0]-1; x++) {
*get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
if(periodic) {
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, 0, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, 0, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
}else{
*get_density_storage_at_(x, y, z, del_phi) += 4.0*inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
}
}
}
/*Set the y faces*/
y=0;
#ifdef OPENMP
#pragma omp parallel for private(x, z)
#endif
for(z=1; z<del_phi->N[2]-1; z++) {
for(x=1; x<del_phi->N[0]-1; x++) {
*get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
if(periodic) {
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, del_phi->N[1]-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, del_phi->N[1]-1, z, phi));
}else{
*get_density_storage_at_(x, y, z, del_phi) += 4.0*inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y, z, phi));
}
}
}
y=del_phi->N[1]-1;
#ifdef OPENMP
#pragma omp parallel for private(x, z)
#endif
for(z=1; z<del_phi->N[2]-1; z++) {
for(x=1; x<del_phi->N[0]-1; x++) {
*get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
if(periodic) {
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_y*log((*get_density_storage_at_c_(x, 0, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, 0, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
}else{
*get_density_storage_at_(x, y, z, del_phi) += 4.0*inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
}
}
}
/*Finally, the x faces*/
x=0;
#ifdef OPENMP
#pragma omp parallel for private(y, z)
#endif
for(z=1; z<del_phi->N[2]-1; z++) {
for(y=1; y<del_phi->N[1]-1; y++) {
*get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
if(periodic) {
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(del_phi->N[0]-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(del_phi->N[0]-1, y, z, phi));
}else{
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x, y, z, phi));
}
}
}
x=del_phi->N[0]-1;
#ifdef OPENMP
#pragma omp parallel for private(y, z)
#endif
for(z=1; z<del_phi->N[2]-1; z++) {
for(y=1; y<del_phi->N[1]-1; y++) {
*get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
if(periodic) {
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_x*log((*get_density_storage_at_c_(0, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(0, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
}else{
*get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_x*log((*get_density_storage_at_c_(x, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
}
}
}
}
Here is the call graph for this function: Here is the caller graph for this function:
1.7.3 
