From 57157247d4bcb65ab21d1a30cd9d9b26203ddb73 Mon Sep 17 00:00:00 2001 From: JanThorbecke Date: Thu, 18 Jul 2024 11:10:32 +0200 Subject: [PATCH] Adding C-kernel --- .gitignore | 4 + Makefile | 3 + empymod/Makefile | 32 ++++ empymod/_ckernel.py | 133 +++++++++++++ empymod/fields.c | 284 +++++++++++++++++++++++++++ empymod/greenfct.c | 414 ++++++++++++++++++++++++++++++++++++++++ empymod/reflections.c | 144 ++++++++++++++ empymod/transform.py | 33 +++- empymod/utils.py | 4 +- empymod/wavenumber.c | 125 ++++++++++++ tests/test_kernel.py | 26 ++- tests/test_model.py | 4 + tests/test_transform.py | 2 +- 13 files changed, 1190 insertions(+), 18 deletions(-) create mode 100644 empymod/Makefile create mode 100644 empymod/_ckernel.py create mode 100644 empymod/fields.c create mode 100644 empymod/greenfct.c create mode 100644 empymod/reflections.c create mode 100644 empymod/wavenumber.c diff --git a/.gitignore b/.gitignore index 7091d0a..fa28ea4 100644 --- a/.gitignore +++ b/.gitignore @@ -26,3 +26,7 @@ empymod/version.py build/ dist/ empymod.egg-info/ + +# C +*.o +*.so diff --git a/Makefile b/Makefile index 599cedc..7ca7163 100644 --- a/Makefile +++ b/Makefile @@ -14,9 +14,11 @@ help: @echo "" install: + cd empymod && make && cd .. python -m pip install -e . dev-install: + cd empymod && make && cd .. python -m pip install -e .[all] .ONESHELL: @@ -46,6 +48,7 @@ linkcheck: cd docs && make linkcheck clean: + cd empymod && make deepclean && cd .. python -m pip uninstall empymod -y rm -rf build/ dist/ .eggs/ empymod.egg-info/ empymod/version.py # build rm -rf */__pycache__/ */*/__pycache__/ # python cache diff --git a/empymod/Makefile b/empymod/Makefile new file mode 100644 index 0000000..f802170 --- /dev/null +++ b/empymod/Makefile @@ -0,0 +1,32 @@ +#CC= /opt/local/bin/gcc-mp-13 + +all: fields.so reflections.so greenfct.so wavenumber.so + +fields.so: fields.o + gcc -shared -o $@ fields.o -lm + +reflections.so: reflections.o + gcc -shared -o $@ reflections.o -lm + +greenfct.so: greenfct.o reflections.o fields.o + gcc -shared -o $@ $^ -lm + +wavenumber.so: wavenumber.o greenfct.o reflections.o fields.o + gcc -shared -o $@ $^ -lm + +#CFLAGS= -fPIC -Ofast -ftree-loop-vectorize -fopt-info-optimized -fopt-info-vec-missed -march=znver2 +#-fopt-info-optimized -fopenmp +CFLAGS= -fPIC -Ofast -ftree-loop-vectorize -ffast-math -fopt-info-optimized -funroll-all-loops +CFLAGS= -fPIC -O3 -ftree-loop-vectorize -funroll-all-loops +#CFLAGS= -fPIC -O0 -g + +.SUFFIXES : .o .c +.c.o : + $(CC) -c $(CFLAGS) $(OPTC) $< + +clean: + rm -f *.o + +deepclean: + rm -rf *.o *.so __pycache__ + diff --git a/empymod/_ckernel.py b/empymod/_ckernel.py new file mode 100644 index 0000000..478dc4e --- /dev/null +++ b/empymod/_ckernel.py @@ -0,0 +1,133 @@ +import os +import ctypes as ct + +import numpy as np + + +C_DOUBLEP = ct.POINTER(ct.c_double) +path = os.path.dirname(os.path.realpath(__file__)) +cwavenumber = np.ctypeslib.load_library("wavenumber", path) +cgreenfct = np.ctypeslib.load_library("greenfct", path) +creflections = np.ctypeslib.load_library("reflections", path) +cfields = np.ctypeslib.load_library("fields", path) + + +def wavenumber(zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH, zetaV, lambd, + ab, xdirect, msrc, mrec): + """C-Wrapper for wavenumber.""" + nfreq, nlayer = etaH.shape + noff, nlambda = lambd.shape + PJ0 = np.zeros((nfreq, noff, nlambda), dtype=complex) + PJ1 = np.zeros((nfreq, noff, nlambda), dtype=complex) + PJ0b = np.zeros((nfreq, noff, nlambda), dtype=complex) + cwavenumber.wavenumber( + int(nfreq), + int(noff), + int(nlayer), + int(nlambda), + ct.c_double(zsrc), + ct.c_double(zrec), + int(lsrc), + int(lrec), + depth.ctypes.data_as(C_DOUBLEP), + etaH.ctypes.data_as(C_DOUBLEP), + etaV.ctypes.data_as(C_DOUBLEP), + zetaH.ctypes.data_as(C_DOUBLEP), + zetaV.ctypes.data_as(C_DOUBLEP), + lambd.ctypes.data_as(C_DOUBLEP), + int(ab), + int(xdirect), + int(msrc), + int(mrec), + PJ0.ctypes.data_as(C_DOUBLEP), + PJ1.ctypes.data_as(C_DOUBLEP), + PJ0b.ctypes.data_as(C_DOUBLEP), + ) + if ab not in [11, 22, 24, 15, 33]: + PJ0 = None + if ab not in [11, 12, 21, 22, 14, 24, 15, 25]: + PJ0b = None + if ab in [33, ]: + PJ1 = None + return PJ0, PJ1, PJ0b + + +def greenfct(zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH, zetaV, lambd, + ab, xdirect, msrc, mrec): + """C-Wrapper for greenfct.""" + nfreq, nlayer = etaH.shape + noff, nlambda = lambd.shape + GTM = np.zeros((nfreq, noff, nlambda), dtype=complex) + GTE = np.zeros((nfreq, noff, nlambda), dtype=complex) + cgreenfct.greenfct( + int(nfreq), + int(noff), + int(nlayer), + int(nlambda), + ct.c_double(zsrc), + ct.c_double(zrec), + int(lsrc), + int(lrec), + depth.ctypes.data_as(C_DOUBLEP), + etaH.ctypes.data_as(C_DOUBLEP), + etaV.ctypes.data_as(C_DOUBLEP), + zetaH.ctypes.data_as(C_DOUBLEP), + zetaV.ctypes.data_as(C_DOUBLEP), + lambd.ctypes.data_as(C_DOUBLEP), + int(ab), + int(xdirect), + int(msrc), + int(mrec), + GTM.ctypes.data_as(C_DOUBLEP), + GTE.ctypes.data_as(C_DOUBLEP), + ) + return GTM, GTE + + +def reflections(depth, e_zH, Gam, lrec, lsrc): + """C-Wrapper for reflections.""" + nfreq, noff, nlayer, nlambda = Gam.shape + maxl = max([lrec, lsrc]) + minl = min([lrec, lsrc]) + nl = maxl-minl+1 + Rp = np.zeros((nfreq, noff, nl, nlambda), dtype=complex) + Rm = np.zeros((nfreq, noff, nl, nlambda), dtype=complex) + creflections.reflections( + int(nfreq), + int(noff), + int(nlayer), + int(nlambda), + depth.ctypes.data_as(C_DOUBLEP), + e_zH.ctypes.data_as(C_DOUBLEP), + Gam.ctypes.data_as(C_DOUBLEP), + int(lrec), + int(lsrc), + Rp.ctypes.data_as(C_DOUBLEP), + Rm.ctypes.data_as(C_DOUBLEP), + ) + return Rp, Rm + + +def fields(depth, Rp, Rm, Gam, lrec, lsrc, zsrc, ab, TM): + """C-Wrapper for fields.""" + nfreq, noff, nlayer, nlambda = Rp.shape + Pu = np.zeros((nfreq, noff, nlambda), dtype=complex) + Pd = np.zeros((nfreq, noff, nlambda), dtype=complex) + cfields.fields( + int(nfreq), + int(noff), + int(nlayer), + int(nlambda), + depth.ctypes.data_as(C_DOUBLEP), + Rp.ctypes.data_as(C_DOUBLEP), + Rm.ctypes.data_as(C_DOUBLEP), + Gam.ctypes.data_as(C_DOUBLEP), + int(lrec), + int(lsrc), + ct.c_double(zsrc), + int(ab), + int(TM), + Pu.ctypes.data_as(C_DOUBLEP), + Pd.ctypes.data_as(C_DOUBLEP), + ) + return Pu, Pd diff --git a/empymod/fields.c b/empymod/fields.c new file mode 100644 index 0000000..206bb95 --- /dev/null +++ b/empymod/fields.c @@ -0,0 +1,284 @@ +#include +#include +#include +#include +#include +#include +#include + +bool isina(int ab, const int *a, int n) { + int i; + bool ok; + + ok=false; + for (i=0;i lsrc)) { + memset(Pu,0,nfreq*noff*nlambda*sizeof(double complex)); + continue; + } + // No downgoing field if rec is in first layer or above src + if (up==0 && (lrec==0 || lrec < lsrc)) { + memset(Pd,0,nfreq*noff*nlambda*sizeof(double complex)); + continue; + } + + // Swaps if up=True + if (up==1){ + if (!last_layer) { + ftmp=dp; + dp=dm; + dm=ftmp; + } + else{ + dp = dm; + } + // reference + Rmp = Rp; + Rpm = Rm; + itmp=first_layer; + first_layer=last_layer; + last_layer=itmp; + rsrcl = nlsr-1; // src-layer in refl. (Rp/Rm), last (nlsr-1) if up + izstart=0; + izend=nlsr-2; + isr = lrec; + last = 0; + pup = 1; + if (!plus) mupm = -1; + P = Pu; + } + else{ + // reference + Rmp = Rm; + Rpm = Rp; + P = Pd; + } + + // Calculate Pu+, Pu-, Pd+, Pd- + if (lsrc == lrec) { // rec in src layer; Eqs 81/82, A-8/A-9 + if (last_layer) { // If src/rec are in top (up) or bottom (down) layer + for (i=0;i 2){ + for (iz=izstart;iz +#include +#include +#include +#include +#include +#include + +#define MAX(x,y) ((x) > (y) ? (x) : (y)) +#define MIN(x,y) ((x) < (y) ? (x) : (y)) + +bool isina(int ab, const int *a, int n); + +// Define the reflections function +void reflections(int nfreq, int noff, int nlayer, int nlambda, double *depth, double complex *e_zH, double complex *Gam, int lrec, int lsrc, double complex *Rp, double complex *Rm) ; + +// Define the fields function +void fields(int nfreq, int noff, int nlayer, int nlambda, double *depth, double complex *Rp, double complex *Rm, double complex *Gam, int lrec, int lsrc, double zsrc, int ab, bool TM, double complex *Pu, double complex *Pd); + +// Define the greenfct function +// Calculate Green's function for TM and TE. +void greenfct(int nfreq, int noff, int nlayer, int nlambda, double zsrc, double zrec, int lsrc, int lrec, double* depth, double complex* etaH, double complex* etaV, double complex* zetaH, double complex* zetaV, double* lambd, int ab, int xdirect, int msrc, int mrec, double complex* GTM, double complex* GTE) +{ + int g1, g2, g3, r1, r2, n1, n2, n3, nlay; + int i, ii, iv, iz, dsign, minl, maxl; + double complex *Gam, *Rp, *Rm, *gamTM, *gamTE; + double complex *Wu, *Wd, *Pu, *Pd, *green; + double complex h_div_v, h_times_h, fexp, fact; + double complex *e_zH; + double complex *e_zV; + double complex *z_eH; + double complex *letaH, *letaV, *lzetaH, *lzetaV, *ptmp, directf; + double temp, ddepth, dfsign, pmw, l2; + const int minus_ab[] = {11, 12, 13, 14, 15, 21, 22, 23, 24, 25}; + + // GTM/GTE have shape (frequency, offset, lambda). + // gamTM/gamTE have shape (frequency, offset, layer, lambda): + + g1 = nlambda; + g2 = g1*nlayer; + g3 = g2*noff; + + r1=nlambda; + r2=noff*r1; + + maxl = MAX(lrec, lsrc); + minl = MIN(lrec, lsrc); + nlay = (maxl-minl+1); + n1 = nlambda; + n2 = n1*nlay; + n3 = n2*noff; + + //fprintf(stderr,"nfreq=%d noff=%d nlayer=%d nlambda=%d lrec=%d lsrc=%d mrec=%d msrc=%d nlay=%d\n",nfreq, noff, nlayer, nlambda, lrec, lsrc, mrec, msrc, nlay); + //fprintf(stderr,"r2=%d r1=%d g3=%d g2=%d g1=%d\n",r2, r1, g3, g2, g1); + + gamTM = (double complex *) calloc(nfreq * noff * nlayer * nlambda , sizeof(double complex)); + gamTE = (double complex *) calloc(nfreq * noff * nlayer * nlambda , sizeof(double complex)); + Rp = (double complex *) calloc(nfreq * noff * nlay * nlambda , sizeof(double complex)); + Rm = (double complex *) calloc(nfreq * noff * nlay * nlambda , sizeof(double complex)); + Wu = (double complex *) malloc(nfreq * noff * nlambda * sizeof(double complex)); + Wd = (double complex *) malloc(nfreq * noff * nlambda * sizeof(double complex)); + Pu = (double complex *) malloc(nfreq * noff * nlambda * sizeof(double complex)); + Pd = (double complex *) malloc(nfreq * noff * nlambda * sizeof(double complex)); + + //green = (double complex *) malloc(nfreq * noff * nlambda * sizeof(double complex)); + //fexp = (double complex *) malloc(nfreq * noff * nlambda * sizeof(double complex)); + letaH = (double complex *) malloc(nfreq * nlayer * sizeof(double complex)); + letaV = (double complex *) malloc(nfreq * nlayer * sizeof(double complex)); + lzetaH = (double complex *) malloc(nfreq * nlayer * sizeof(double complex)); + lzetaV = (double complex *) malloc(nfreq * nlayer * sizeof(double complex)); + +// TODO use pointers to assign GTM to green in TM loop +// TODO use pointers to assign Gam to GamTM? + + for (int i = 0; i < nfreq; i++) { + for (int j = 0; j < nlayer; j++) { //TODO can this be done with pointers and use the - in equations? + letaH[i*nlayer+j] = etaH[i*nlayer+j]; + lzetaH[i*nlayer+j] = zetaH[i*nlayer+j]; + letaV[i*nlayer+j] = etaV[i*nlayer+j]; + lzetaV[i*nlayer+j] = zetaV[i*nlayer+j]; + } + } + // Reciprocity switches for magnetic receivers + if (mrec) { + if (msrc) { //If src is also magnetic, switch eta and zeta (MM => EE). + // G^mm_ab(s, r, e, z) = -G^ee_ab(s, r, -z, -e) + for (int i = 0; i < nfreq; i++) { + for (int j = 0; j < nlayer; j++) { //TODO can this be done with pointers and use the - in equations? + //fprintf(stderr," i=%d j=%d etaH=%e %e\n", i, j, crealf(etaH[i*nlayer+j]), cimagf(etaH[i*nlayer+j])); + //fprintf(stderr," i=%d j=%d etaV=%e %e\n", i, j, crealf(etaV[i*nlayer+j]), cimagf(etaV[i*nlayer+j])); + letaH[i*nlayer+j] = -zetaH[i*nlayer+j]; + lzetaH[i*nlayer+j] = -etaH[i*nlayer+j]; + letaV[i*nlayer+j] = -zetaV[i*nlayer+j]; + lzetaV[i*nlayer+j] = -etaV[i*nlayer+j]; + } + } + } else { //If src is electric, swap src and rec (ME => EM). + // G^me_ab(s, r, e, z) = -G^em_ba(r, s, e, z) + temp = zsrc; + zsrc = zrec; + zrec = temp; + temp = lsrc; + lsrc = lrec; + lrec = temp; + } + } + + for (int TM = 0; TM < 2; TM++) { + // Continue if Green's function not required + if (TM && (ab == 16 || ab == 26)) { + continue; + } else if (!TM && (ab == 13 || ab == 23 || ab == 31 || ab == 32 || ab == 33 || ab == 34 || ab == 35)) { + continue; + } + + // Define eta/zeta depending if TM or TE + if (TM) { + e_zH = letaH; + e_zV = letaV; + z_eH = lzetaH; + green = GTM; + Gam = gamTM; +//TODO set pointer for green and Gam + } else { + e_zH = lzetaH; + e_zV = lzetaV; + z_eH = letaH; + green = GTE; + Gam = gamTE; + } + + // Uppercase gamma + for (int i = 0; i < nfreq; i++) { + for (int iii = 0; iii < nlayer; iii++) { + h_div_v = e_zH[i*nlayer+iii] / e_zV[i*nlayer+iii]; + h_times_h = z_eH[i*nlayer+iii] * e_zH[i*nlayer+iii]; + //fprintf(stderr," i=%d iii=%d h_div_v=%e %e\n", i, iii, crealf(h_div_v), cimagf(h_div_v)); + //fprintf(stderr," i=%d iii=%d h_times_h=%e %e\n", i, iii, crealf(h_times_h), cimagf(h_times_h)); + for (int ii = 0; ii < noff; ii++) { + for (int iv = 0; iv < nlambda; iv++) { + l2 = lambd[ii*nlambda+iv] * lambd[ii*nlambda+iv]; + //fprintf(stderr," ii=%d iv=%d l2=%e %e\n", ii, iv, crealf(l2), cimagf(l2)); + Gam[i*g3+ii*g2+iii*g1+iv] = csqrt(h_div_v * l2 + h_times_h); + //fprintf(stderr,"Gam i=%d ii=%d iii=%d iv=%d %e %e\n", i, ii, iii, iv, crealf(Gam[i*g3+ii*g2+iii*g1+iv]), cimagf(Gam[i*g3+ii*g2+iii*g1+iv])); + } + } + } + } + + // Gamma in receiver layer + + // Reflection (coming from below (Rp) and above (Rm) rec) + if (nlayer > 1) { + // TODO set Rp Rm to zero if needed + reflections(nfreq, noff, nlayer, nlambda, depth, e_zH, Gam, lrec, lsrc, Rp, Rm); + // Field propagators + // (Up- (Wu) and downgoing (Wd), in rec layer); Eq 74 + + if (lrec != nlayer - 1) { + ddepth = depth[lrec + 1] - zrec; + for (int i = 0; i < nfreq; i++) { + for (int ii = 0; ii < noff; ii++) { + for (int iv = 0; iv < nlambda; iv++) { + Wu[i*r2+ii*r1+iv] = cexp(-Gam[i*g3+ii*g2+lrec*g1+iv] * ddepth); + } + } + } + } + + if (lrec != 0) { + ddepth = zrec - depth[lrec]; + for (int i = 0; i < nfreq; i++) { + for (int ii = 0; ii < noff; ii++) { + for (int iv = 0; iv < nlambda; iv++) { + Wd[i*r2+ii*r1+iv] = cexp(-Gam[i*g3+ii*g2+lrec*g1+iv] * ddepth); + } + } + } + } + + // Field at rec level (coming from below (Pu) and above (Pd) rec) + // set Pu Pd to zero + //memset(Pu,0,nfreq*noff*nlambda*sizeof(double complex)); + //memset(Pd,0,nfreq*noff*nlambda*sizeof(double complex)); + fields(nfreq, noff, nlayer, nlambda, depth, Rp, Rm, Gam, lrec, lsrc, zsrc, ab, TM, Pu, Pd); + } + + // Green's functions + + if (lsrc == lrec) { //Rec in src layer; Eqs 108, 109, 110, 117, 118, 122 + // Green's function depending on + // (If only one layer, no reflections/fields) + if (nlayer > 1 && (ab == 13 || ab == 23 || ab == 31 || ab == 32 || ab == 14 || ab == 24 || ab == 15 || ab == 25)) { + for (int i = 0; i < nfreq; i++) { + for (int ii = 0; ii < noff; ii++) { + for (int iv = 0; iv < nlambda; iv++) { + green[i*r2+ii*r1+iv] = Pu[i*r2+ii*r1+iv]*Wu[i*r2+ii*r1+iv]; + green[i*r2+ii*r1+iv] -= Pd[i*r2+ii*r1+iv]*Wd[i*r2+ii*r1+iv]; + } + } + } + } else if (nlayer > 1) { + for (int i = 0; i < nfreq; i++) { + for (int ii = 0; ii < noff; ii++) { + for (int iv = 0; iv < nlambda; iv++) { + green[i*r2+ii*r1+iv] = Pu[i*r2+ii*r1+iv]*Wu[i*r2+ii*r1+iv]; + green[i*r2+ii*r1+iv] += Pd[i*r2+ii*r1+iv]*Wd[i*r2+ii*r1+iv]; + } + } + } + } + + // Direct field, if it is computed in the wavenumber domain + if (!xdirect) { + ddepth = abs(zsrc - zrec); + if ((zrec - zsrc) < 0 ) dsign = -1; + else dsign = 1; + + // Swap TM for certain + dfsign = 1; + if (TM && isina(ab, minus_ab, 10)) { + dfsign = -1; + } + //# Multiply by zrec-zsrc-sign for certain + if ((ab == 11 || ab == 12 || ab == 13 || ab == 14 || ab == 15 || ab == 21 || ab == 22 || ab == 23 || ab == 24 || ab == 25)) { + dfsign *= dsign; + } + + for (int i = 0; i < nfreq; i++) { + for (int ii = 0; ii < noff; ii++) { + for (int iv = 0; iv < nlambda; iv++) { + // Direct field + directf = dfsign*cexp(-Gam[i*g3+ii*g2+lrec*g1+iv]*ddepth); + + // Add direct field to Green's function + green[i*r2+ii*r1+iv] += directf; + } + } + } + + // Implementation + } + } else { + // Calculate exponential factor + if (lrec == nlayer-1) { + ddepth = 0; + } + else { + ddepth = depth[lrec+1] - depth[lrec]; + } + +/* replaced by scalar in loops + for (int i = 0; i < nfreq; i++) { + for (int ii = 0; ii < noff; ii++) { + for (int iv = 0; iv < nlambda; iv++) { + fexp[i*r2+ii*r1+iv] = cexp(-Gam[i*g3+ii*g2+lrec*g1+iv]*ddepth); + } + } + } +*/ + + // Sign-switch for Green calculation + if (TM && isina(ab, minus_ab, 10)) { + pmw = -1; + } + else { + pmw = 1; + } + + if (lrec < lsrc) { // Rec above src layer: Pd not used + // Eqs 89-94, A18-A23, B13-B15 + for (int i = 0; i < nfreq; i++) { + for (int ii = 0; ii < noff; ii++) { + for (int iv = 0; iv < nlambda; iv++) { + fexp = cexp(-Gam[i*g3+ii*g2+lrec*g1+iv]*ddepth); + green[i*r2+ii*r1+iv] = Pu[i*r2+ii*r1+iv]*( + Wu[i*r2+ii*r1+iv] + pmw*Rm[i*n3+ii*n2+0*n1+iv] * + fexp*Wd[i*r2+ii*r1+iv]); + } + } + } + } + else if (lrec > lsrc) { // rec below src layer: Pu not used + // Eqs 97-102 A26-A30, B16-B18 + for (int i = 0; i < nfreq; i++) { + for (int ii = 0; ii < noff; ii++) { + for (int iv = 0; iv < nlambda; iv++) { + fexp = cexp(-Gam[i*g3+ii*g2+lrec*g1+iv]*ddepth); + green[i*r2+ii*r1+iv] = Pd[i*r2+ii*r1+iv]*( + pmw*Wd[i*r2+ii*r1+iv] + + Rp[i*n3+ii*n2+abs(lsrc-lrec)*n1+iv] * + fexp*Wu[i*r2+ii*r1+iv]); + //fprintf(stderr,"Rp*fexpWu i=%d ii=%d iv=%d %e %e\n", i, ii, iv, crealf(Rp[i*g3+ii*g2+abs(lsrc-lrec)*g1+iv]*fexp*Wu[i*r2+ii*r1+iv]), cimagf(Rp[i*g3+ii*g2+abs(lsrc-lrec)*g1+iv]*fexp*Wu[i*r2+ii*r1+iv])); + //fprintf(stderr,"Rp i=%d ii=%d iv=%d %e %e\n", i, ii, iv, crealf(Rp[i*n3+ii*n2+abs(lsrc-lrec)*n1+iv]), cimagf(Rp[i*n3+ii*n2+abs(lsrc-lrec)*n1+iv])); + } + } + } + } + } + // Store in corresponding variable + // TODO: => done to check + //if (TM) { + //gamTM = Gam; + //memcpy(GTM, green, nfreq*noff*nlambda*sizeof(double complex)); + //gamTM, GTM = Gam, green + //} + //else{ + //gamTE = Gam; + //memcpy(GTE, green, nfreq*noff*nlambda*sizeof(double complex)); + //gamTE, GTE = Gam, green + //} + + } // end of TM loop + + // ** AB-SPECIFIC FACTORS AND CALCULATION OF PTOT'S + // These are the factors inside the integrals + // Eqs 105-107, 111-116, 119-121, 123-128 + // + + if (ab == 11 || ab == 12 || ab == 21 || ab == 22) { + for (int i = 0; i < nfreq; ++i) { + fact = 1.0 / letaH[i*nlayer+lrec]; + for (int ii = 0; ii < noff; ++ii) { + for (int iv = 0; iv < nlambda; ++iv) { + //fprintf(stderr,"GTM i=%d ii=%d iv=%d %e %e\n", i, ii, iv, crealf(GTM[i*r2+ii*r1+iv]), cimagf(GTM[i*r2+ii*r1+iv])); + GTM[i*r2+ii*r1+iv] *= fact * gamTM[i*g3+ii*g2+lrec*g1+iv]; + GTE[i*r2+ii*r1+iv] *= lzetaH[i*nlayer+lsrc] / gamTE[i*g3+ii*g2+lsrc*g1+iv]; + //fprintf(stderr,"gamTE i=%d ii=%d iv=%d %e %e\n", i, ii, iv, crealf(gamTE[i*g3+ii*g2+lsrc*g1+iv]), cimagf(gamTE[i*g3+ii*g2+lsrc*g1+iv])); + } + } + } + } else if (ab == 14 || ab == 15 || ab == 24 || ab == 25) { + for (int i = 0; i < nfreq; ++i) { + fact = letaH[i*nlayer+lsrc] / letaH[i*nlayer+lrec]; + for (int ii = 0; ii < noff; ++ii) { + for (int iv = 0; iv < nlambda; ++iv) { + GTM[i*r2+ii*r1+iv] *= fact * gamTM[i*g3+ii*g2+lrec*g1+iv]; + GTM[i*r2+ii*r1+iv] /= gamTM[i*g3+ii*g2+lsrc*g1+iv]; + } + } + } + } else if (ab == 13 || ab == 23) { + memset(GTE,0,nfreq*noff*nlambda*sizeof(double complex)); + for (int i = 0; i < nfreq; ++i) { + fact = letaH[i*nlayer+lsrc] / letaH[i*nlayer+lrec] / letaV[i*nlayer+lsrc]; + for (int ii = 0; ii < noff; ++ii) { + for (int iv = 0; iv < nlambda; ++iv) { + GTM[i*r2+ii*r1+iv] *= -fact * gamTM[i*g3+ii*g2+lrec*g1+iv]; + GTM[i*r2+ii*r1+iv] /= gamTM[i*g3+ii*g2+lsrc*g1+iv]; + } + } + } + } else if (ab == 31 || ab == 32) { + memset(GTE,0,nfreq*noff*nlambda*sizeof(double complex)); + for (int i = 0; i < nfreq; ++i) { + fact= 1.0/letaV[i*nlayer+lrec]; + for (int ii = 0; ii < noff; ++ii) { + for (int iv = 0; iv < nlambda; ++iv) { + GTM[i*r2+ii*r1+iv] *= fact; + } + } + } + } else if (ab == 34 || ab == 35) { + memset(GTE,0,nfreq*noff*nlambda*sizeof(double complex)); + for (int i = 0; i < nfreq; ++i) { + fact = letaH[i*nlayer+lsrc]/letaV[i*nlayer+lrec]; + for (int ii = 0; ii < noff; ++ii) { + for (int iv = 0; iv < nlambda; ++iv) { + GTM[i*r2+ii*r1+iv] *= fact/gamTM[i*g3+ii*g2+lsrc*g1+iv]; + } + } + } + } else if (ab == 16 || ab == 26) { + memset(GTM,0,nfreq*noff*nlambda*sizeof(double complex)); + for (int i = 0; i < nfreq; ++i) { + fact = lzetaH[i*nlayer+lsrc]/lzetaV[i*nlayer+lsrc]; + for (int ii = 0; ii < noff; ++ii) { + for (int iv = 0; iv < nlambda; ++iv) { + GTE[i*r2+ii*r1+iv] *= fact/gamTE[i*g3+ii*g2+lsrc*g1+iv]; + } + } + } + + } else if (ab == 33) { + memset(GTE,0,nfreq*noff*nlambda*sizeof(double complex)); + for (int i = 0; i < nfreq; ++i) { + fact = letaH[i*nlayer+lsrc]/letaV[i*nlayer+lsrc]/letaV[i*nlayer+lrec]; + for (int ii = 0; ii < noff; ++ii) { + for (int iv = 0; iv < nlambda; ++iv) { + GTM[i*r2+ii*r1+iv] *= fact/gamTM[i*g3+ii*g2+lsrc*g1+iv]; + } + } + } + } + + free(letaH); + free(letaV); + free(lzetaH); + free(lzetaV); + free(Rp); + free(Rm); + free(Wu); + free(Wd); + free(Pu); + free(Pd); + free(gamTM); + free(gamTE); + + // Return Green's functions GTM and GTE + return; + +} + + diff --git a/empymod/reflections.c b/empymod/reflections.c new file mode 100644 index 0000000..23517f7 --- /dev/null +++ b/empymod/reflections.c @@ -0,0 +1,144 @@ +#include +#include +#include +#include +#include +#include +#include + +#define MAX(x,y) ((x) > (y) ? (x) : (y)) +#define MIN(x,y) ((x) < (y) ? (x) : (y)) + +void reflections(int nfreq, int noff, int nlayer, int nlambda, double *depth, double complex *e_zH, double complex *Gam, int lrec, int lsrc, double complex *Rp, double complex *Rm) +{ + int maxl, minl, pm, izout, minmax, izout0; + int *layer_count, lcount; + bool shiftplus, shiftminus; + int g1, g2, g3, n1, n2, n3, nlay; + int i, ii, iv, iz; + double complex rloc, *tRef, *out, *Ref; + double complex ra, rb, rloca, rlocb, term; + double ddepth; + + g1 = nlambda; + g2 = g1*nlayer; + g3 = g2*noff; + + // Get numbers and max/min layer. + maxl = MAX(lrec, lsrc); + minl = MIN(lrec, lsrc); + nlay = (maxl-minl+1); + n1 = nlambda; + n2 = n1*nlay; + n3 = n2*noff; + + //fprintf(stderr,"nfreq=%d noff=%d nlayer=%d nlambda=%d lrec=%d lsrc=%d nlay=%d\n",nfreq, noff, nlayer, nlambda, lrec, lsrc,nlay); + //fprintf(stderr,"r2=%d r1=%d n3=%d n2=%d n1=%d\n",r2, r1, n3, n2, n1); + // Pre-allocate tRef + // TODO eliminate tRef array + tRef = (double complex *)malloc(nlambda*sizeof(double complex)); + layer_count = (int *)malloc(sizeof(int) * nlayer); + + // Loop over Rp, Rm + for (int plus = 0; plus <= 1; plus++) { + + // Switches depending if plus or minus + if (plus==1) { + pm = 1; + //layer_count = np.arange(depth.size-2, minl-1, -1) + for (int i = nlayer - 2; i > minl-1; i--) { + layer_count[nlayer-2 - i] = i; +// fprintf(stderr,"plus=%d i=%d layer_count[%d] = %d \n",plus, i, nlayer - 2 - i, layer_count[nlayer - 2 - i]); + } + lcount = nlayer-2-(minl-1); + izout = abs(lsrc-lrec); + minmax = pm*maxl; + out = Rp; + } + else { + pm = -1; + //layer_count = np.arange(1, maxl+1, 1) + for (int i = 0; i < maxl; i++) { + layer_count[i] = i + 1; +// fprintf(stderr,"plus=%d i=%d layer_count[%d] = %d \n",plus, i, i, layer_count[i]); + } + lcount = maxl; + izout = 0; + minmax = pm*minl; + out = Rm; + } + + // If rec in last and rec below src (plus) or + // if rec in first and rec above src (minus), shift izout + if ( (lreclsrc) && (lrec==nlayer-1) && (plus==1)) {shiftminus = true;} + else {shiftminus = false;} + if (shiftplus || shiftminus) {izout -= pm;} + + izout0=izout; + + // Calculate the reflection + // Eqs 65, A-12 + for (i=0;i 0 && il == lcount-1 ) { + for (int iv = 0; iv < nlambda; iv++) { + out[i*n3 + ii*n2 + iv] = tRef[iv]; + } + } + + if ((lrec != lsrc) && (pm*iz <= minmax)) { + izout -= pm; + } + + } // end of lcount layer loop + } // end of noff loop + } // end of nfreq loop + + } // end of plus loop + + free(layer_count); + free(tRef); + + // Return reflections (minus and plus) + return; +} + diff --git a/empymod/transform.py b/empymod/transform.py index 4235f51..758eaac 100644 --- a/empymod/transform.py +++ b/empymod/transform.py @@ -26,16 +26,15 @@ # WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the # License for the specific language governing permissions and limitations under # the License. - - import numpy as np import scipy as sp -from empymod import kernel +from empymod import kernel, _ckernel + -__all__ = ['hankel_dlf', 'hankel_qwe', 'hankel_quad', 'fourier_dlf', - 'fourier_qwe', 'fourier_fftlog', 'fourier_fft', 'dlf', 'qwe', - 'get_dlf_points', 'get_fftlog_input'] +__all__ = ['hankel_dlf', 'hankel_cdlf', 'hankel_qwe', 'hankel_quad', + 'fourier_dlf', 'fourier_qwe', 'fourier_fftlog', 'fourier_fft', + 'dlf', 'qwe', 'get_dlf_points', 'get_fftlog_input'] def __dir__(): @@ -119,6 +118,28 @@ def hankel_dlf(zsrc, zrec, lsrc, lrec, off, ang_fact, depth, ab, etaH, etaV, return fEM, 1, True +def hankel_cdlf(zsrc, zrec, lsrc, lrec, off, ang_fact, depth, ab, etaH, etaV, + zetaH, zetaV, xdirect, htarg, msrc, mrec): + r"""Hankel Transform using the Digital Linear Filter method. + + This function is the same as :func:`hankel_dlf`, but uses a C-version of + the kernel. + """ + + # Compute required lambdas for given Hankel-filter-base + lambd, int_pts = get_dlf_points(htarg['dlf'], off, htarg['pts_per_dec']) + + # Call the kernel + PJ = _ckernel.wavenumber(zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH, + zetaV, lambd, ab, xdirect, msrc, mrec) + + # Carry out the dlf + fEM = dlf(PJ, lambd, off, htarg['dlf'], htarg['pts_per_dec'], + ang_fact=ang_fact, ab=ab, int_pts=int_pts) + + return fEM, 1, True + + def hankel_qwe(zsrc, zrec, lsrc, lrec, off, ang_fact, depth, ab, etaH, etaV, zetaH, zetaV, xdirect, htarg, msrc, mrec): r"""Hankel Transform using Quadrature-With-Extrapolation. diff --git a/empymod/utils.py b/empymod/utils.py index 9247205..36c8984 100644 --- a/empymod/utils.py +++ b/empymod/utils.py @@ -513,7 +513,7 @@ def check_hankel(ht, htarg, verb): targ = {} args = copy.deepcopy(htarg) - if ht == 'dlf': # DLF + if ht in ['dlf', 'cdlf']: # DLF # If filter is a name (str), get it targ['dlf'] = args.pop('dlf', filters.Hankel().key_201_2009) @@ -883,7 +883,7 @@ def check_loop(loop, ht, htarg, verb): # Define if to loop over frequencies or over offsets lagged_splined_dlf = False - if ht == 'dlf': + if ht in ['dlf', 'cdlf']: if htarg['pts_per_dec'] != 0: lagged_splined_dlf = True diff --git a/empymod/wavenumber.c b/empymod/wavenumber.c new file mode 100644 index 0000000..37b2983 --- /dev/null +++ b/empymod/wavenumber.c @@ -0,0 +1,125 @@ +#include +#include +#include +#include +#include +#include +#include + +// Function prototypes +void greenfct(int nfreq, int noff, int nlayer, int nlambda, double zsrc, double zrec, int lsrc, int lrec, double* depth, double complex* etaH, double complex* etaV, double complex* zetaH, double complex* zetaV, double* lambd, int ab, int xdirect, int msrc, int mrec, double complex* GTM, double complex* GTE); + +void wavenumber(int nfreq, int noff, int nlayer, int nlambda, + double zsrc, double zrec, int lsrc, int lrec, double *depth, + double complex *etaH, double complex *etaV, double complex *zetaH, double complex *zetaV, + double *lambd, int ab, int xdirect, int msrc, int mrec, + double complex *PJ0, double complex *PJ1, double complex *PJ0b) { + int g1, g2, g3, r1, r2; + int i, ii, iv; + double fourpi, eightpi, sign, dlambd, tlambd; + double complex *PTM, *PTE, Ptot; + + // Calculate Green's functions + g1 = nlambda; + g2 = g1*nlayer; + g3 = g2*noff; + + r1=nlambda; + r2=noff*r1; + + PTM = (double complex *) calloc(nfreq * noff * nlambda , sizeof(double complex)); + PTE = (double complex *) calloc(nfreq * noff * nlambda , sizeof(double complex)); + + greenfct(nfreq, noff, nlayer, nlambda, zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH, zetaV, lambd, + ab, xdirect, msrc, mrec, PTM, PTE); + + // Pre-allocate output + if (ab == 11 || ab == 22 || ab == 24 || ab == 15 || ab == 33) { + PJ0 = PJ0; + } else { + PJ0 = NULL; + } + if (ab == 11 || ab == 12 || ab == 21 || ab == 22 || ab == 14 || ab == 24 || + ab == 15 || ab == 25) { + PJ0b = PJ0b; + } else { + PJ0b = NULL; + } + if (ab != 33) { + PJ1 = PJ1; + } else { + PJ1 = NULL; + } + + + //fourpi = 1.0 /(4 * 3.14159265358979323846); + fourpi = 1.0 /(4 * M_PI); + // If rec is magnetic switch sign (reciprocity MM/ME => EE/EM) + if (mrec) { + sign = -1.0; + } else { + sign = 1.0; + } + + // Group into PJ0 and PJ1 for J0/J1 Hankel Transform + if (ab == 11 || ab == 12 || ab == 21 || ab == 22 || ab == 14 || ab == 24 || + ab == 15 || ab == 25) { + if (ab == 14 || ab == 22) { + sign *= -1; + } + + for (i = 0; i < nfreq; i++) { + for (ii = 0; ii < noff; ii++) { + for (iv = 0; iv < nlambda; iv++) { + Ptot = (PTM[i*r2+ii*r1+iv] + PTE[i*r2+ii*r1+iv]) * fourpi; + PJ0b[i*r2+ii*r1+iv] = sign * 0.5 * Ptot * lambd[ii*nlambda+iv]; + PJ1[i*r2+ii*r1+iv] = -1.0*sign * Ptot; + } + } + } + + if (ab == 11 || ab == 22 || ab == 24 || ab == 15) { + if (ab == 22 || ab == 24) { + sign *= -1; + } + + eightpi = sign / (8 * M_PI); + for (i = 0; i < nfreq; i++) { + for (ii = 0; ii < noff; ii++) { + for (iv = 0; iv < nlambda; iv++) { + PJ0[i*r2+ii*r1+iv] = (PTM[i*r2+ii*r1+iv] - PTE[i*r2+ii*r1+iv])*lambd[ii*nlambda+iv] * eightpi; + } + } + } + } + + } else if (ab == 13 || ab == 23 || ab == 31 || ab == 32 || ab == 34 || + ab == 35 || ab == 16 || ab == 26) { + if (ab == 34 || ab == 26) { + sign *= -1; + } + for (i = 0; i < nfreq; i++) { + for (ii = 0; ii < noff; ii++) { + for (iv = 0; iv < nlambda; iv++) { + dlambd = lambd[ii*nlambda+iv] * lambd[ii*nlambda+iv]; + Ptot = (PTM[i*r2+ii*r1+iv] + PTE[i*r2+ii*r1+iv]) * fourpi; + PJ1[i*r2+ii*r1+iv] = sign * Ptot * dlambd; + } + } + } + } else if (ab == 33) { + for (i = 0; i < nfreq; i++) { + for (ii = 0; ii < noff; ii++) { + for (iv = 0; iv < nlambda; iv++) { + tlambd = lambd[ii*nlambda+iv] * lambd[ii*nlambda+iv] * lambd[ii*nlambda+iv]; + Ptot = (PTM[i*r2+ii*r1+iv] + PTE[i*r2+ii*r1+iv]) * fourpi; + PJ0[i*r2+ii*r1+iv] = sign * Ptot * tlambd; + } + } + } + } + free(PTM); + free(PTE); + return; +} + diff --git a/tests/test_kernel.py b/tests/test_kernel.py index fcf458b..479692b 100644 --- a/tests/test_kernel.py +++ b/tests/test_kernel.py @@ -3,7 +3,7 @@ from os.path import join, dirname from numpy.testing import assert_allclose -from empymod import kernel +from empymod import kernel, _ckernel from empymod import bipole # No input checks are carried out in kernel, by design. Input checks are @@ -21,9 +21,11 @@ allow_pickle=True) -@pytest.mark.parametrize("njit", [True, False]) +@pytest.mark.parametrize("njit", [True, False, 'c']) def test_wavenumber(njit): # 1. wavenumber - if njit: + if njit == 'c': + wavenumber = _ckernel.wavenumber + elif njit: wavenumber = kernel.wavenumber else: wavenumber = kernel.wavenumber.py_func @@ -48,9 +50,11 @@ def test_wavenumber(njit): # 1. wavenumber assert out[2] is None -@pytest.mark.parametrize("njit", [True, False]) +@pytest.mark.parametrize("njit", [True, False, 'c']) def test_greenfct(njit): # 2. greenfct - if njit: + if njit == 'c': + greenfct = _ckernel.greenfct + elif njit: greenfct = kernel.greenfct else: greenfct = kernel.greenfct.py_func @@ -66,9 +70,11 @@ def test_greenfct(njit): # 2. greenfct assert_allclose(out[1], val[i+1][1]) -@pytest.mark.parametrize("njit", [True, False]) +@pytest.mark.parametrize("njit", [True, False, 'c']) def test_reflections(njit): # 3. reflections - if njit: + if njit == 'c': + reflections = _ckernel.reflections + elif njit: reflections = kernel.reflections else: reflections = kernel.reflections.py_func @@ -80,9 +86,11 @@ def test_reflections(njit): # 3. reflections assert_allclose(Rm, val[2]) -@pytest.mark.parametrize("njit", [True, False]) +@pytest.mark.parametrize("njit", [True, False, 'c']) def test_fields(njit): # 4. fields - if njit: + if njit == 'c': + fields = _ckernel.fields + elif njit: fields = kernel.fields else: fields = kernel.fields.py_func diff --git a/tests/test_model.py b/tests/test_model.py index 3a077d6..c7b03a7 100644 --- a/tests/test_model.py +++ b/tests/test_model.py @@ -12,6 +12,10 @@ from empymod.model import gpr, dipole_k, fem, tem from empymod.kernel import fullspace, halfspace +# NOTE: To test the c-kernel, uncomment the following two lines. +# from empymod import transform +# transform.hankel_dlf = transform.hankel_cdlf + # These are kind of macro-tests, as they check the final results. # I try to use different parameters for each test, to cover a wide range of # possibilities. It won't be possible to check all the possibilities though. diff --git a/tests/test_transform.py b/tests/test_transform.py index 2ccb08e..146295e 100644 --- a/tests/test_transform.py +++ b/tests/test_transform.py @@ -16,7 +16,7 @@ allow_pickle=True) -@pytest.mark.parametrize("htype", ['dlf', 'qwe', 'quad']) +@pytest.mark.parametrize("htype", ['cdlf', 'dlf', 'qwe', 'quad']) def test_hankel(htype): # 1. DLF / 2. QWE / 3. QUAD # Compare wavenumber-domain calculation / DLF with analytical # frequency-domain fullspace solution