potential_1d.c

Go to the documentation of this file.

/*
 *Gets the potential from a 1D charge density through a number of methods.

 Copyright (C) 2010 Joseph Pingenot

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU Affero General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Affero General Public License for more details.

 You should have received a copy of the GNU Affero General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <math.h>
#include <libpostproc/potential_1d.h>
#include <libpostproc/unit_conversion.h>
#ifdef OPENMP
#include <omp.h>
#endif

      
//#define DEBUG_POTENTIAL_1D_C

/*Argument to the GSL integration workspace allocation calls*/
#define MAX_GSL_INTEGRATION_INTERVALS 4096

#define POTENTIAL_1D_GERROR_DOMAIN 73

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) {
  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;
}

struct potential_integration_method_integrand_params {
  gsl_spline* dfield;
  gsl_interp_accel* dfield_accel;
  gsl_spline* epsilon;
  gsl_interp_accel* epsilon_accel;
#ifdef DEBUG_POTENTIAL_1D_C
  long int call;
  FILE* file;
#endif
};

static double potential_integration_method_integrand(double x, void* privdat) {
  struct potential_integration_method_integrand_params* p = (struct potential_integration_method_integrand_params*)privdat;
  /*Need this minus sign; D = \epsilon -\nabla V*/
  double value=-gsl_spline_eval(p->dfield, x, p->dfield_accel)/gsl_spline_eval(p->epsilon, x, p->epsilon_accel);
#ifdef DEBUG_POTENTIAL_1D_C
  fprintf(p->file, "%ld\t%g\t%g\n", p->call++, x, value);
#endif
  return value;
}

struct dfield_integrand_params {
  gsl_spline* charge;
  gsl_interp_accel* charge_accel;
#ifdef DEBUG_POTENTIAL_1D_C
  long int call;
  FILE* file;
#endif
};
static double dfield_integrand(double x, void* privdat) {
  struct dfield_integrand_params* p = (struct dfield_integrand_params*)privdat;
  double value= gsl_spline_eval(p->charge, x, p->charge_accel)*oo_epsilon0;
#ifdef DEBUG_POTENTIAL_1D_C
  fprintf(p->file, "%ld\t%g\t%g\n", p->call++, x, value);
#endif
  return value;
}


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* 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;
}

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) {
  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;
}

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) {
  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;
}