|
funkalicious 0.1
|
|
funkalicious 0.1
|
#include <libpostproc/lp_wavefunction.h>#include <libpostproc/density_1d.h>
Include dependency graph for fftw_wavefunc.h: This graph shows which files directly or indirectly include this file:Go to the source code of this file.
Functions | |
| 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 * | get_potential_for_density_inhomog_multigrid (const struct density *rho, const struct density *epsilon, double close_enough, GError **e) |
| 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_1d * | get_potential_for_density_1d (const struct density_1d *rho, const struct density_1d *epsilon, struct density_1d **rhoprime_out, GError **e) |
| 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:
1.7.3 
