Skip to content

Migrate KiLCA/QL-Balance quadrature to fortnum (3/8) #142

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
marjohma merged 1 commit into from
Jul 1, 2026
Merged

Migrate KiLCA/QL-Balance quadrature to fortnum (3/8) #142

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion KiLCA/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,7 @@ set (kilca_eigparam_main progs/main_eig_param.cpp)
set (libs_root_path ${CMAKE_BINARY_DIR}/external-install)
set (gsl_include_path ${libs_root_path}/gsl/include)

set(EXTERNAL_LIBS slatec lapack blas SUNDIALS::cvode SUNDIALS::nvecserial gsl gslcblas)
set(EXTERNAL_LIBS slatec lapack blas SUNDIALS::cvode SUNDIALS::nvecserial fortnum gsl gslcblas)

# configure a header file to pass some of the CMake settings to the source code
configure_file(
Expand All @@ -122,6 +122,8 @@ configure_file(
include_directories("${CMAKE_CURRENT_BINARY_DIR}")

include_directories ("${gsl_include_path}")
# fortnum C ABI header (fortnum.h); the Fortran fortnum target is linked above.
include_directories ("${fortnum_SOURCE_DIR}/include")
include_directories ("${PROJECT_BINARY_DIR}")
include_directories ("${CMAKE_CURRENT_SOURCE_DIR}/antenna")
include_directories ("${CMAKE_CURRENT_SOURCE_DIR}/background")
Expand Down Expand Up @@ -283,6 +285,7 @@ target_include_directories(kilca_lib PUBLIC
# External dependencies
${sundials_include_path}
${gsl_include_path}
${fortnum_SOURCE_DIR}/include
)

add_executable (kilca_normal_exe ${kilca_normal_main})
Expand Down
20 changes: 3 additions & 17 deletions KiLCA/flre/conductivity/tests/calc_I_array_drift.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,7 @@
#include <complex>
#include <cmath>

#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include "fortnum.h"

#include "constants.h"

Expand Down Expand Up @@ -95,8 +94,6 @@ int calc_imnl_quad_ (int * m_max, int * n_max, int * l_max,
double * x4_re, double * x4_im,
double * Imnl)
{
gsl_set_error_handler_off();

complex<double> x1(*x1_re, *x1_im);
complex<double> x2(*x2_re, *x2_im);
complex<double> x3(*x3_re, *x3_im);
Expand All @@ -117,15 +114,8 @@ x4 = x4 / ep2a;

quad_params qp = {x1, x2, x3, x4, 0, gamma, 0, 0, 0};

size_t limit = 1024;
double epsabs = 1.0e-14, epsrel = 1.0e-14, err;

gsl_integration_workspace * w = gsl_integration_workspace_alloc(limit);

gsl_function F;
F.function = &subint;
F.params = &qp;

double x = 0.0;
double y = 0.0;

Expand Down Expand Up @@ -154,13 +144,11 @@ for(int l = 0; l <= *l_max; ++l)
{
double I_re;
qp.part = 0;
//gsl_integration_qagiu(&F, 0.0e0, epsabs, epsrel, limit, w, &I_re, &err);
gsl_integration_qag(&F, t1, t2, epsabs, epsrel, limit, GSL_INTEG_GAUSS21, w, &I_re, &err);
fortnum_integrate_qag(&subint, t1, t2, epsabs, epsrel, 21, &I_re, &err, &qp);

double I_im;
qp.part = 1;
//gsl_integration_qagiu(&F, 0.0e0, epsabs, epsrel, limit, w, &I_im, &err);
gsl_integration_qag(&F, t1, t2, epsabs, epsrel, limit, GSL_INTEG_GAUSS21, w, &I_im, &err);
fortnum_integrate_qag(&subint, t1, t2, epsabs, epsrel, 21, &I_im, &err, &qp);

InT = (I_re + I_im * I) / pow(ep2a, 0.5e0*(3 + m + n));

Expand All @@ -172,8 +160,6 @@ for(int l = 0; l <= *l_max; ++l)
}
}

gsl_integration_workspace_free(w);

return 0;
}

Expand Down
15 changes: 3 additions & 12 deletions KiLCA/math/fourier/four_transf.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,7 @@
#include "spline.h"
#include "constants.h"

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include "fortnum.h"

/*-----------------------------------------------------------------*/

Expand Down Expand Up @@ -72,25 +70,19 @@ spline_alloc_ (N, 1, dimx, x, C, &sid);

spline_calc_ (sid, yt, 0, 1, NULL, &ierr);

gsl_integration_workspace *W = gsl_integration_workspace_alloc (1000);

struct quad_func_params P = {0, 0, sid, pi/delta, x0};

gsl_function F;
F.function = &func;
F.params = &P;

double re, im, err_re, err_im;

for (k=-M; k<=M; k++)
{
P.k = k;

P.part = 0;
gsl_integration_qag (&F, x0-delta, x0+delta, 1.0e-8, 1.0e-8, 1000, 3, W, &re, &err_re);
fortnum_integrate_qag (&func, x0-delta, x0+delta, 1.0e-8, 1.0e-8, 31, &re, &err_re, &P);

P.part = 1;
gsl_integration_qag (&F, x0-delta, x0+delta, 1.0e-8, 1.0e-8, 1000, 3, W, &im, &err_im);
fortnum_integrate_qag (&func, x0-delta, x0+delta, 1.0e-8, 1.0e-8, 31, &im, &err_im, &P);

//fprintf (stdout, "\nk = %d:\tresult = %lg%+lgi\terror = %lg%+lgi", k, re, im, err_re, err_im);

Expand All @@ -102,7 +94,6 @@ delete [] C;
delete [] yt;

spline_free_ (sid);
gsl_integration_workspace_free (W);

return 0;
}
Expand Down
4 changes: 4 additions & 0 deletions QL-Balance/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,9 @@ set_source_files_properties(${base_balance_src}/kim_wave_code_adapter.f90
# Include central slatec math sources via cache variable defined by common/math
target_include_directories(ql-balance_lib PUBLIC ${SLATEC_INCLUDE_DIR})

# fortnum C ABI header (fortnum.h) for vel_integral.cpp.
target_include_directories(ql-balance_lib PUBLIC ${fortnum_SOURCE_DIR}/include)

find_package(HDF5 REQUIRED COMPONENTS C HL Fortran Fortran_HL)
message(" ====== HDF5 was found! ======")

Expand All @@ -150,6 +153,7 @@ target_link_libraries(ql-balance_lib PUBLIC
SUNDIALS::cvode
SUNDIALS::nvecserial
blas
fortnum
gsl
gslcblas
hdf5::hdf5
Expand Down
18 changes: 5 additions & 13 deletions QL-Balance/src/base/vel_integral.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,7 @@

#include "constants.h"

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include "fortnum.h"

/*-----------------------------------------------------------------*/

Expand Down Expand Up @@ -65,14 +63,8 @@ return exp(-x*x)/(pow((P->a)*x + P->omE, 2)+(P->nu)*(P->nu))*real(field_fac)*ind

void calc_velocity_integral_ (int ind, double vT, double ks, double kp, double omE, double nu, double *Es, double *Ep, double *res)
{
gsl_integration_workspace *W = gsl_integration_workspace_alloc (10000);

struct quad_func_params P;

gsl_function F;
F.function = &func;
F.params = &P;

//set structure:
P.ind = ind;
P.a = sqrt(2.0)*vT*kp;
Expand Down Expand Up @@ -100,10 +92,10 @@ for (i=0; i<10000; i++)

double err;

gsl_integration_qagi (&F, 1.0e-6, 1.0e-6, 10000, W, res, &err);
//gsl_integration_qag (&F, -10.0, 10.0, 1.0e-6, 1.0e-6, 10000, 3, W, res, &err);

gsl_integration_workspace_free (W);
// gsl_integration_qagi over (-inf, +inf): fortnum QAGIU with inf=2 splits at
// the bound and sums the two semi-infinite transforms (bound ignored except
// as the split point, matching QAGI's split at 0).
fortnum_integrate_qagiu (&func, 0.0, 2, 1.0e-6, 1.0e-6, res, &err, &P);
}

/*-----------------------------------------------------------------*/
Loading