#include <glib-2.0/glib.h>
#include <libpostproc/density_1d.h>

Include dependency graph for potential_1d.h:

This graph shows which files directly or indirectly include this file:

Enumerations
enum	density_1d_potential_calculation_method { INTEGRATION_METHOD_NP, INTEGRATION_METHOD_PN }
Functions
struct density_1d *	potential_1d_get_potential (double const_efield, double const_background, struct density_1d charge_density, struct density_1d epsilon, enum density_1d_potential_calculation_method method, double max_pot_int_err, double max_efield_int_err, double relerr, double abserr, GError **err)
struct density_1d *	potential_1d_get_dfield_integrate (double const_background, const struct density_1d charge_density, short int direction_pn, double max_int_err, double relerr, double abserr, GError **err)
struct density_1d *	potential_1d_get_potential_integrate (double dfield_background, double pot_background, const struct density_1d charge_density, const struct density_1d epsilon, short int direction_pn, double max_pot_int_err, double max_dfield_int_error, double relerr, double abserr, GError **err)
struct density_1d *	potential_1d_get_potential_from_dfield (double background, const struct density_1d Dfield, const struct density_1d epsilon, short int direction_pn, double max_int_err, double relerr, double abserr, GError *err)

Enumeration Type Documentation

enum density_1d_potential_calculation_method

Enumerator:

INTEGRATION_METHOD_NP
INTEGRATION_METHOD_PN

Definition at line 32 of file potential_1d.h.

                                               {
    INTEGRATION_METHOD_NP,
    INTEGRATION_METHOD_PN
  };

Function Documentation

struct density_1d* potential_1d_get_dfield_integrate	(	double	const_background,
		const struct density_1d *	charge_density,
		short int	direction_pn,
		double *	max_int_err,
		double	relerr,
		double	abserr,
		GError **	err
	)		`[read]`

Gets the 1D displacement field from a 1D charge density.

Parameters:

charge_density	the charge density (C/m^-3)
direction	0 to go negative-positive; 1 to go positive-negative
err	pointer to a pointer to a GError to store the routine's error information, if any

Returns:: density_1d* the potential, or NULL if an error occurred

err was set in the loop before calling errored

Definition at line 95 of file potential_1d.c.

References dfield_integrand(), density_1d::dx, e, free_density_1d(), MAX_GSL_INTEGRATION_INTERVALS, density_1d::N, new_density_1d(), POTENTIAL_1D_GERROR_DOMAIN, density_1d::storage, and str.

Referenced by potential_1d_get_potential_integrate().

                                                                                                                                                                                                               {

  struct density_1d* Efield = new_density_1d(charge_density->N, charge_density->dx, const_background);
  if(Efield == NULL) {
    *err = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 1, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_efield_integrate: Error allocating the new density_1d struct for the potential.");
    return NULL;
  }
  int n_threads = 1;
  #ifdef OPENMP
  n_threads = omp_get_max_threads();
  #endif
  double* x = (double*)malloc(sizeof(double)*charge_density->N);
  for(unsigned i=0; i<charge_density->N; i++) x[i] = i*charge_density->dx;
  gsl_interp_accel *charge_accel[n_threads];
  gsl_spline *charge[n_threads];
  struct dfield_integrand_params parms[n_threads];
  gsl_function gsl_func[n_threads];
  for(int i=0; i<n_threads; i++) {
    charge_accel[i] = gsl_interp_accel_alloc();
    charge[i] = gsl_spline_alloc(gsl_interp_akima, charge_density->N);
    gsl_spline_init(charge[i], x, charge_density->storage, charge_density->N);
    parms[i].charge = charge[i];
    parms[i].charge_accel = charge_accel[i];
#ifdef DEBUG_POTENTIAL_1D_C
    parms[i].call = 0;
    parms[i].file = NULL;
#endif
    gsl_func[i].function = &dfield_integrand;
    gsl_func[i].params = &(parms[i]);
  }
#ifdef DEBUG_POTENTIAL_1D_C
  FILE* foo = fopen("charge_density.pdata", "w");
  fprintf(foo, "#X\tx\tq\tq'\n");
  for(unsigned i=0; i<charge_density->N; i++) fprintf(foo, "%u\t%g\t%g\t%g\n", i, x[i], charge_density->storage[i], gsl_spline_eval(charge[0], x[i], charge_accel[0]));
  fclose(foo);
#endif
  int errored=0;
  gsl_integration_workspace** w = (gsl_integration_workspace**)malloc(sizeof(gsl_integration_workspace*)*n_threads);
  if(w == NULL) {
    *err = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 1, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_efield_integrate: Error allocating the workspaces for the integration.");
    return NULL;
  }
  for(int i=0; i<n_threads; i++) {
    w[i] = gsl_integration_workspace_alloc(MAX_GSL_INTEGRATION_INTERVALS);
  }
  int status;
  double abserr_ = abserr;
  double max_abs_ierr[n_threads];
  for(int i=0; i<n_threads; i++) max_abs_ierr[i]=-1;
  double max_ierr[n_threads];
  double ierr;
  if(direction_pn) {
    errored=0;
    int i;
    GString* str;
    #ifdef OPENMP
#pragma omp parallel for private(i, status, parms, ierr, str, abserr_)
    #endif
    for(i=charge_density->N-1; i>=0; i--) {
      int thread=1;
      #ifdef OPENMP
      thread = omp_get_thread_num();
      #endif
      abserr_ = abserr;
      do {
#ifdef DEBUG_POTENTIAL_1D_C
  str = g_string_new("charge_");
  g_string_append_printf(str, "thread=%d_i=%d_abserr=%g.pdata", thread, i, abserr_);
  parms[thread].call = 0;
  parms[thread].file=fopen(str->str, "w");
  g_string_free(str, TRUE);
#else
  str = NULL;
#endif
  status = gsl_integration_qag(&(gsl_func[thread]), (charge_density->N-1)*charge_density->dx, i*charge_density->dx, abserr_, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS61, w[thread], &(Efield->storage[i]), &ierr);
#ifdef DEBUG_POTENTIAL_1D_C
  fclose(parms[thread].file);
#endif
  if(abs(ierr) > max_abs_ierr[thread]) {
    max_ierr[thread] = ierr;
    max_abs_ierr[thread] = abs(ierr);
  }
  if(status == GSL_EROUND) abserr_ *= 10;
      }while(status == GSL_EROUND);
      if(status != 0) {
  #ifdef OPENMP
  #pragma omp critical
  #endif
  {
  errored = 1;
  *err = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 1, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_efield_integrate: Got error %d (%s) integrating from 0 to %g (grid site %u).", status, gsl_strerror(status), i*charge_density->dx, i);
  }
    #ifndef OPENMP
  break;
        #endif
      }
    }
    *max_int_err = -1;
    /*Find the maximum error from all the threads*/
    for(int i=0; i<n_threads; i++) if(max_abs_ierr[i] > abs(*max_int_err)) *max_int_err = max_ierr[i];
    if(errored) {
      /**err was set in the loop before calling errored*/
      free_density_1d(Efield);
      for(int i=0; i<n_threads; i++) {
  gsl_spline_free(charge[i]);
  gsl_interp_accel_free(charge_accel[i]);
  gsl_integration_workspace_free(w[i]);
      }
      return NULL;
    }
  }else{
    errored=0;
    int i;
    GString* str;
    #ifdef OPENMP
#pragma omp parallel for private(i, status, ierr, str, abserr_)
    #endif
    for(i=0; i<charge_density->N; i++) {
      int thread=1;
      #ifdef OPENMP
        thread = omp_get_thread_num();
      #endif
      abserr_ = abserr;
      do {
#ifdef DEBUG_POTENTIAL_1D_C
  str = g_string_new("charge_");
  g_string_append_printf(str, "thread=%d_i=%d_abserr=%g.pdata", thread, i, abserr_);
  parms[thread].call = 0;
  parms[thread].file=fopen(str->str, "w");
  g_string_free(str, TRUE);
#else
  str = NULL;
#endif
  status = gsl_integration_qag(&(gsl_func[thread]), 0, i*charge_density->dx, abserr_, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS61, w[thread], &(Efield->storage[i]), &ierr);
#ifdef DEBUG_POTENTIAL_1D_C
  fclose(parms[thread].file);
#endif
  if(abs(ierr) > max_abs_ierr[thread]) {
    max_ierr[thread] = ierr;
    max_abs_ierr[thread] = abs(ierr);
  }
  if(status == GSL_EROUND) abserr_ *= 10;
      }while(status == GSL_EROUND);
      if(status != 0) {
  #ifdef OPENMP
  #pragma omp critical
  #endif
  {
    errored = 1;
    if(*err == NULL) {
      *err = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 1, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_efield_integrate: Got error %d (%s) integrating from 0 to %g (grid site %u) (abserr_ is %g).", status, gsl_strerror(status), i*charge_density->dx, i, abserr_);
    }else{
      GError *e = NULL;
      g_propagate_prefixed_error(&e, *err, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_efield_integrate: Got error %d (%s) integrating from 0 to %g (grid site %u) (abserr_ is %g).\n", status, gsl_strerror(status), i*charge_density->dx, i, abserr_);
      *err = e;
    }
  }
    #ifndef OPENMP
  break;
        #endif
      }
    }
    *max_int_err = -1;
    /*Find the maximum error from all the threads*/
    for(int i=0; i<n_threads; i++) if(max_abs_ierr[i] > abs(*max_int_err)) *max_int_err = max_ierr[i];
    if(errored) {
      /**err was set in the loop before calling errored*/
      free_density_1d(Efield);
      for(int i=0; i<n_threads; i++) {
  gsl_spline_free(charge[i]);
  gsl_interp_accel_free(charge_accel[i]);
  gsl_integration_workspace_free(w[i]);
      }
      return NULL;
    }
  }
  for(int i=0; i<n_threads; i++) {
    gsl_spline_free(charge[i]);
    gsl_interp_accel_free(charge_accel[i]);
    gsl_integration_workspace_free(w[i]);
  }
  free(x);
  return Efield;
}

Here is the call graph for this function:

Here is the caller graph for this function:

struct density_1d* potential_1d_get_potential	(	double	const_efield,
		double	const_background,
		struct density_1d *	charge_density,
		struct density_1d *	epsilon,
		enum density_1d_potential_calculation_method	method,
		double *	max_pot_int_err,
		double *	max_efield_int_err,
		double	relerr,
		double	abserr,
		GError **	err
	)		`[read]`

Gets the 1D potential from a 1D charge density.

Parameters:

charge_density	the charge density (C/m^-3)
epsilon	dielectric constant across the structure
method	method to use to calculate the potential
err	pointer to a pointer to a GError to store the routine's error information, if any

Returns:: density_1d* the potential, or NULL if an error occurred

Definition at line 40 of file potential_1d.c.

References INTEGRATION_METHOD_NP, INTEGRATION_METHOD_PN, POTENTIAL_1D_GERROR_DOMAIN, and potential_1d_get_potential_integrate().

                                                                                                                                                                                                                                                                                                                 {
  switch(method) {
  case INTEGRATION_METHOD_NP:
    return potential_1d_get_potential_integrate(const_efield, const_background, charge_density, epsilon, 0, max_pot_int_err, max_efield_int_err, relerr, abserr, err);
    break;
  case INTEGRATION_METHOD_PN:
    return potential_1d_get_potential_integrate(const_efield, const_background, charge_density,epsilon, 1, max_pot_int_err, max_efield_int_err, relerr, abserr, err);
    break;
  default:
    *err = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 1, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_potential: no such potential calculation method: %d\n", method);
    return NULL;
  }
  *err = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 1, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_potential: Error allocating the new density_1d struct for the potential.");
  return NULL;
}

Here is the call graph for this function:

struct density_1d* potential_1d_get_potential_from_dfield	(	double	background,
		const struct density_1d *	Dfield,
		const struct density_1d *	epsilon,
		short int	direction_pn,
		double *	max_int_err,
		double	relerr,
		double	abserr,
		GError **	err
	)		`[read]`

Definition at line 306 of file potential_1d.c.

References density_1d::dx, free_density_1d(), MAX_GSL_INTEGRATION_INTERVALS, density_1d::N, new_density_1d(), POTENTIAL_1D_GERROR_DOMAIN, potential_integration_method_integrand(), density_1d::storage, and str.

Referenced by potential_1d_get_potential_integrate().

                                                                                                                                                                                                                                          {
  struct density_1d* phi = new_density_1d(Dfield->N, Dfield->dx, background);
  if(phi == NULL) {
    *err = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 1, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_potential: Error allocating the new density_1d struct for the potential.");
    return NULL;
  }
  int n_threads = 1;
  #ifdef OPENMP
  n_threads = omp_get_max_threads();
  #endif
  double* x = (double*)malloc(sizeof(double)*Dfield->N);
  for(unsigned i=0; i<Dfield->N; i++) x[i] = i*Dfield->dx;
  gsl_interp_accel *dfield_accel[n_threads];
  gsl_spline *dfield[n_threads];
  gsl_interp_accel *eps_accel[n_threads];
  gsl_spline *eps[n_threads];
  struct potential_integration_method_integrand_params parms[n_threads];
  gsl_function gsl_func[n_threads];
  for (int i=0; i<n_threads; i++) {
    dfield_accel[i] = gsl_interp_accel_alloc();
    dfield[i] = gsl_spline_alloc(gsl_interp_akima, Dfield->N);
    gsl_spline_init(dfield[i], x, Dfield->storage, Dfield->N);
    eps_accel[i] = gsl_interp_accel_alloc();
    eps[i] = gsl_spline_alloc(gsl_interp_akima, Dfield->N);
    gsl_spline_init(eps[i], x, epsilon->storage, epsilon->N);
    parms[i].dfield = dfield[i];
    parms[i].dfield_accel = dfield_accel[i];
    parms[i].epsilon = eps[i];
    parms[i].epsilon_accel = eps_accel[i];
#ifdef DEBUG_POTENTIAL_1D_C
    parms[i].call = 0;
    parms[i].file = NULL;
#endif
    gsl_func[i].function = &potential_integration_method_integrand;
    gsl_func[i].params = &(parms[i]);
  };
#ifdef DEBUG_POTENTIAL_1D_C
  FILE* foo = fopen("dfield.pdata", "w");
  fprintf(foo, "#X\tx\tq\tq'\teps\teps'\tq/eps\tq'/eps'\n");
  for(unsigned i=0; i<Dfield->N; i++) fprintf(foo, "%u\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", i, x[i], Dfield->storage[i], gsl_spline_eval(dfield[0], x[i], dfield_accel[0]), epsilon->storage[i], gsl_spline_eval(eps[0], x[i], eps_accel[0]), Dfield->storage[i]/epsilon->storage[i], gsl_spline_eval(dfield[0], x[i], dfield_accel[0])/gsl_spline_eval(eps[0], x[i], eps_accel[0]));
  fclose(foo);
#endif
  int errored=0;
  gsl_integration_workspace** w = (gsl_integration_workspace**)malloc(sizeof(gsl_integration_workspace*)*n_threads);
  if(w == NULL) {
    *err = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 1, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_efield_integrate: Error allocating the workspaces for the integration.");
    return NULL;
  }
  for(int i=0; i<n_threads; i++) {
    w[i] = gsl_integration_workspace_alloc(MAX_GSL_INTEGRATION_INTERVALS);
  }
  int status;
  double abserr_ = abserr;
  double max_abs_ierr[n_threads];
  for(int i=0; i<n_threads; i++) max_abs_ierr[i]=-1;
  double max_ierr[n_threads];
  double ierr;
  if(direction_pn) {
    errored=0;
    int i;
    GString* str;
    #ifdef OPENMP
#pragma omp parallel for private(i, status, ierr, str, abserr_)
    #endif
    for(i=Dfield->N-1; i>=0; i--) {
      int thread=1;
      #ifdef OPENMP
      thread = omp_get_thread_num();
      #endif
      abserr_ = abserr;
      do{
#ifdef DEBUG_POTENTIAL_1D_C
  str = g_string_new("charge_");
  g_string_append_printf(str, "thread=%d_i=%d_abserr=%g.pdata", thread, i, abserr_);
  parms[thread].call = 0;
  parms[thread].file=fopen(str->str, "w");
  g_string_free(str, TRUE);
#else
  str = NULL;
#endif
  status = gsl_integration_qag(&(gsl_func[thread]), (Dfield->N-1)*Dfield->dx, i*Dfield->dx, abserr_, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS61, w[thread], &(phi->storage[i]), &ierr);
#ifdef DEBUG_POTENTIAL_1D_C
  fclose(parms[thread].file);
#endif
  if(abs(ierr) > max_abs_ierr[thread]) {
    max_ierr[thread] = ierr;
    max_abs_ierr[thread] = abs(ierr);
  }
  if(status == GSL_EROUND) abserr_ *= 10;
      }while(status == GSL_EROUND);
      if(status != 0) {
  #ifdef OPENMP
  #pragma omp critical
  #endif
  {
  errored=1;
  *err = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 1, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_potential: Got error %d (%s) integrating from 0 to %g (grid site %u).", status, gsl_strerror(status), i*Dfield->dx, i);
  }
  #ifndef OPENMP
  break;
        #endif
      }
    }
    *max_int_err = -1;
    /*Find the maximum error from all the threads*/
    for(int i=0; i<n_threads; i++) if(max_abs_ierr[i] > abs(*max_int_err)) *max_int_err = max_ierr[i];
    if(errored) {
  free_density_1d(phi);
  for(int i=0; i<n_threads; i++) {
    gsl_spline_free(dfield[i]);
    gsl_interp_accel_free(dfield_accel[i]);
    gsl_spline_free(eps[i]);
    gsl_interp_accel_free(eps_accel[i]);
    gsl_integration_workspace_free(w[i]);
  }
  free(x);
  return NULL;
    }
  }else{
    errored=0;
    int i;
    GString* str;
    #ifdef OPENMP
#pragma omp parallel for private(i, status, ierr, str, abserr_)
    #endif
    for(i=0; i<Dfield->N; i++) {
      int thread=1;
      #ifdef OPENMP
      thread = omp_get_thread_num();
      #endif
      abserr_ = abserr;
      do{
#ifdef DEBUG_POTENTIAL_1D_C
  str = g_string_new("charge_");
  g_string_append_printf(str, "thread=%d_i=%d_abserr=%g.pdata", thread, i, abserr_);
  parms[thread].call = 0;
  parms[thread].file=fopen(str->str, "w");
  g_string_free(str, TRUE);
#else
  str = NULL;
#endif
  status = gsl_integration_qag(&(gsl_func[thread]), 0, i*Dfield->dx, abserr_, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS61, w[thread], &(phi->storage[i]), &ierr);
#ifdef DEBUG_POTENTIAL_1D_C
  fclose(parms[thread].file);
#endif
  if(abs(ierr) > max_abs_ierr[thread]) {
    max_ierr[thread] = ierr;
    max_abs_ierr[thread] = abs(ierr);
  }
  if(status == GSL_EROUND) abserr_ *= 10;
      }while(status == GSL_EROUND);
      if(status != 0) {
  #ifdef OPENMP
  #pragma omp critical
  #endif
  {
  errored=1;
  *err = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 1, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_potential: Got error %d (%s) integrating from 0 to %g (grid site %u).", status, gsl_strerror(status), i*Dfield->dx, i);
  }
  #ifndef OPENMP
  break;
        #endif
      }
    }
    *max_int_err = -1;
    /*Find the maximum error from all the threads*/
    for(int i=0; i<n_threads; i++) if(max_abs_ierr[i] > abs(*max_int_err)) *max_int_err = max_ierr[i];
    if(errored) {
  free_density_1d(phi);
  for(int i=0; i<n_threads; i++) {
    gsl_spline_free(dfield[i]);
    gsl_interp_accel_free(dfield_accel[i]);
    gsl_spline_free(eps[i]);
    gsl_interp_accel_free(eps_accel[i]);
    gsl_integration_workspace_free(w[i]);
  }
  free(x);
  return NULL;
    }
  }
  for(int i=0; i<n_threads; i++) {
    gsl_spline_free(dfield[i]);
    gsl_interp_accel_free(dfield_accel[i]);
    gsl_spline_free(eps[i]);
    gsl_interp_accel_free(eps_accel[i]);
    gsl_integration_workspace_free(w[i]);
  }
  free(x);
  return phi;
}

Here is the call graph for this function:

Here is the caller graph for this function:

struct density_1d* potential_1d_get_potential_integrate	(	double	dfield_background,
		double	pot_background,
		const struct density_1d *	charge_density,
		const struct density_1d *	epsilon,
		short int	direction_pn,
		double *	max_pot_int_err,
		double *	max_dfield_int_error,
		double	relerr,
		double	abserr,
		GError **	err
	)		`[read]`

Definition at line 280 of file potential_1d.c.

References density_1d_add_const(), e, free_density_1d(), density_1d::N, POTENTIAL_1D_GERROR_DOMAIN, potential_1d_get_dfield_integrate(), potential_1d_get_potential_from_dfield(), and density_1d::storage.

Referenced by potential_1d_get_potential().

                                                                                                                                                                                                                                                                                                               {
  GError *e = NULL;
  struct density_1d* Dfield = potential_1d_get_dfield_integrate(dfield_background, charge_density, direction_pn, max_dfield_int_error, relerr, abserr, &e);
  if((Dfield == NULL) || (e != NULL)) {
    if(e == NULL) {
      e = g_error_new(POTENTIAL_1D_GERROR_DOMAIN, 3, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_potential: Error getting displacement field in the structure (got NULL, but no error was reported).");
    }
    g_propagate_prefixed_error(err, e, "funkalicious::libpostproc::potential_1d.c::potential_1d_get_potential: Error getting displacement field in the structure:");
    return NULL;
  }
  /*Add the half the difference between the left- and right-hand sides to even out the calculation
   *  in case of residual charge imbalance
   */
#ifdef DEBUG_POTENTIAL_1D_C
  fprintf(stderr, "Dfield difference between positive and negative sides: %g-%g=%g ->", Dfield->storage[Dfield->N-1], Dfield->storage[0], Dfield->storage[Dfield->N-1] + Dfield->storage[0]);
#endif
  density_1d_add_const(Dfield, -(Dfield->storage[Dfield->N-1] + Dfield->storage[0])/2.0);
#ifdef DEBUG_POTENTIAL_1D_C
  fprintf(stderr, " %g-%g=%g\n", Dfield->storage[Dfield->N-1], Dfield->storage[0], Dfield->storage[Dfield->N-1] + Dfield->storage[0]);
#endif
  struct density_1d* phi = potential_1d_get_potential_from_dfield(pot_background, Dfield, epsilon, direction_pn, max_pot_int_err, relerr, abserr, err);
  /*Free the dfield before returning; we don't care about it if the caller used this function*/
  free_density_1d(Dfield);
  return phi;
}

Here is the call graph for this function:

Here is the caller graph for this function:

potential_1d.h File Reference

Enumerations

Functions

Enumeration Type Documentation

Function Documentation