C Calls

In some production workflows, it is necessary to introduce additional apparatus to an FD kernel to carry out functions such as Fourier transforms or status updates. Due to the wide variety of possibilities, DevitoPRO features an API for introducing arbitrary C calls to the generated code, enabling operators to be customized as needed.

In this notebook we show how to add custom C calls to an Operator specified in DevitoPRO.

#NBVAL_IGNORE_OUTPUT
import os

import matplotlib.pyplot as plt
import numpy as np

from scipy.signal import hilbert
from IPython.display import Code

from devito import Eq, Grid, Operator, Function, TimeFunction, Dimension, solve, ConditionalDimension
from devito.types import Symbol, Array
from devito.tools import CustomDtype
from devitopro import *
from devitopro.types.enriched import NoLayers

from examples.seismic import TimeAxis, RickerSource, demo_model
Submodule devito version mismatch. Installed: 4.8.14+40.g0c9f18b21.dirty, Required: 4.8.20+62.g74903be47

Inserting simple print statements

In this first example, we will insert some print statements within the kernel to report the timestep and index being calculated.

To insert a custom C call, we use the CCall class. Specifically, we subclass it, adding name and header attributes, corresponding to the function to be called, and the header file from which it is to be included.

class FPrintF(CCall):
    name = 'fprintf'
    header = 'stdio.h'

To minimise the amount of output, we will define a minimal Grid.

grid = Grid(shape=(3,), extent=(1.,))
u = TimeFunction(name='u', grid=grid, dtype=np.int32)

Now we will define a trivial operator, which simply adds 1 to the field values at each iteration. We also insert print statements to report the timestep and index before each update.

Note that in this case we are printing to the stderr to simplify output routing in the Jupyter notebook; stdout is not always written to the notebook.

time = grid.time_dim
x = grid.dimensions[0]

# Our C call statements are included as if they were Eq, albeit passing a list of arguments
# Note the '\\n' for newline. The additional escape slash is needed for correct printing.
op = Operator([FPrintF(['stderr', '"Out: Timestep %d\\n"', time]),
               # implicit_dims=time used to nest this loop within the t, x loop nest
               # If it is not known to be time-dependent, then the Devito compiler will extract it
               # into its own loop
               FPrintF(['stderr', '"Out: Index %d\\n"', x], implicit_dims=time),
               Eq(u.forward, u + 1.)])

Printing the generated C code, we observe that the specificed print statements have been inserted.

#NBVAL_IGNORE_OUTPUT
Code(str(op.ccode), language='C')
/* NON-COMMERCIAL, NON-EXCLUSIVE, NON-TRANSFERABLE SOFTWARE LICENSE

Copyright (c) 2021-24, Devito Codes Ltd, 1st Floor, Commerce House Raven Road,
London E18 1HB, United Kingdom.

Definitions:

* "Devito Codes": Refers to Devito Codes Ltd., a company incorporated and operating under the laws of England and Wales, and encompasses all its divisions, affiliates, and subsidiaries involved in the development and licensing of software products.
* "Software": Includes all software products developed by Devito Codes, such as Devito PRO, and any other software titles offered by Devito Codes, along with their associated materials, documentation, and software generated by their respective compilers or runtime environments.
* "Licensee": Denotes any individual, company, or entity that has legally acquired a copy or license of the Software under the terms specified in an Agreement with Devito Codes. This term also encompasses all employees, contractors, and agents acting within the scope of their engagement by the Licensee.
* "Agreement": Refers to the legally binding license agreement executed between the Licensee and Devito Codes, which outlines the terms and conditions under which the Licensee is authorized to use the Software. This Agreement may be specific to a single software product or may encompass multiple products offered by Devito Codes.
* "Protected Data/Software": Encompasses all data, software, algorithms, source and object code, documentation, proprietary methodologies, and intellectual property owned or developed by Devito Codes. This includes, without limitation, all versions and iterations of Devito PRO software, materials, and any output generated by Devito Codes' compilers or software, as well as any other proprietary software and tools developed by Devito Codes.
* "Unauthorized Use": Constitutes any action involving the Protected Data/Software that is not expressly permitted under the terms of the Agreement. This includes, but is not limited to, copying, modifying, distributing, selling, leasing, reverse engineering, decompiling, or incorporating the Protected Data/Software into any external system or software, including data retrieval systems, databases, generative AI models, machine learning models, Retrieval Augmented Generation (RAG) applications, or any form of technology or AI system, without the explicit written consent of Devito Codes.
* "Non-Commercial Use" refers to activities not intended for or directed towards commercial advantage or monetary compensation. This includes, but is not limited to, academic research, educational purposes, and personal projects.

Devito Codes provides the Software on the following terms:

1. Scope of License: The Licensee is hereby granted a non-commercial, non-exclusive, non-transferable license without the right to sublicense the Protected Data/Software for the purpose of internal evaluation, testing, and development within the Licensee's organization. This license explicitly excludes any commercial applications, product development for sale, or external consulting services.
2. Prohibition on Unauthorized Use: The Licensee shall not use, adapt, modify, or incorporate the Protected Data/Software into any software or hardware that is not explicitly authorized by this Agreement. This includes prohibitions against the use of the Protected Data/Software in the development, training, or operation of generative AI, machine learning models, RAG applications, or any form of AI that utilizes external data sources for augmentation without the prior written consent of Devito Codes.
3. Distribution and External Use: The Licensee agrees not to distribute, share, lease, or otherwise make the Protected Data/Software available to any third party, including through cloud services, Software as a Service (SaaS), or any form of redistribution mechanism. The Licensee also agrees not to use the Protected Data/Software to provide any service to external parties that involves the processing of data or computation using the Protected Data/Software.
4. The Licensee is prohibited from reverse engineering, decompiling, disassembling, or otherwise attempting to discover the source code of the Protected Data/Software. This prohibition extends to any form of reverse engineering or hacking intended to access or create unauthorized derivatives of the Software.
5. The Software is experimental and licensed "as is." The license of the Software does not include technical support.
6. The Licensee agrees that it will not license or sell the Software or any other Software, information, or data that incorporate any part of the Software, including derivative works ("Derivatives"), to any other parties.
7. Devito Codes will consider all requests for a commercial license but shall be under no obligation to grant such a license
8. Licensee agrees that any person within the Licensee utilizing the Software will be advised of, and is subject to, the conditions in the License.
9. THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR NON-INFRINGEMENT. Devito Codes Ltd expressly disclaims all warranties, whether express, implied, statutory, or otherwise, with respect to the software provided hereunder, including any warranty of reliability, accuracy, or the absence of errors or defects, and the implied warranty arising from course of performance, course of dealing, or usage of trade. Devito Codes Ltd shall not be liable for any direct, indirect, incidental, special, consequential, or exemplary damages, including but not limited to, damages for loss of profits, goodwill, use, data, or other intangible losses (even if Devito Codes Ltd has been advised of the possibility of such damages), resulting from the use or the inability to use the Software. The Licensee acknowledges and agrees that the entire risk arising out of the use or performance of the Software remains with the Licensee. Notwithstanding any damages that the Licensee might incur for any reason whatsoever (including, without limitation, all damages referenced above and all direct or general damages), the entire liability of Devito Codes Ltd and any of its suppliers under any provision of this License and the Licensee's exclusive remedy for all of the foregoing shall be limited to the amount actually paid by the Licensee for the Software. The foregoing limitations, exclusions, and disclaimers shall apply to the maximum extent permitted by applicable law, even if any remedy fails its essential purpose.
10. All rights, title, and interest in and to all data, information, and inventions that result from the Software by the Licensee shall vest in and belong to the Licensee. Licensee agrees to provide Devito Codes with copies of publications referencing the Software and acknowledge Devito Codes in those publications.
11. The Licensee agrees to notify Devito Codes in writing of any known or suspected distribution or use of the Software not in compliance with the License requirements and to enforce the terms of the License.
12. The Licensee agrees to indemnify and hold harmless Devito Codes against any claims or damages arising from the Licensee's use of the Software, exceeding the scope of this License. This includes any third-party claims related to the Software's use.
13. Intellectual Property Rights: All intellectual property rights in the Protected Data/Software and any derivatives thereof shall remain with Devito Codes. The Licensee acknowledges that no ownership rights are transferred by this Agreement.
14. Compliance with Future Technologies: The Licensee commits to continuously ensuring that the use of the Software, including in the context of emerging technologies and computational techniques, adheres to the terms outlined in this License. Devito Codes will provide guidelines for compliance as new technologies develop, facilitating ongoing lawful use of the Software.
15. Notification of Unauthorized Use: The Licensee shall promptly notify Devito Codes of any unauthorized use, copying, or distribution of the Protected Data/Software, or any breach of this License, of which it becomes aware.
16. Where a clause in this License and the Agreement conflict, the clause in the Agreement will take precedence.
 */

#define _POSIX_C_SOURCE 200809L
#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);
#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;
#define uL0(t,x) u[(t)*x_stride0 + (x)]

#include "stdlib.h"
#include "math.h"
#include "sys/time.h"
#include "xmmintrin.h"
#include "pmmintrin.h"
#include "stdio.h"
#include "omp.h"

struct dataobj
{
  void *restrict data;
  int * size;
  unsigned long nbytes;
  unsigned long * npsize;
  unsigned long * dsize;
  int * hsize;
  int * hofs;
  int * oofs;
  void * dmap;
} ;

struct profiler
{
  double section0;
} ;


int Kernel(struct dataobj *restrict u_vec, const int time_M, const int time_m, const int x_M, const int x_m, const int nthreads, struct profiler * timers)
{
  int *u __attribute__ ((aligned (64))) = (int *) u_vec->data;

  const int x_fsz0 = u_vec->size[1];

  const int x_stride0 = x_fsz0;

  /* Flush denormal numbers to zero in hardware */
  _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
  _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);

  for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
  {
    fprintf(stderr,"Out: Timestep %d\n",time);

    START(section0)
    #pragma omp parallel num_threads(nthreads)
    {
      #pragma omp for schedule(static,1)
      for (int x = x_m; x <= x_M; x += 1)
      {
        fprintf(stderr,"Out: Index %d\n",x);
        uL0(t1, x + 1) = uL0(t0, x + 1) + 1.0;
      }
    }
    STOP(section0,timers)
  }

  return 0;
}

Running this Operator, we observe the timesteps and indices outputted to stderr. The indices are out of order due to OpenMP parallelism.

#NBVAL_IGNORE_OUTPUT
op(time_M=5)
Out: Timestep 0
Out: Index 0
Operator `Kernel` ran in 0.01 s
Out: Index 1
Out: Index 2
Out: Timestep 1
Out: Index 0
Out: Index 1
Out: Index 2
Out: Timestep 2
Out: Index 0
Out: Index 2
Out: Index 1
Out: Timestep 3
Out: Index 0
Out: Index 2
Out: Index 1
Out: Timestep 4
Out: Index 0
Out: Index 1
Out: Index 2
Out: Timestep 5
Out: Index 0
Out: Index 1
Out: Index 2
PerformanceSummary([(PerfKey(name='section0', rank=None),
                     PerfEntry(time=0.0009200000000000001, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])

Inserting a call from a custom header

We will now consider the insertion of a C call from a custom header file.

In this case, the header file (fizzbuzz.h) contains a simple “FizzBuzz” function, which takes an int, and returns an appropriate string. We will include this in the Devito kernel by setting a Function equal to the sum of the time and space indices, then feeding the Function value at each point to the FizzBuzz. The string will be assigned to a value which is then printed.

First, we need to make the C call corresponding to our function and specify the correct directory to search.

class FizzBuzz(CCall):
    name = 'fizzBuzz'
    header = 'fizzbuzz.h'
    # Search the directory containing this notebook
    header_dirs = (os.path.abspath(""),)

Note that the function fizzBuzz consists of the following:

#NBVAL_IGNORE_OUTPUT
contents = !cat fizzbuzz.h
Code('\n'.join(contents), language='C')
#include <stdio.h>

// Takes an integer and returns the corresponding FizzBuzz
const char * fizzBuzz(int v)
{
    if (v % 3 == 0)
    {
        if (v % 5 == 0)
        {
            return "fizzbuzz";
        }
        return "fizz";
    } else if (v % 5 == 0)
    {
        return "buzz";
    }
    char out[3];
    sprintf(out, "%d", v);
    const char * ptr = out;
    return ptr;
}

Now we build our Operator. As a C string is an array of char, we will need to create a CustomDtype for a pointer to a char array. We can then create a Symbol using this dtype, to which we can assign the return value of fizzBuzz.

# Make a char * dtype
charstar = CustomDtype(name='char', modifier=' *')

u = TimeFunction(name='u', grid=grid, dtype=np.int32)
s = Symbol(name='s', dtype=charstar, is_const=True)

op1 = Operator([Eq(u, time + x),
                FizzBuzz([u], s),
                # Note the ordering of implicit dims needs to match the grid dimension order
                FPrintF(['stderr', '"Out: %s\\n"', s], implicit_dims=(time, x))])

Printing the generated C code, we can see both the include at the top, and the function call within the loop.

#NBVAL_IGNORE_OUTPUT
Code(str(op1.ccode), language='C')
/* NON-COMMERCIAL, NON-EXCLUSIVE, NON-TRANSFERABLE SOFTWARE LICENSE

Copyright (c) 2021-24, Devito Codes Ltd, 1st Floor, Commerce House Raven Road,
London E18 1HB, United Kingdom.

Definitions:

* "Devito Codes": Refers to Devito Codes Ltd., a company incorporated and operating under the laws of England and Wales, and encompasses all its divisions, affiliates, and subsidiaries involved in the development and licensing of software products.
* "Software": Includes all software products developed by Devito Codes, such as Devito PRO, and any other software titles offered by Devito Codes, along with their associated materials, documentation, and software generated by their respective compilers or runtime environments.
* "Licensee": Denotes any individual, company, or entity that has legally acquired a copy or license of the Software under the terms specified in an Agreement with Devito Codes. This term also encompasses all employees, contractors, and agents acting within the scope of their engagement by the Licensee.
* "Agreement": Refers to the legally binding license agreement executed between the Licensee and Devito Codes, which outlines the terms and conditions under which the Licensee is authorized to use the Software. This Agreement may be specific to a single software product or may encompass multiple products offered by Devito Codes.
* "Protected Data/Software": Encompasses all data, software, algorithms, source and object code, documentation, proprietary methodologies, and intellectual property owned or developed by Devito Codes. This includes, without limitation, all versions and iterations of Devito PRO software, materials, and any output generated by Devito Codes' compilers or software, as well as any other proprietary software and tools developed by Devito Codes.
* "Unauthorized Use": Constitutes any action involving the Protected Data/Software that is not expressly permitted under the terms of the Agreement. This includes, but is not limited to, copying, modifying, distributing, selling, leasing, reverse engineering, decompiling, or incorporating the Protected Data/Software into any external system or software, including data retrieval systems, databases, generative AI models, machine learning models, Retrieval Augmented Generation (RAG) applications, or any form of technology or AI system, without the explicit written consent of Devito Codes.
* "Non-Commercial Use" refers to activities not intended for or directed towards commercial advantage or monetary compensation. This includes, but is not limited to, academic research, educational purposes, and personal projects.

Devito Codes provides the Software on the following terms:

1. Scope of License: The Licensee is hereby granted a non-commercial, non-exclusive, non-transferable license without the right to sublicense the Protected Data/Software for the purpose of internal evaluation, testing, and development within the Licensee's organization. This license explicitly excludes any commercial applications, product development for sale, or external consulting services.
2. Prohibition on Unauthorized Use: The Licensee shall not use, adapt, modify, or incorporate the Protected Data/Software into any software or hardware that is not explicitly authorized by this Agreement. This includes prohibitions against the use of the Protected Data/Software in the development, training, or operation of generative AI, machine learning models, RAG applications, or any form of AI that utilizes external data sources for augmentation without the prior written consent of Devito Codes.
3. Distribution and External Use: The Licensee agrees not to distribute, share, lease, or otherwise make the Protected Data/Software available to any third party, including through cloud services, Software as a Service (SaaS), or any form of redistribution mechanism. The Licensee also agrees not to use the Protected Data/Software to provide any service to external parties that involves the processing of data or computation using the Protected Data/Software.
4. The Licensee is prohibited from reverse engineering, decompiling, disassembling, or otherwise attempting to discover the source code of the Protected Data/Software. This prohibition extends to any form of reverse engineering or hacking intended to access or create unauthorized derivatives of the Software.
5. The Software is experimental and licensed "as is." The license of the Software does not include technical support.
6. The Licensee agrees that it will not license or sell the Software or any other Software, information, or data that incorporate any part of the Software, including derivative works ("Derivatives"), to any other parties.
7. Devito Codes will consider all requests for a commercial license but shall be under no obligation to grant such a license
8. Licensee agrees that any person within the Licensee utilizing the Software will be advised of, and is subject to, the conditions in the License.
9. THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR NON-INFRINGEMENT. Devito Codes Ltd expressly disclaims all warranties, whether express, implied, statutory, or otherwise, with respect to the software provided hereunder, including any warranty of reliability, accuracy, or the absence of errors or defects, and the implied warranty arising from course of performance, course of dealing, or usage of trade. Devito Codes Ltd shall not be liable for any direct, indirect, incidental, special, consequential, or exemplary damages, including but not limited to, damages for loss of profits, goodwill, use, data, or other intangible losses (even if Devito Codes Ltd has been advised of the possibility of such damages), resulting from the use or the inability to use the Software. The Licensee acknowledges and agrees that the entire risk arising out of the use or performance of the Software remains with the Licensee. Notwithstanding any damages that the Licensee might incur for any reason whatsoever (including, without limitation, all damages referenced above and all direct or general damages), the entire liability of Devito Codes Ltd and any of its suppliers under any provision of this License and the Licensee's exclusive remedy for all of the foregoing shall be limited to the amount actually paid by the Licensee for the Software. The foregoing limitations, exclusions, and disclaimers shall apply to the maximum extent permitted by applicable law, even if any remedy fails its essential purpose.
10. All rights, title, and interest in and to all data, information, and inventions that result from the Software by the Licensee shall vest in and belong to the Licensee. Licensee agrees to provide Devito Codes with copies of publications referencing the Software and acknowledge Devito Codes in those publications.
11. The Licensee agrees to notify Devito Codes in writing of any known or suspected distribution or use of the Software not in compliance with the License requirements and to enforce the terms of the License.
12. The Licensee agrees to indemnify and hold harmless Devito Codes against any claims or damages arising from the Licensee's use of the Software, exceeding the scope of this License. This includes any third-party claims related to the Software's use.
13. Intellectual Property Rights: All intellectual property rights in the Protected Data/Software and any derivatives thereof shall remain with Devito Codes. The Licensee acknowledges that no ownership rights are transferred by this Agreement.
14. Compliance with Future Technologies: The Licensee commits to continuously ensuring that the use of the Software, including in the context of emerging technologies and computational techniques, adheres to the terms outlined in this License. Devito Codes will provide guidelines for compliance as new technologies develop, facilitating ongoing lawful use of the Software.
15. Notification of Unauthorized Use: The Licensee shall promptly notify Devito Codes of any unauthorized use, copying, or distribution of the Protected Data/Software, or any breach of this License, of which it becomes aware.
16. Where a clause in this License and the Agreement conflict, the clause in the Agreement will take precedence.
 */

#define _POSIX_C_SOURCE 200809L
#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);
#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;
#define uL0(t,x) u[(t)*x_stride0 + (x)]

#include "stdlib.h"
#include "math.h"
#include "sys/time.h"
#include "xmmintrin.h"
#include "pmmintrin.h"
#include "fizzbuzz.h"
#include "stdio.h"
#include "omp.h"

struct dataobj
{
  void *restrict data;
  int * size;
  unsigned long nbytes;
  unsigned long * npsize;
  unsigned long * dsize;
  int * hsize;
  int * hofs;
  int * oofs;
  void * dmap;
} ;

struct profiler
{
  double section0;
} ;


int Kernel(struct dataobj *restrict u_vec, const int time_M, const int time_m, const int x_M, const int x_m, const int nthreads, struct profiler * timers)
{
  int *u __attribute__ ((aligned (64))) = (int *) u_vec->data;

  const int x_fsz0 = u_vec->size[1];

  const int x_stride0 = x_fsz0;

  /* Flush denormal numbers to zero in hardware */
  _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
  _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);

  for (int time = time_m, t0 = (time)%(2); time <= time_M; time += 1, t0 = (time)%(2))
  {
    START(section0)
    #pragma omp parallel num_threads(nthreads)
    {
      #pragma omp for schedule(static,1)
      for (int x = x_m; x <= x_M; x += 1)
      {
        uL0(t0, x + 1) = x + time;
        const char * s = fizzBuzz(uL0(t0, x + 1));
        fprintf(stderr,"Out: %s\n",s);
      }
    }
    STOP(section0,timers)
  }

  return 0;
}

Running the Operator, we obtain the FizzBuzz output. Note that numbers are repeated due to the use of x + time as the input for the call, combined with the use of #pragma omp for schedule(static,1) over x, which causes output to not increase monotonically, as in previous examples.

#NBVAL_IGNORE_OUTPUT
op1(time_M=15)
Out: fizzbuzz
Out: (null)
Out: (null)
Out: fizz
Out: (null)
Out: (null)
Out: fizz
Out: (null)
Out: (null)
Out: buzz
Out: (null)
Out: fizz
Out: buzz
Out: fizz
Out: (null)
Out: buzz
Out: (null)
Out: fizz
Out: fizz
Out: (null)
Out: (null)
Out: fizz
Out: (null)
Out: (null)
Out: (null)
Out: buzz
Out: fizz
Out: fizz
Out: (null)
Out: buzz
Out: buzz
Out: (null)
Out: fizz
Out: (null)
Out: (null)
Out: fizz
Out: fizz
Out: (null)
Out: (null)
Out: (null)
Out: (null)
Out: fizzbuzz
Out: fizzbuzz
Out: (null)
Out: (null)
Out: fizzbuzz
Out: (null)
Out: (null)
Operator `Kernel` ran in 0.01 s
PerformanceSummary([(PerfKey(name='section0', rank=None),
                     PerfEntry(time=0.0005709999999999999, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])

Inserting a call from a library: Part 1

The CCall API also allows the calling of functions from libraries such as BLAS. To demonstrate a simple example, we will modify a trivial kernel setting values in two Functions, to print the dot product of the two resultant underlying arrays to stderr.

For this purpose, we will use cblas_sdot, which takes a single-precision dot product of two arrays. Note that the default dtype for Functions is np.float32. The BLAS operator needs to be matched to the precision of the values to which it applied.

# Find path to BLAS (assumed to be OpenBlas, change as needed or manually specify path)
# Selecting the correct directory is left as an exercise for the user, due to the wide range
# of installation paths and operating systems
# blasdir=''

# On ubntu installed via `apt install libopenblas-dev` and similarly on rhel/centos with yum

blasdir = '/usr'
# For brew install
# blasdir = !brew --prefix openblas
# blasdir = blasdir[0]  # Required for MacOS command
class BlasCall(CCall):
    """Class for calls to BLAS"""
    header = 'cblas.h'
    header_dirs = (f'{blasdir}/include/',)
    libs = ('openblas',)
    lib_dirs = (f'{blasdir}/lib',)


class SDot(BlasCall):
    name = 'cblas_sdot'

As before, we insert our custom C calls into the Operator.

u = Function(name='u', grid=grid)
v = Function(name='v', grid=grid)
s = Symbol(name='s', dtype=np.float32)

op2 = Operator([Eq(u, x),
                Eq(v, 2*x + 1),
                # .base gives the whole array, rather than an indexed value in the array 
                SDot([u.shape_with_halo[0], u.base, 1, v.base, 1], s),
                # Note the ordering of implicit dims needs to match the grid dimension order
                FPrintF(['stderr', '"Out: %f\\n"', s])])

Printing the C code, we observe both the BLAS call and the print statement.

#NBVAL_IGNORE_OUTPUT
Code(str(op2.ccode), language='C')
/* NON-COMMERCIAL, NON-EXCLUSIVE, NON-TRANSFERABLE SOFTWARE LICENSE

Copyright (c) 2021-24, Devito Codes Ltd, 1st Floor, Commerce House Raven Road,
London E18 1HB, United Kingdom.

Definitions:

* "Devito Codes": Refers to Devito Codes Ltd., a company incorporated and operating under the laws of England and Wales, and encompasses all its divisions, affiliates, and subsidiaries involved in the development and licensing of software products.
* "Software": Includes all software products developed by Devito Codes, such as Devito PRO, and any other software titles offered by Devito Codes, along with their associated materials, documentation, and software generated by their respective compilers or runtime environments.
* "Licensee": Denotes any individual, company, or entity that has legally acquired a copy or license of the Software under the terms specified in an Agreement with Devito Codes. This term also encompasses all employees, contractors, and agents acting within the scope of their engagement by the Licensee.
* "Agreement": Refers to the legally binding license agreement executed between the Licensee and Devito Codes, which outlines the terms and conditions under which the Licensee is authorized to use the Software. This Agreement may be specific to a single software product or may encompass multiple products offered by Devito Codes.
* "Protected Data/Software": Encompasses all data, software, algorithms, source and object code, documentation, proprietary methodologies, and intellectual property owned or developed by Devito Codes. This includes, without limitation, all versions and iterations of Devito PRO software, materials, and any output generated by Devito Codes' compilers or software, as well as any other proprietary software and tools developed by Devito Codes.
* "Unauthorized Use": Constitutes any action involving the Protected Data/Software that is not expressly permitted under the terms of the Agreement. This includes, but is not limited to, copying, modifying, distributing, selling, leasing, reverse engineering, decompiling, or incorporating the Protected Data/Software into any external system or software, including data retrieval systems, databases, generative AI models, machine learning models, Retrieval Augmented Generation (RAG) applications, or any form of technology or AI system, without the explicit written consent of Devito Codes.
* "Non-Commercial Use" refers to activities not intended for or directed towards commercial advantage or monetary compensation. This includes, but is not limited to, academic research, educational purposes, and personal projects.

Devito Codes provides the Software on the following terms:

1. Scope of License: The Licensee is hereby granted a non-commercial, non-exclusive, non-transferable license without the right to sublicense the Protected Data/Software for the purpose of internal evaluation, testing, and development within the Licensee's organization. This license explicitly excludes any commercial applications, product development for sale, or external consulting services.
2. Prohibition on Unauthorized Use: The Licensee shall not use, adapt, modify, or incorporate the Protected Data/Software into any software or hardware that is not explicitly authorized by this Agreement. This includes prohibitions against the use of the Protected Data/Software in the development, training, or operation of generative AI, machine learning models, RAG applications, or any form of AI that utilizes external data sources for augmentation without the prior written consent of Devito Codes.
3. Distribution and External Use: The Licensee agrees not to distribute, share, lease, or otherwise make the Protected Data/Software available to any third party, including through cloud services, Software as a Service (SaaS), or any form of redistribution mechanism. The Licensee also agrees not to use the Protected Data/Software to provide any service to external parties that involves the processing of data or computation using the Protected Data/Software.
4. The Licensee is prohibited from reverse engineering, decompiling, disassembling, or otherwise attempting to discover the source code of the Protected Data/Software. This prohibition extends to any form of reverse engineering or hacking intended to access or create unauthorized derivatives of the Software.
5. The Software is experimental and licensed "as is." The license of the Software does not include technical support.
6. The Licensee agrees that it will not license or sell the Software or any other Software, information, or data that incorporate any part of the Software, including derivative works ("Derivatives"), to any other parties.
7. Devito Codes will consider all requests for a commercial license but shall be under no obligation to grant such a license
8. Licensee agrees that any person within the Licensee utilizing the Software will be advised of, and is subject to, the conditions in the License.
9. THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR NON-INFRINGEMENT. Devito Codes Ltd expressly disclaims all warranties, whether express, implied, statutory, or otherwise, with respect to the software provided hereunder, including any warranty of reliability, accuracy, or the absence of errors or defects, and the implied warranty arising from course of performance, course of dealing, or usage of trade. Devito Codes Ltd shall not be liable for any direct, indirect, incidental, special, consequential, or exemplary damages, including but not limited to, damages for loss of profits, goodwill, use, data, or other intangible losses (even if Devito Codes Ltd has been advised of the possibility of such damages), resulting from the use or the inability to use the Software. The Licensee acknowledges and agrees that the entire risk arising out of the use or performance of the Software remains with the Licensee. Notwithstanding any damages that the Licensee might incur for any reason whatsoever (including, without limitation, all damages referenced above and all direct or general damages), the entire liability of Devito Codes Ltd and any of its suppliers under any provision of this License and the Licensee's exclusive remedy for all of the foregoing shall be limited to the amount actually paid by the Licensee for the Software. The foregoing limitations, exclusions, and disclaimers shall apply to the maximum extent permitted by applicable law, even if any remedy fails its essential purpose.
10. All rights, title, and interest in and to all data, information, and inventions that result from the Software by the Licensee shall vest in and belong to the Licensee. Licensee agrees to provide Devito Codes with copies of publications referencing the Software and acknowledge Devito Codes in those publications.
11. The Licensee agrees to notify Devito Codes in writing of any known or suspected distribution or use of the Software not in compliance with the License requirements and to enforce the terms of the License.
12. The Licensee agrees to indemnify and hold harmless Devito Codes against any claims or damages arising from the Licensee's use of the Software, exceeding the scope of this License. This includes any third-party claims related to the Software's use.
13. Intellectual Property Rights: All intellectual property rights in the Protected Data/Software and any derivatives thereof shall remain with Devito Codes. The Licensee acknowledges that no ownership rights are transferred by this Agreement.
14. Compliance with Future Technologies: The Licensee commits to continuously ensuring that the use of the Software, including in the context of emerging technologies and computational techniques, adheres to the terms outlined in this License. Devito Codes will provide guidelines for compliance as new technologies develop, facilitating ongoing lawful use of the Software.
15. Notification of Unauthorized Use: The Licensee shall promptly notify Devito Codes of any unauthorized use, copying, or distribution of the Protected Data/Software, or any breach of this License, of which it becomes aware.
16. Where a clause in this License and the Agreement conflict, the clause in the Agreement will take precedence.
 */

#define _POSIX_C_SOURCE 200809L
#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);
#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;

#include "stdlib.h"
#include "math.h"
#include "sys/time.h"
#include "xmmintrin.h"
#include "pmmintrin.h"
#include "cblas.h"
#include "stdio.h"
#include "omp.h"

struct dataobj
{
  void *restrict data;
  int * size;
  unsigned long nbytes;
  unsigned long * npsize;
  unsigned long * dsize;
  int * hsize;
  int * hofs;
  int * oofs;
  void * dmap;
} ;

struct profiler
{
  double section0;
} ;


int Kernel(struct dataobj *restrict u_vec, struct dataobj *restrict v_vec, const int x_M, const int x_m, const int nthreads, struct profiler * timers)
{
  float *restrict u __attribute__ ((aligned (64))) = (float (*)) u_vec->data;
  float *restrict v __attribute__ ((aligned (64))) = (float (*)) v_vec->data;

  /* Flush denormal numbers to zero in hardware */
  _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
  _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);

  START(section0)
  #pragma omp parallel num_threads(nthreads)
  {
    #pragma omp for schedule(static,1)
    for (int x = x_m; x <= x_M; x += 1)
    {
      u[x + 1] = x;
      v[x + 1] = 2*x + 1;
    }
  }
  STOP(section0,timers)

  float s = cblas_sdot(5,u,1,v,1);
  fprintf(stderr,"Out: %f\n",s);

  return 0;
}

And applying our kernel yields the dot product of the two discrete Functions:

\(0(0) + 1(3) + 2(5) = 13\)

#NBVAL_IGNORE_OUTPUT
op2.apply()
Operator `Kernel` ran in 0.01 s
Out: 13.000000
PerformanceSummary([(PerfKey(name='section0', rank=None),
                     PerfEntry(time=0.000261, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])

Inserting a call from a library: Part 2

In this next example, we will use BLAS to implement a 1D diffusion kernel using cblas_sgemv (single-precision general matrix-vector multiplication), and compare it to a stencil-based version.

class Gemv(BlasCall):
    name = 'cblas_sgemv'

Setting up the matrices and the matrix-vector multiplication is a little more complex than previous applications.

grid = Grid(shape=(5,), extent=(1.,))
time = grid.time_dim
dt = time.spacing

x = grid.dimensions[0]
dx = x.spacing

u = TimeFunction(name='u', grid=grid, space_order=2)
v = TimeFunction(name='v', grid=grid, space_order=2)

# Initialise fields
u.data[:, 2:-2] = 1
v.data[:, 2:-2] = 1

# Array temporaries needed in the kernel
# These are typically only generated by the Devito compiler. However, they
# may be necessary or useful for developing custom C Calls and associated
# machinery
tmp0 = Array(name='tmp0', dimensions=grid.dimensions)
tmp1 = Array(name='tmp1', dimensions=grid.dimensions)

# Matrix operator
n = Dimension('n')
m = Dimension('m')
A = Function(name='A', dimensions=(n, m),
             shape=(v.shape[-1], v.shape[-1]),
             space_order=0)

# Fill matrix diagonal
A.data[:] = np.diag(np.ones(v.shape[-1]-1), k=-1) \
    + np.diag(-2*np.ones(v.shape[-1])) \
    + np.diag(np.ones(v.shape[-1]-1), k=1)

The matrix we are using for the 2nd-derivative operator is as follows:

A.data
Data([[-2.,  1.,  0.,  0.,  0.],
      [ 1., -2.,  1.,  0.,  0.],
      [ 0.,  1., -2.,  1.,  0.],
      [ 0.,  0.,  1., -2.,  1.],
      [ 0.,  0.,  0.,  1., -2.]], dtype=float32)

Now we build our Operator. Note the use of temporaries to store the input and output of the matrix-vector multiplication.

op3 = Operator([Eq(tmp0, v),
                # Uses writes=tmp1, as tmp1 is modified in-place.
                # This prevents possible reordering by the Devito compiler, as it now knows the call
                # will modify tmp1.
                Gemv(['CblasRowMajor', 'CblasNoTrans',
                      A.shape[0], A.shape[1],
                      1, A.base, A.shape[0],
                      tmp0.base, 1, 0, tmp1.base, 1], implicit_dims=time,
                     writes=tmp1),
                Eq(v.forward, v + dt*tmp1/dx**2),
                Eq(u.forward, u + dt*u.dx2)])

Printing the C code, we see the addition of the cblas_sgemv call.

#NBVAL_IGNORE_OUTPUT
Code(str(op3.ccode), language='C')
/* NON-COMMERCIAL, NON-EXCLUSIVE, NON-TRANSFERABLE SOFTWARE LICENSE

Copyright (c) 2021-24, Devito Codes Ltd, 1st Floor, Commerce House Raven Road,
London E18 1HB, United Kingdom.

Definitions:

* "Devito Codes": Refers to Devito Codes Ltd., a company incorporated and operating under the laws of England and Wales, and encompasses all its divisions, affiliates, and subsidiaries involved in the development and licensing of software products.
* "Software": Includes all software products developed by Devito Codes, such as Devito PRO, and any other software titles offered by Devito Codes, along with their associated materials, documentation, and software generated by their respective compilers or runtime environments.
* "Licensee": Denotes any individual, company, or entity that has legally acquired a copy or license of the Software under the terms specified in an Agreement with Devito Codes. This term also encompasses all employees, contractors, and agents acting within the scope of their engagement by the Licensee.
* "Agreement": Refers to the legally binding license agreement executed between the Licensee and Devito Codes, which outlines the terms and conditions under which the Licensee is authorized to use the Software. This Agreement may be specific to a single software product or may encompass multiple products offered by Devito Codes.
* "Protected Data/Software": Encompasses all data, software, algorithms, source and object code, documentation, proprietary methodologies, and intellectual property owned or developed by Devito Codes. This includes, without limitation, all versions and iterations of Devito PRO software, materials, and any output generated by Devito Codes' compilers or software, as well as any other proprietary software and tools developed by Devito Codes.
* "Unauthorized Use": Constitutes any action involving the Protected Data/Software that is not expressly permitted under the terms of the Agreement. This includes, but is not limited to, copying, modifying, distributing, selling, leasing, reverse engineering, decompiling, or incorporating the Protected Data/Software into any external system or software, including data retrieval systems, databases, generative AI models, machine learning models, Retrieval Augmented Generation (RAG) applications, or any form of technology or AI system, without the explicit written consent of Devito Codes.
* "Non-Commercial Use" refers to activities not intended for or directed towards commercial advantage or monetary compensation. This includes, but is not limited to, academic research, educational purposes, and personal projects.

Devito Codes provides the Software on the following terms:

1. Scope of License: The Licensee is hereby granted a non-commercial, non-exclusive, non-transferable license without the right to sublicense the Protected Data/Software for the purpose of internal evaluation, testing, and development within the Licensee's organization. This license explicitly excludes any commercial applications, product development for sale, or external consulting services.
2. Prohibition on Unauthorized Use: The Licensee shall not use, adapt, modify, or incorporate the Protected Data/Software into any software or hardware that is not explicitly authorized by this Agreement. This includes prohibitions against the use of the Protected Data/Software in the development, training, or operation of generative AI, machine learning models, RAG applications, or any form of AI that utilizes external data sources for augmentation without the prior written consent of Devito Codes.
3. Distribution and External Use: The Licensee agrees not to distribute, share, lease, or otherwise make the Protected Data/Software available to any third party, including through cloud services, Software as a Service (SaaS), or any form of redistribution mechanism. The Licensee also agrees not to use the Protected Data/Software to provide any service to external parties that involves the processing of data or computation using the Protected Data/Software.
4. The Licensee is prohibited from reverse engineering, decompiling, disassembling, or otherwise attempting to discover the source code of the Protected Data/Software. This prohibition extends to any form of reverse engineering or hacking intended to access or create unauthorized derivatives of the Software.
5. The Software is experimental and licensed "as is." The license of the Software does not include technical support.
6. The Licensee agrees that it will not license or sell the Software or any other Software, information, or data that incorporate any part of the Software, including derivative works ("Derivatives"), to any other parties.
7. Devito Codes will consider all requests for a commercial license but shall be under no obligation to grant such a license
8. Licensee agrees that any person within the Licensee utilizing the Software will be advised of, and is subject to, the conditions in the License.
9. THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR NON-INFRINGEMENT. Devito Codes Ltd expressly disclaims all warranties, whether express, implied, statutory, or otherwise, with respect to the software provided hereunder, including any warranty of reliability, accuracy, or the absence of errors or defects, and the implied warranty arising from course of performance, course of dealing, or usage of trade. Devito Codes Ltd shall not be liable for any direct, indirect, incidental, special, consequential, or exemplary damages, including but not limited to, damages for loss of profits, goodwill, use, data, or other intangible losses (even if Devito Codes Ltd has been advised of the possibility of such damages), resulting from the use or the inability to use the Software. The Licensee acknowledges and agrees that the entire risk arising out of the use or performance of the Software remains with the Licensee. Notwithstanding any damages that the Licensee might incur for any reason whatsoever (including, without limitation, all damages referenced above and all direct or general damages), the entire liability of Devito Codes Ltd and any of its suppliers under any provision of this License and the Licensee's exclusive remedy for all of the foregoing shall be limited to the amount actually paid by the Licensee for the Software. The foregoing limitations, exclusions, and disclaimers shall apply to the maximum extent permitted by applicable law, even if any remedy fails its essential purpose.
10. All rights, title, and interest in and to all data, information, and inventions that result from the Software by the Licensee shall vest in and belong to the Licensee. Licensee agrees to provide Devito Codes with copies of publications referencing the Software and acknowledge Devito Codes in those publications.
11. The Licensee agrees to notify Devito Codes in writing of any known or suspected distribution or use of the Software not in compliance with the License requirements and to enforce the terms of the License.
12. The Licensee agrees to indemnify and hold harmless Devito Codes against any claims or damages arising from the Licensee's use of the Software, exceeding the scope of this License. This includes any third-party claims related to the Software's use.
13. Intellectual Property Rights: All intellectual property rights in the Protected Data/Software and any derivatives thereof shall remain with Devito Codes. The Licensee acknowledges that no ownership rights are transferred by this Agreement.
14. Compliance with Future Technologies: The Licensee commits to continuously ensuring that the use of the Software, including in the context of emerging technologies and computational techniques, adheres to the terms outlined in this License. Devito Codes will provide guidelines for compliance as new technologies develop, facilitating ongoing lawful use of the Software.
15. Notification of Unauthorized Use: The Licensee shall promptly notify Devito Codes of any unauthorized use, copying, or distribution of the Protected Data/Software, or any breach of this License, of which it becomes aware.
16. Where a clause in this License and the Agreement conflict, the clause in the Agreement will take precedence.
 */

#define _POSIX_C_SOURCE 200809L
#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);
#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;
#define uL0(t,x) u[(t)*x_stride0 + (x)]
#define vL0(t,x) v[(t)*x_stride0 + (x)]

#include "stdlib.h"
#include "math.h"
#include "sys/time.h"
#include "xmmintrin.h"
#include "pmmintrin.h"
#include "cblas.h"
#include "omp.h"

struct dataobj
{
  void *restrict data;
  int * size;
  unsigned long nbytes;
  unsigned long * npsize;
  unsigned long * dsize;
  int * hsize;
  int * hofs;
  int * oofs;
  void * dmap;
} ;

struct profiler
{
  double section0;
  double section1;
} ;


int Kernel(struct dataobj *restrict A_vec, struct dataobj *restrict u_vec, struct dataobj *restrict v_vec, const float dt, const float h_x, const int time_M, const int time_m, const int x_M, const int x_m, const int nthreads, const int x_size, struct profiler * timers)
{
  float *restrict tmp0_vec __attribute__ ((aligned (64)));
  posix_memalign((void**)(&tmp0_vec),64,sizeof(float)*(long)x_size);
  float *restrict tmp1_vec __attribute__ ((aligned (64)));
  posix_memalign((void**)(&tmp1_vec),64,sizeof(float)*(long)x_size);

  float *A __attribute__ ((aligned (64))) = (float *) A_vec->data;
  float *restrict tmp0 __attribute__ ((aligned (64))) = (float (*)) tmp0_vec;
  float *restrict tmp1 __attribute__ ((aligned (64))) = (float (*)) tmp1_vec;
  float *u __attribute__ ((aligned (64))) = (float *) u_vec->data;
  float *v __attribute__ ((aligned (64))) = (float *) v_vec->data;

  const int x_fsz0 = v_vec->size[1];

  const int x_stride0 = x_fsz0;

  /* Flush denormal numbers to zero in hardware */
  _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
  _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);

  float r0 = 1.0F/(h_x*h_x);

  for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
  {
    START(section0)
    #pragma omp parallel num_threads(nthreads)
    {
      #pragma omp for schedule(static,1)
      for (int x = x_m; x <= x_M; x += 1)
      {
        tmp0[x] = vL0(t0, x + 2);
      }
    }
    STOP(section0,timers)

    cblas_sgemv(CblasRowMajor,CblasNoTrans,5,5,1,A,5,tmp0,1,0,tmp1,1);

    START(section1)
    #pragma omp parallel num_threads(nthreads)
    {
      #pragma omp for schedule(dynamic,1)
      for (int x = x_m; x <= x_M; x += 1)
      {
        vL0(t1, x + 2) = dt*r0*tmp1[x] + vL0(t0, x + 2);
        uL0(t1, x + 2) = dt*(r0*uL0(t0, x + 1) - 2.0F*r0*uL0(t0, x + 2) + r0*uL0(t0, x + 3)) + uL0(t0, x + 2);
      }
    }
    STOP(section1,timers)
  }

  free(tmp0_vec);
  free(tmp1_vec);

  return 0;
}

Running the kernel, we see that the results of both approaches to diffusion are the same.

#NBVAL_SKIP
op3.apply(time_M=5, dt=0.01)
print(u.data[0])
print(v.data[0])
Operator `Kernel` ran in 0.01 s
[0.09589221 0.2207028  0.2906519  0.2207028  0.09589221]
[0.09589221 0.2207028  0.29065192 0.2207028  0.09589221]

Up-down separation with FFTW

As a final example, we will demonstrate inserting calls to FFTW3 to perform up-down separation on snapshots for the calculation of imaging conditions. This will implement the approach detailed by Wang et al. 2016 (https://doi.org/10.1190/geo2015-0456.1), which explicitly separates the wavefield into two separate fields.

Note that if one only aims to obtain a cross-correlation imaging condition, the approach of Liu et al. 2011 (https://doi.org/10.1190/1.3533914) is notably more efficient, as it does not require propagation of a complex wavefield.

The approach of Wang et al. 2016 is as follows: by propagating a complex source time series obtained via the Hilbert-transform of the true source time series, the wavefield can be separated into upgoing and downgoing components by taking a Fourier transform along the vertical axis, zeroing the transformed field over either the positive or negative wavenumbers, and then inverse transforming.

We use a callback to perform this operation on the wavefield before the separated wavefields are saved as snapshots.

First we initialize the model:

#NBVAL_IGNORE_OUTPUT
so = 8
nbl = 30
model = demo_model(preset='layers-isotropic', nlayers=5,
                   shape=(301, 301), spacing=(10., 10.),
                   space_order=so, nbl=nbl)
Operator `initdamp` ran in 0.01 s

The complex source is then set up:

#NBVAL_IGNORE_OUTPUT
t0, tn = -100., 1200.
dt = model.critical_dt  # Time step from model grid spacing

time_range = TimeAxis(start=t0, stop=tn, step=dt)

f0 = 0.020  # Source peak frequency is 20Hz (0.020 kHz)

# NOTE: For the approach of Liu et al. 2011, one would not need to use a complex dtype
# here or take a Hilbert transform of the source
src = RickerSource(name='src', grid=model.grid, f0=f0,
                   npoint=1, time_range=time_range, dtype=np.complex64)

src.data[:] = hilbert(np.real(src.data), axis=0)

# Plot the sources against one another
plt.plot(time_range.time_values, np.real(src.data), label='real')
plt.plot(time_range.time_values, np.imag(src.data), label='imag')
plt.legend()
plt.xlabel("Time (ms)")
plt.ylabel("Amplitude")
plt.show()

# First, position source centrally in all dimensions, then set depth
src.coordinates.data[0, :] = np.array(model.domain_size) * .5
src.coordinates.data[0, -1] = 20.  # Depth is 20m

We can then define our TimeFunction.

# NOTE: For the approach of Liu et al. 2011, one would not need to use a complex dtype
# here, as the source and wavefield would be real. However, one would want to either have
# a cast to complex or use the `r2c` version of the plan in the callback.
u = TimeFunction(name='u', grid=model.grid, time_order=2, space_order=so, dtype=np.complex64)

pde = model.m * u.dt2 - u.laplace + model.damp * u.dt
stencil = Eq(u.forward, solve(pde, u.forward))

src_term = src.inject(field=u.forward, expr=src*dt**2/model.m)

We then construct our C callback to perform the separation, along with the snapshotting. The example callback is located in the header file updown.h which consists of the following:

#NBVAL_IGNORE_OUTPUT
contents = !cat updown.h
Code('\n'.join(contents), language='C')
#define tmpuAccess(t,x,y) tmpu[(t)*x_stride + (x)*y_stride + (y)]
#define tmpdAccess(t,x,y) tmpd[(t)*x_stride + (x)*y_stride + (y)]

#include "complex.h"
// NOTE: FFTW3 can be compiled to run with OpenMP (the threaded FFT API
// differs from that illustrated in this example however)
#include "fftw3.h"


void upDown(float _Complex *u, float _Complex *tmpu, float _Complex *tmpd,
  int x_size, int y_size, int t)
{
  const int y_stride = y_size;
  const int x_stride = x_size*y_stride;

  int t0 = t % 3;

  int n[1] = {y_size};

  // NOTE: In practice one may not want to set up and tear down the FFT
  // plan at every timestep. Using additional callbacks to generate
  // the plan at the start of the kernel and tear it down at the end
  // would be more efficient.
  fftwf_plan fft_u = fftwf_plan_many_dft(1, n, x_size,
    &u[t0*x_stride], n, 1, y_stride,
    tmpu, n, 1, y_stride,
    FFTW_FORWARD, FFTW_ESTIMATE);

  fftwf_plan fft_d = fftwf_plan_many_dft(1, n, x_size,
    &u[t0*x_stride], n, 1, y_stride,
    tmpd, n, 1, y_stride,
    FFTW_FORWARD, FFTW_ESTIMATE);

  fftwf_plan ifft_u = fftwf_plan_many_dft(1, n, x_size,
    tmpu, n, 1, y_stride,
    tmpu, n, 1, y_stride,
    FFTW_BACKWARD, FFTW_ESTIMATE);

  fftwf_plan ifft_d = fftwf_plan_many_dft(1, n, x_size,
    tmpd, n, 1, y_stride,
    tmpd, n, 1, y_stride,
    FFTW_BACKWARD, FFTW_ESTIMATE);

  fftwf_execute(fft_u);
  fftwf_execute(fft_d);

  // NOTE: In practice, this loop should be parallelised with OpenMP. One could
  // also FFT u only into tmpu, and then copy the required entries into tmpd.
  for (int x = 0; x <= x_size - 1 ; x += 1)
  {
    for (int y = 0; y <= y_size/2; y += 1)
    {
      tmpdAccess(0, x, y) = 0.;
    }
    for (int y = y_size/2; y <= y_size - 1; y += 1)
    {
      tmpuAccess(0, x, y) = 0.;
    }
  }

  fftwf_execute(ifft_u);
  fftwf_execute(ifft_d);

  fftwf_destroy_plan(fft_u);
  fftwf_destroy_plan(fft_d);
  fftwf_destroy_plan(ifft_u);
  fftwf_destroy_plan(ifft_d);
}

It is generally advisable to write a header for inserting complex C calls into Devito operators, rather than inserting a large number of smaller C calls into the kernel.

We then set up the callback using the following:

# Expects fftw to be installed an in a standard path. When installed from apt (libdfftw3-dev) it should automatically
# be found. If installed elsewhere, set fftwdir to the location of your FFTW3 installation
# Assume installed via apt
fftwdir = "/usr"

class FFTWCall(CCall):
    """Abstract class for C calls using FFTW3"""
    # Search the directory containing this notebook plus the
    # FFTW dir. Note we are using the single-precision version of
    # FFTW3, this being fftw3f.
    header_dirs = (os.path.abspath(""), f'{fftwdir}/include/')
    libs = ('fftw3f',)
    lib_dirs = (f'{fftwdir}/lib',)


class UpDownSnapshot(FFTWCall):
    """Save snapshots of the upward- and downward-propagating wavefields"""
    name = 'upDown'
    header = 'updown.h'


# Snapshotting is set up broadly as per usual
nsnaps = 30
factor = int(np.ceil(src.nt/nsnaps))
time = model.grid.time_dim

time_subsampled = ConditionalDimension('t_sub', parent=time, factor=factor)

# Upgoing and downgoing wavefield snapshots
usave_u = TimeFunction(name='usaveu', grid=model.grid, time_order=0, space_order=so,
                       save=nsnaps, time_dim=time_subsampled, dtype=np.complex64,
                       layers=NoLayers)
usave_d = TimeFunction(name='usaved', grid=model.grid, time_order=0, space_order=so,
                       save=nsnaps, time_dim=time_subsampled, dtype=np.complex64,
                       layers=NoLayers)

# Temporaries as destination for upward and downward portions of
# Fourier-transformed wavefield
tmp_u = TimeFunction(name='tmpu', grid=model.grid, space_order=so, time_order=0, dtype=np.complex64)
tmp_d = TimeFunction(name='tmpd', grid=model.grid, space_order=so, time_order=0, dtype=np.complex64)

# Note that the C call is separate from the snapshotting Eqs. This is to ensure that the
# Devito compiler generates the correct snapshotting code.
# implicit_dims=time_subsampled ensures the callback is nested inside the same conditional
save = [UpDownSnapshot([u.base, tmp_u.base, tmp_d.base,
                        u.symbolic_shape[1], u.symbolic_shape[2], model.grid.time_dim], 
                       implicit_dims=time_subsampled),
        Eq(usave_u, tmp_u, implicit_dims=time_subsampled),
        Eq(usave_d, tmp_d, implicit_dims=time_subsampled)]

We then assemble our operator:

op_split = Operator([stencil] + src_term + save)

The operator is run:

#NBVAL_IGNORE_OUTPUT
op_split(dt=dt)
Operator `Kernel` ran in 0.22 s
PerformanceSummary([(PerfKey(name='section0', rank=None),
                     PerfEntry(time=0.12772799999999995, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
                    (PerfKey(name='section1', rank=None),
                     PerfEntry(time=0.006003000000000002, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
                    (PerfKey(name='section2', rank=None),
                     PerfEntry(time=0.004634, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])

We are then able to plot the snapshots of both the upgoing and downgoing components of the wavefield.

#NBVAL_IGNORE_OUTPUT
def plot_field(field, model, nbl, title):
    clip = 0.5  # Clip value
    vmax = clip*np.amax(np.abs(field[-1]))
    vmin = -vmax

    plt.imshow(model.vp.data[nbl:-nbl, nbl:-nbl].T, cmap='Greys')
    plt.imshow(field[-1, nbl:-nbl, nbl:-nbl].T, cmap='seismic',
               vmin=vmin, vmax=vmax, alpha=0.7)
    plt.title(title)
    plt.show()

plot_field(np.real(usave_u.data), model, nbl, "Upgoing")

plot_field(np.real(usave_d.data), model, nbl, "Downgoing")

from devito import norm

assert np.isclose(norm(usave_d), 190828.75, rtol=1e-3)
assert np.isclose(norm(usave_u), 85467.01, rtol=1e-3)
Back to top