Repository: hpc4cmb/toast (language: python)

Issue:
Parallel conviqt tests hung during figure output.



Discussion / hints:

Fix patch:
diff --git a/src/toast/mpi.py b/src/toast/mpi.py
index de15e07d5..989e8bc4b 100644
--- a/src/toast/mpi.py
+++ b/src/toast/mpi.py
@@ -6,6 +6,10 @@
 
 use_mpi = None
 MPI = None
+# traitlets require the MPI communicator type to be an actual class,
+# even if None is a valid value
+class MPI_Comm:
+    pass
 
 if use_mpi is None:
     # See if the user has explicitly disabled MPI.
@@ -26,6 +30,7 @@
                 import mpi4py.MPI as MPI
 
                 use_mpi = True
+                MPI_Comm = MPI.Comm
             except:
                 # There could be many possible exceptions raised...
                 from ._libtoast import Logger
diff --git a/src/toast/ops/CMakeLists.txt b/src/toast/ops/CMakeLists.txt
index 1abfc0df6..d5a90c1c3 100644
--- a/src/toast/ops/CMakeLists.txt
+++ b/src/toast/ops/CMakeLists.txt
@@ -31,5 +31,6 @@ install(FILES
     madam.py
     madam_utils.py
     conviqt.py
+    totalconvolve.py
     DESTINATION ${PYTHON_SITE}/toast/ops
 )
diff --git a/src/toast/ops/__init__.py b/src/toast/ops/__init__.py
index 2c975bc04..13f2ca6c1 100644
--- a/src/toast/ops/__init__.py
+++ b/src/toast/ops/__init__.py
@@ -57,3 +57,4 @@
 from .madam import Madam
 
 from .conviqt import SimConviqt, SimWeightedConviqt
+from .totalconvolve import SimTotalconvolve
diff --git a/src/toast/ops/conviqt.py b/src/toast/ops/conviqt.py
index 1dc17efd9..bb9e20fe3 100644
--- a/src/toast/ops/conviqt.py
+++ b/src/toast/ops/conviqt.py
@@ -9,7 +9,7 @@
 import numpy as np
 import traitlets
 
-from ..mpi import MPI, use_mpi, Comm
+from ..mpi import MPI, MPI_Comm, use_mpi, Comm
 from .operator import Operator
 from .. import qarray as qa
 from ..timing import function_timer
@@ -41,7 +41,7 @@ class SimConviqt(Operator):
     API = Int(0, help="Internal interface version for this operator")
 
     comm = Instance(
-        klass=MPI.Comm,
+        klass=MPI_Comm,
         allow_none=True,
         help="MPI communicator to use for the convolution. libConviqt does "
         "not work without MPI.",
diff --git a/src/toast/ops/totalconvolve.py b/src/toast/ops/totalconvolve.py
new file mode 100644
index 000000000..e1fcde956
--- /dev/null
+++ b/src/toast/ops/totalconvolve.py
@@ -0,0 +1,822 @@
+# Copyright (c) 2015-2020 by the parties listed in the AUTHORS file.
+# All rights reserved.  Use of this source code is governed by
+# a BSD-style license that can be found in the LICENSE file.
+
+import os
+import warnings
+
+from astropy import units as u
+import numpy as np
+import traitlets
+import healpy as hp
+
+from ..mpi import MPI, MPI_Comm, use_mpi, Comm
+from .operator import Operator
+from .. import qarray as qa
+from ..timing import function_timer
+from ..traits import trait_docs, Int, Float, Unicode, Bool, Dict, Quantity, Instance
+from ..utils import Logger, Environment, Timer, GlobalTimers, dtype_to_aligned
+
+
+totalconvolve = None
+
+try:
+    import ducc0.totalconvolve as totalconvolve
+except ImportError:
+    totalconvolve = None
+
+
+def available():
+    """(bool): True if ducc0.totalconvolve is found in the library search path."""
+    global totalconvolve
+    return totalconvolve is not None
+
+
+@trait_docs
+class SimTotalconvolve(Operator):
+    """Operator which uses ducc0.totalconvolve to generate beam-convolved timestreams."""
+
+    # Class traits
+
+    API = Int(0, help="Internal interface version for this operator")
+
+    comm = Instance(
+        klass=MPI_Comm,
+        allow_none=True,
+        help="MPI communicator to use for the convolution.",
+    )
+
+    detector_pointing = Instance(
+        klass=Operator,
+        allow_none=True,
+        help="Operator that translates boresight pointing into detector frame",
+    )
+
+    det_flags = Unicode(
+        None, allow_none=True, help="Observation detdata key for flags to use"
+    )
+
+    det_flag_mask = Int(0, help="Bit mask value for optional detector flagging")
+
+    shared_flags = Unicode(
+        None, allow_none=True, help="Observation shared key for telescope flags to use"
+    )
+
+    shared_flag_mask = Int(0, help="Bit mask value for optional shared flagging")
+
+    view = Unicode(
+        None, allow_none=True, help="Use this view of the data in all observations"
+    )
+
+    det_data = Unicode(
+        "signal",
+        allow_none=False,
+        help="Observation detdata key for accumulating convolved timestreams",
+    )
+
+    calibrate = Bool(
+        True,
+        allow_none=False,
+        help="Calibrate intensity to 1.0, rather than (1 + epsilon) / 2. "
+        "Calibrate has no effect if the beam is found to be normalized rather than "
+        "scaled with the leakage factor.",
+    )
+
+    dxx = Bool(
+        True,
+        allow_none=False,
+        help="The beam frame is either Dxx or Pxx. Pxx includes the rotation to "
+        "polarization sensitive basis, Dxx does not. When Dxx=True, detector "
+        "orientation from attitude quaternions is corrected for the polarization "
+        "angle.",
+    )
+
+    pol = Bool(
+        True,
+        allow_none=False,
+        help="Toggle simulated signal polarization",
+    )
+
+    mc = Int(
+        None,
+        allow_none=True,
+        help="Monte Carlo index used in synthesizing the input file names.",
+    )
+
+    @traitlets.validate("mc")
+    def _check_mc(self, proposal):
+        check = proposal["value"]
+        if check is not None and check < 0:
+            raise traitlets.TraitError("MC index cannot be negative")
+        return check
+
+    beammmax = Int(
+        -1,
+        allow_none=False,
+        help="Beam maximum m.  Actual resolution in the Healpix FITS file may differ. "
+        "If not set, will use the maximum expansion order from file.",
+    )
+
+    oversampling_factor = Float(
+        1.8,
+        allow_none=False,
+        help="Oversampling factor for total convolution (useful range is 1.5-2.0)",
+    )
+
+    epsilon = Float(
+        1e-5,
+        allow_none=False,
+        help="Relative accuracy of the interpolation step",
+    )
+
+    lmax = Int(
+        -1,
+        allow_none=False,
+        help="Maximum ell (and m).  Actual resolution in the Healpix FITS file may "
+        "differ.  If not set, will use the maximum expansion order from file.",
+    )
+
+    verbosity = Int(
+        0,
+        allow_none=False,
+        help="",
+    )
+
+    normalize_beam = Bool(
+        False,
+        allow_none=False,
+        help="Normalize beam to have unit response to temperature monopole.",
+    )
+
+    remove_dipole = Bool(
+        False,
+        allow_none=False,
+        help="Suppress the temperature dipole in sky_file.",
+    )
+
+    remove_monopole = Bool(
+        False,
+        allow_none=False,
+        help="Suppress the temperature monopole in sky_file.",
+    )
+
+    apply_flags = Bool(
+        False,
+        allow_none=False,
+        help="Only synthesize signal for unflagged samples.",
+    )
+
+    fwhm = Quantity(
+        4.0 * u.arcmin,
+        allow_none=False,
+        help="Width of a symmetric gaussian beam already present in the skyfile "
+        "(will be deconvolved away).",
+    )
+
+    sky_file_dict = Dict(
+        None,
+        allow_none=True,
+        help="Dictionary of files containing the sky a_lm expansions. An entry for "
+        "each detector name must be present. If provided, supersedes `sky_file`.",
+    )
+
+    sky_file = Unicode(
+        None,
+        allow_none=True,
+        help="File containing the sky a_lm expansion.  Tag {detector} will be "
+        "replaced with the detector name",
+    )
+
+    beam_file_dict = Dict(
+        None,
+        allow_none=True,
+        help="Dictionary of files containing the beam a_lm expansions. An entry for "
+        "each detector name must be present. If provided, supersedes `beam_file`.",
+    )
+
+    beam_file = Unicode(
+        None,
+        allow_none=True,
+        help="File containing the beam a_lm expansion.  Tag {detector} will be "
+        "replaced with the detector name.",
+    )
+
+    @traitlets.validate("shared_flag_mask")
+    def _check_shared_flag_mask(self, proposal):
+        check = proposal["value"]
+        if check < 0:
+            raise traitlets.TraitError("Shared flag mask should be a positive integer")
+        return check
+
+    @traitlets.validate("det_flag_mask")
+    def _check_det_flag_mask(self, proposal):
+        check = proposal["value"]
+        if check < 0:
+            raise traitlets.TraitError("Det flag mask should be a positive integer")
+        return check
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        return
+
+    @property
+    def available(self):
+        """Return True if ducc0.totalconvolve is found in the library search path."""
+        return totalconvolve is not None
+
+    hwp_angle = Unicode(
+        None, allow_none=True, help="Observation shared key for HWP angle"
+    )
+
+    @function_timer
+    def _exec(self, data, detectors=None, **kwargs):
+        if not self.available:
+            raise RuntimeError("ducc0.totalconvolve is not available")
+
+        if self.detector_pointing is None:
+            raise RuntimeError("detector_pointing cannot be None.")
+
+        log = Logger.get()
+
+        timer = Timer()
+        timer.start()
+
+        env = Environment.get()
+        nthreads = env.max_threads()
+
+        verbose = self.verbosity > 0
+        if use_mpi:
+            self.comm.barrier()
+            verbose = verbose and self.comm.rank == 0
+            if self.comm.size > 1 and self.comm.rank == 0:
+                log.warning(
+                    "communicator size>1: totalconvolve will work, "
+                    "but will waste memory. To be fixed in future releases."
+                )
+
+        all_detectors = self._get_all_detectors(data, detectors)
+
+        for det in all_detectors:
+            # Expand detector pointing
+            self.detector_pointing.apply(data, detectors=[det])
+
+            if det in self.sky_file_dict:
+                sky_file = self.sky_file_dict[det]
+            else:
+                sky_file = self.sky_file.format(detector=det, mc=self.mc)
+
+            if det in self.beam_file_dict:
+                beam_file = self.beam_file_dict[det]
+            else:
+                beam_file = self.beam_file.format(detector=det, mc=self.mc)
+
+            lmax, mmax = self.get_lmmax(sky_file, beam_file)
+            sky = self.get_sky(sky_file, lmax, det, verbose)
+            beam = self.get_beam(beam_file, lmax, mmax, det, verbose)
+
+            theta, phi, psi, psi_pol = self.get_pointing(data, det, verbose)
+            pnt = self.get_buffer(theta, phi, psi, det, verbose)
+            del theta, phi, psi
+            if self.hwp_angle is None:
+                psi_pol = None
+
+            convolved_data = self.convolve(
+                sky, beam, lmax, mmax, pnt, psi_pol, det, nthreads, verbose
+            )
+            del psi_pol
+
+            self.calibrate_signal(data, det, beam, convolved_data, verbose)
+            self.save(data, det, convolved_data, verbose)
+
+            del pnt, beam, sky
+
+            if verbose:
+                timer.report_clear(f"totalconvolve process detector {det}")
+
+        return
+
+    def _get_all_detectors(self, data, detectors):
+        """Assemble a list of detectors across all processes and
+        observations in `self._comm`.
+        """
+        my_dets = set()
+        for obs in data.obs:
+            # Get the detectors we are using for this observation
+            obs_dets = obs.select_local_detectors(detectors)
+            for det in obs_dets:
+                my_dets.add(det)
+            # Make sure detector data output exists
+            obs.detdata.ensure(self.det_data, detectors=detectors)
+        if use_mpi:
+            all_dets = self.comm.gather(my_dets, root=0)
+            if self.comm.rank == 0:
+                for some_dets in all_dets:
+                    my_dets.update(some_dets)
+                my_dets = sorted(my_dets)
+            all_dets = self.comm.bcast(my_dets, root=0)
+        else:
+            all_dets = my_dets
+        return all_dets
+
+    def _get_psi_pol(self, focalplane, det):
+        """Parse polarization angle in radians from the focalplane
+        dictionary.  The angle is relative to the Pxx basis.
+        """
+        if det not in focalplane:
+            raise RuntimeError(f"focalplane does not include {det}")
+        props = focalplane[det]
+        if "psi_pol" in props.colnames:
+            psi_pol = props["pol_angle"].to_value(u.radian)
+        elif "pol_angle" in props.colnames:
+            warnings.warn(
+                "Use psi_pol and psi_uv rather than pol_angle", DeprecationWarning
+            )
+            psi_pol = props["pol_angle"].to_value(u.radian)
+        else:
+            raise RuntimeError(f"focalplane[{det}] does not include psi_pol")
+        return psi_pol
+
+    def _get_psi_uv(self, focalplane, det):
+        """Parse Pxx basis angle in radians from the focalplane
+        dictionary.  The angle is measured from Dxx to Pxx basis.
+        """
+        if det not in focalplane:
+            raise RuntimeError(f"focalplane does not include {det}")
+        props = focalplane[det]
+        if "psi_uv_deg" in props.colnames:
+            psi_uv = props["psi_uv"].to_value(u.radian)
+        else:
+            raise RuntimeError(f"focalplane[{det}] does not include psi_uv")
+        return psipol
+
+    def _get_epsilon(self, focalplane, det):
+        """Parse polarization leakage (epsilon) from the focalplane
+        object or dictionary.
+        """
+        if det not in focalplane:
+            raise RuntimeError(f"focalplane does not include {det}")
+        props = focalplane[det]
+        if "pol_leakage" in props.colnames:
+            epsilon = focalplane[det]["pol_leakage"]
+        else:
+            # Assume zero polarization leakage
+            epsilon = 0
+        return epsilon
+
+    def get_lmmax(self, skyfile, beamfile):
+        """Determine the actual lmax and beammmax to use for the convolution
+        from class parameters and values in the files.
+        """
+        ncomp = 3 if self.pol else 1
+        slmax, blmax, bmmax = -1, -1, -1
+        for i in range(ncomp):
+            # for sky and beam respectively, lmax is the max of all components
+            alm_tmp, mmax_tmp = hp.fitsfunc.read_alm(
+                skyfile, hdu=i + 1, return_mmax=True
+            )
+            lmax_tmp = hp.sphtfunc.Alm.getlmax(alm_tmp.shape[0], mmax_tmp)
+            slmax = max(slmax, lmax_tmp)
+            alm_tmp, mmax_tmp = hp.fitsfunc.read_alm(
+                beamfile, hdu=i + 1, return_mmax=True
+            )
+            lmax_tmp = hp.sphtfunc.Alm.getlmax(alm_tmp.shape[0], mmax_tmp)
+            blmax = max(blmax, lmax_tmp)
+            # for the beam, determine also the largest mmax present
+            bmmax = max(bmmax, mmax_tmp)
+        # no need to go higher than the lower of the lmax from sky and beam
+        lmax_out = min(slmax, blmax)
+        mmax_out = bmmax
+        # if parameters are lower than the detected values, reduce even further
+        if self.lmax != -1:
+            lmax_out = min(lmax_out, self.lmax)
+        if self.beammmax != -1:
+            mmax_out = min(mmax_out, self.beammmax)
+        return lmax_out, mmax_out
+
+    def load_alm(self, file, lmax, mmax):
+        def read_comp(file, comp, out):
+            almX, mmaxX = hp.fitsfunc.read_alm(file, hdu=comp + 1, return_mmax=True)
+            lmaxX = hp.sphtfunc.Alm.getlmax(almX.shape[0], mmaxX)
+
+            ofs1, ofs2 = 0, 0
+            mylmax = min(lmax, lmaxX)
+            for m in range(0, min(mmax, mmaxX) + 1):
+                out[comp, ofs1 : ofs1 + mylmax - m + 1] = almX[
+                    ofs2 : ofs2 + mylmax - m + 1
+                ]
+                ofs1 += lmax - m + 1
+                ofs2 += lmaxX - m + 1
+
+        ncomp = 3 if self.pol else 1
+        res = np.zeros(
+            (ncomp, hp.sphtfunc.Alm.getsize(lmax, mmax)), dtype=np.complex128
+        )
+        for i in range(ncomp):
+            read_comp(file, i, res)
+        return res
+
+    def get_sky(self, skyfile, lmax, det, verbose):
+        timer = Timer()
+        timer.start()
+        sky = self.load_alm(skyfile, lmax, lmax)
+        fwhm = self.fwhm.to_value(u.radian)
+        if fwhm != 0:
+            gauss = hp.sphtfunc.gauss_beam(fwhm, lmax, pol=True)
+            for i in range(sky.shape[0]):
+                sky[i] = hp.sphtfunc.almxfl(
+                    sky[i], 1.0 / gauss[:, i], mmax=lmax, inplace=True
+                )
+        if self.remove_monopole:
+            sky[0, 0] = 0
+        if self.remove_dipole:
+            sky[0, 1] = 0
+            sky[0, lmax + 1] = 0
+        if verbose:
+            timer.report_clear(f"initialize sky for detector {det}")
+        return sky
+
+    def get_beam(self, beamfile, lmax, mmax, det, verbose):
+        timer = Timer()
+        timer.start()
+        beam = self.load_alm(beamfile, lmax, mmax)
+        if self.normalize_beam:
+            beam *= 1.0 / (2 * np.sqrt(np.pi) * beam[0, 0])
+        if verbose:
+            timer.report_clear(f"initialize beam for detector {det}")
+        return beam
+
+    def get_pointing(self, data, det, verbose):
+        """Return the detector pointing as ZYZ Euler angles without the
+        polarization sensitive angle.  These angles are to be compatible
+        with Pxx or Dxx frame beam products
+        """
+        # We need the three pointing angles to describe the
+        # pointing.  local_pointing() returns the attitude quaternions.
+        nullquat = np.array([0, 0, 0, 1], dtype=np.float64)
+        timer = Timer()
+        timer.start()
+        all_theta, all_phi, all_psi, all_psi_pol = [], [], [], []
+        for obs in data.obs:
+            if det not in obs.local_detectors:
+                continue
+            focalplane = obs.telescope.focalplane
+            # Loop over views
+            views = obs.view[self.view]
+            for view in range(len(views)):
+                # Get the flags if needed
+                flags = None
+                if self.apply_flags:
+                    if self.shared_flags is not None:
+                        flags = np.array(views.shared[self.shared_flags][view])
+                        flags &= self.shared_flag_mask
+                    if self.det_flags is not None:
+                        detflags = np.array(views.detdata[self.det_flags][view][det])
+                        detflags &= self.det_flag_mask
+                        if flags is not None:
+                            flags |= detflags
+                        else:
+                            flags = detflags
+
+                # Timestream of detector quaternions
+                quats = views.detdata[self.detector_pointing.quats][view][det]
+                if verbose:
+                    timer.report_clear(f"get detector pointing for {det}")
+
+                if flags is not None:
+                    quats = quats.copy()
+                    quats[flags != 0] = nullquat
+                    if verbose:
+                        timer.report_clear(f"initialize flags for detector {det}")
+
+                theta, phi, psi = qa.to_angles(quats)
+                # Polarization angle in the Pxx basis
+                psi_pol = self._get_psi_pol(focalplane, det)
+                if self.dxx:
+                    # Add angle between Dxx and Pxx
+                    psi_pol += self._get_psi_uv(focalplane, det)
+                psi -= psi_pol
+                psi_pol = np.ones(psi.size) * psi_pol
+                if self.hwp_angle is not None:
+                    hwp_angle = views.shared[self.hwp_angle][view]
+                    psi_pol += 2 * hwp_angle
+                all_theta.append(theta)
+                all_phi.append(phi)
+                all_psi.append(psi)
+                all_psi_pol.append(psi_pol)
+
+        if len(all_theta) > 0:
+            all_theta = np.hstack(all_theta)
+            all_phi = np.hstack(all_phi)
+            all_psi = np.hstack(all_psi)
+            all_psi_pol = np.hstack(all_psi_pol)
+
+        if verbose:
+            timer.report_clear(f"compute pointing angles for detector {det}")
+        return all_theta, all_phi, all_psi, all_psi_pol
+
+    def get_buffer(self, theta, phi, psi, det, verbose):
+        """Pack the pointing into the pointing array"""
+        timer = Timer()
+        timer.start()
+        pnt = np.empty((len(theta), 3))
+        pnt[:, 0] = theta
+        pnt[:, 1] = phi
+        pnt[:, 2] = psi
+        pnt[:, 2] += np.pi  # FIXME: not clear yet why this is necessary
+        if verbose:
+            timer.report_clear(f"pack input array for detector {det}")
+        return pnt
+
+    # simple approach when there is only one task
+    def conv_and_interpol_serial(
+        self, skycomp, beamcomp, lmax, mmax, pnt, nthreads, t_conv, t_inter
+    ):
+        t_conv.start()
+        plan = totalconvolve.ConvolverPlan(
+            lmax=lmax,
+            kmax=mmax,
+            sigma=self.oversampling_factor,
+            epsilon=self.epsilon,
+            nthreads=nthreads,
+        )
+        cube = np.empty((plan.Npsi(), plan.Ntheta(), plan.Nphi()), dtype=np.float64)
+        cube[()] = 0
+
+        # convolution part
+        for icomp in range(skycomp.shape[0]):
+            plan.getPlane(skycomp[icomp, :], beamcomp[icomp, :], 0, cube[0:1])
+            for mbeam in range(1, mmax + 1):
+                plan.getPlane(
+                    skycomp[icomp, :],
+                    beamcomp[icomp, :],
+                    mbeam,
+                    cube[2 * mbeam - 1 : 2 * mbeam + 1],
+                )
+
+        plan.prepPsi(cube)
+        t_conv.stop()
+        t_inter.start()
+        res = np.empty(pnt.shape[0], dtype=np.float64)
+        plan.interpol(cube, 0, 0, pnt[:, 0], pnt[:, 1], pnt[:, 2], res)
+        t_inter.stop()
+        return res
+
+    # MPI version storing the full data cube at every MPI task (wasteful)
+    def conv_and_interpol_mpi(
+        self, skycomp, beamcomp, lmax, mmax, pnt, nthreads, t_conv, t_inter
+    ):
+        t_conv.start()
+        plan = totalconvolve.ConvolverPlan(
+            lmax=lmax,
+            kmax=mmax,
+            sigma=self.oversampling_factor,
+            epsilon=self.epsilon,
+            nthreads=nthreads,
+        )
+        myrank, nranks = self.comm.rank, self.comm.size
+        cube = np.empty((plan.Npsi(), plan.Ntheta(), plan.Nphi()), dtype=np.float64)
+        cube[()] = 0
+
+        # convolution part
+        # the work in this nested loop can be distributed over skycomp.shape[0]*(mmax+1) tasks
+        for icomp in range(skycomp.shape[0]):
+            if (icomp * (mmax + 1)) % nranks == myrank:
+                plan.getPlane(skycomp[icomp, :], beamcomp[icomp, :], 0, cube[0:1])
+            for mbeam in range(1, mmax + 1):
+                if (icomp * (mmax + 1) + mbeam) % nranks == myrank:
+                    plan.getPlane(
+                        skycomp[icomp, :],
+                        beamcomp[icomp, :],
+                        mbeam,
+                        cube[2 * mbeam - 1 : 2 * mbeam + 1],
+                    )
+        if nranks > 1:  # broadcast the results
+            for icomp in range(skycomp.shape[0]):
+                self.comm.Bcast(
+                    [cube[0:1], MPI.DOUBLE], root=(icomp * (mmax + 1)) % nranks
+                )
+            for mbeam in range(1, mmax + 1):
+                self.comm.Bcast(
+                    [cube[2 * mbeam - 1 : 2 * mbeam + 1], MPI.DOUBLE],
+                    root=(icomp * (mmax + 1) + mbeam) % nranks,
+                )
+
+        plan.prepPsi(cube)
+        t_conv.stop()
+        t_inter.start()
+        res = np.empty(pnt.shape[0], dtype=np.float64)
+        plan.interpol(cube, 0, 0, pnt[:, 0], pnt[:, 1], pnt[:, 2], res)
+        t_inter.stop()
+        return res
+
+    # MPI version with shared memory tricks, storing the full data cube only
+    # once per node.
+    def conv_and_interpol_mpi_shmem(
+        self, skycomp, beamcomp, lmax, mmax, pnt, nthreads, t_conv, t_inter
+    ):
+        from pshmem import MPIShared
+
+        t_conv.start()
+        plan = totalconvolve.ConvolverPlan(
+            lmax=lmax,
+            kmax=mmax,
+            sigma=self.oversampling_factor,
+            epsilon=self.epsilon,
+            nthreads=nthreads,
+        )
+
+        with MPIShared(
+            (plan.Npsi(), plan.Ntheta(), plan.Nphi()), np.float64, self.comm
+        ) as shm:
+            cube = shm.data
+            # Create a separate communicator on every node.
+            intracomm = self.comm.Split_type(MPI.COMM_TYPE_SHARED)
+            # Create a communicator with all master tasks of the intracomms;
+            # on every other task, intercomm will be MPI.COMM_NULL.
+            color = 0 if intracomm.rank == 0 else MPI.UNDEFINED
+            intercomm = self.comm.Split(color)
+            if intracomm.rank == 0:
+                # We are on the master task of intracomm and all other tasks on
+                # the node will be idle during the next computation step,
+                # so we can hijack all their threads.
+                nodeplan = totalconvolve.ConvolverPlan(
+                    lmax=lmax,
+                    kmax=mmax,
+                    sigma=self.oversampling_factor,
+                    epsilon=self.epsilon,
+                    nthreads=intracomm.size * nthreads,
+                )
+                mynode, nnodes = intercomm.rank, intercomm.size
+                cube[()] = 0.0
+
+                # Convolution part
+                # The skycomp.shape[0]*(mmax+1) work items in this nested loop
+                # are distributed among the nodes in a round-robin fashion.
+                for icomp in range(skycomp.shape[0]):
+                    if (icomp * (mmax + 1)) % nnodes == mynode:
+                        nodeplan.getPlane(
+                            skycomp[icomp, :], beamcomp[icomp, :], 0, cube[0:1]
+                        )
+                    for mbeam in range(1, mmax + 1):
+                        if (icomp * (mmax + 1) + mbeam) % nnodes == mynode:
+                            nodeplan.getPlane(
+                                skycomp[icomp, :],
+                                beamcomp[icomp, :],
+                                mbeam,
+                                cube[2 * mbeam - 1 : 2 * mbeam + 1],
+                            )
+                if nnodes > 1:  # results must be broadcast to all nodes
+                    for icomp in range(skycomp.shape[0]):
+                        intercomm.Bcast(
+                            [cube[0:1], MPI.DOUBLE], root=(icomp * (mmax + 1)) % nnodes
+                        )
+                    for mbeam in range(1, mmax + 1):
+                        intercomm.Bcast(
+                            [cube[2 * mbeam - 1 : 2 * mbeam + 1], MPI.DOUBLE],
+                            root=(icomp * (mmax + 1) + mbeam) % nnodes,
+                        )
+                nodeplan.prepPsi(cube)
+                del nodeplan
+
+            t_conv.stop()
+            t_inter.start()
+            # Interpolation part
+            # No fancy communication is necessary here, since every task has
+            # access to the full data cube.
+            res = np.empty(pnt.shape[0], dtype=np.float64)
+            plan.interpol(cube, 0, 0, pnt[:, 0], pnt[:, 1], pnt[:, 2], res)
+            t_inter.stop()
+            del plan
+            del cube
+        return res
+
+    def conv_and_interpol(
+        self, skycomp, beamcomp, lmax, mmax, pnt, nthreads, t_conv, t_inter
+    ):
+        if (not use_mpi) or (self.comm.size == 1):
+            return self.conv_and_interpol_serial(
+                skycomp, beamcomp, lmax, mmax, pnt, nthreads, t_conv, t_inter
+            )
+        else:
+            return self.conv_and_interpol_mpi_shmem(
+                skycomp, beamcomp, lmax, mmax, pnt, nthreads, t_conv, t_inter
+            )
+
+    def convolve(self, sky, beam, lmax, mmax, pnt, psi_pol, det, nthreads, verbose):
+        t_conv = Timer()
+        t_inter = Timer()
+
+        if self.hwp_angle is None:
+            # simply compute TT+EE+BB
+            convolved_data = self.conv_and_interpol(
+                np.array(sky),
+                np.array(beam),
+                lmax,
+                mmax,
+                pnt,
+                nthreads,
+                t_conv,
+                t_inter,
+            )
+        else:
+            # TT
+            convolved_data = self.conv_and_interpol(
+                np.array([sky[0]]),
+                np.array([beam[0]]),
+                lmax,
+                mmax,
+                pnt,
+                nthreads,
+                t_conv,
+                t_inter,
+            )
+            if self.pol:
+                # EE+BB
+                slm = np.array([sky[1], sky[2]])
+                blm = np.array([beam[1], beam[2]])
+                convolved_data += np.cos(4 * psi_pol) * self.conv_and_interpol(
+                    slm, blm, lmax, mmax, pnt, nthreads, t_conv, t_inter
+                ).reshape((-1,))
+                # -EB+BE
+                blm = np.array([-beam[2], beam[1]])
+                convolved_data += np.sin(4 * psi_pol) * self.conv_and_interpol(
+                    slm, blm, lmax, mmax, pnt, nthreads, t_conv, t_inter
+                ).reshape((-1,))
+
+        if verbose:
+            t_conv.report_clear(f"convolve detector {det}")
+            t_inter.report_clear(f"extract convolved data for {det}")
+
+        convolved_data *= 0.5  # FIXME: not sure where this factor comes from
+        return convolved_data
+
+    def calibrate_signal(self, data, det, beam, convolved_data, verbose):
+        """By default, libConviqt results returns a signal that conforms to
+        TOD = (1 + epsilon) / 2 * intensity + (1 - epsilon) / 2 * polarization.
+
+        When calibrate = True, we rescale the TOD to
+        TOD = intensity + (1 - epsilon) / (1 + epsilon) * polarization
+        """
+        if not self.calibrate:  # or beam.normalized():
+            return
+
+        timer = Timer()
+        timer.start()
+        offset = 0
+        for obs in data.obs:
+            if det not in obs.local_detectors:
+                continue
+            focalplane = obs.telescope.focalplane
+            epsilon = self._get_epsilon(focalplane, det)
+            # Make sure detector data output exists
+            obs.detdata.ensure(self.det_data, detectors=[det])
+            # Loop over views
+            views = obs.view[self.view]
+            for view in views.detdata[self.det_data]:
+                nsample = len(view[det])
+                convolved_data[offset : offset + nsample] *= 2 / (1 + epsilon)
+                offset += nsample
+        if verbose:
+            timer.report_clear(f"calibrate detector {det}")
+        return
+
+    def save(self, data, det, convolved_data, verbose):
+        """Store the convolved data."""
+        timer = Timer()
+        timer.start()
+        offset = 0
+        for obs in data.obs:
+            if det not in obs.local_detectors:
+                continue
+            # Loop over views
+            views = obs.view[self.view]
+            for view in views.detdata[self.det_data]:
+                nsample = len(view[det])
+                view[det] += convolved_data[offset : offset + nsample]
+                offset += nsample
+        if verbose:
+            timer.report_clear(f"save detector {det}")
+        return
+
+    def _finalize(self, data, **kwargs):
+        return
+
+    def _requires(self):
+        req = self.detector_pointing.requires()
+        req["meta"].extend([self.noise_model, self.pixel_dist, self.covariance])
+        req["shared"] = [self.boresight]
+        if self.shared_flags is not None:
+            req["shared"].append(self.shared_flags)
+        if self.det_flags is not None:
+            req["detdata"].append(self.det_flags)
+        if self.view is not None:
+            req["intervals"].append(self.view)
+        return req
+
+    def _provides(self):
+        prov = self.detector_pointing.provides()
+        prob["detdata"].append(self.det_data)
+        return prov
+
+    def _accelerators(self):
+        return list()


Test patch:
diff --git a/src/toast/tests/CMakeLists.txt b/src/toast/tests/CMakeLists.txt
index 09e4aade5..e45056ed5 100644
--- a/src/toast/tests/CMakeLists.txt
+++ b/src/toast/tests/CMakeLists.txt
@@ -23,6 +23,7 @@ install(FILES
     ops_sim_tod_noise.py
     ops_sim_tod_dipole.py
     ops_sim_tod_conviqt.py
+    ops_sim_tod_totalconvolve.py
     ops_mapmaker_utils.py
     ops_mapmaker_binning.py
     ops_mapmaker_solve.py
diff --git a/src/toast/tests/_helpers.py b/src/toast/tests/_helpers.py
index 38b3f8f9d..79d07dc3a 100644
--- a/src/toast/tests/_helpers.py
+++ b/src/toast/tests/_helpers.py
@@ -367,24 +367,38 @@ def uniform_chunks(samples, nchunk=100):
     return chunks
 
 
-def create_fake_sky_alm(lmax=128, fwhm=10 * u.degree, pol=True):
-    # Power spectrum
-    if pol:
-        cl = np.ones(4 * (lmax + 1)).reshape([4, -1])
+def create_fake_sky_alm(lmax=128, fwhm=10 * u.degree, pol=True, pointsources=False):
+    if pointsources:
+        nside = 512
+        while nside < lmax:
+            nside *= 2
+        npix = 12 * nside ** 2
+        m = np.zeros(npix)
+        for lon in np.linspace(-180, 180, 6):
+            for lat in np.linspace(-80, 80, 6):
+                m[hp.ang2pix(nside, lon, lat, lonlat=True)] = 1
+        if pol:
+            m = np.vstack([m, m, m])
+        m = hp.smoothing(m, fwhm=fwhm.to_value(u.radian))
+        a_lm = hp.map2alm(m, lmax=lmax)
     else:
-        cl = np.ones(lmax + 1)
-    # Draw a_lm
-    nside = 2
-    while 4 * nside < lmax:
-        nside *= 2
-    _, a_lm = hp.synfast(
-        cl,
-        nside,
-        alm=True,
-        lmax=lmax,
-        fwhm=fwhm.to_value(u.radian),
-        verbose=False,
-    )
+        # Power spectrum
+        if pol:
+            cl = np.ones(4 * (lmax + 1)).reshape([4, -1])
+        else:
+            cl = np.ones(lmax + 1)
+        # Draw a_lm
+        nside = 2
+        while 4 * nside < lmax:
+            nside *= 2
+        _, a_lm = hp.synfast(
+            cl,
+            nside,
+            alm=True,
+            lmax=lmax,
+            fwhm=fwhm.to_value(u.radian),
+            verbose=False,
+        )
     return a_lm
 
 
@@ -396,8 +410,9 @@ def create_fake_beam_alm(
     pol=True,
     separate_IQU=False,
 ):
+    # pick an nside >= lmax to be sure that the a_lm will be fairly accurate
     nside = 2
-    while 2 * nside < lmax:
+    while nside < lmax:
         nside *= 2
     npix = 12 * nside ** 2
     pix = np.arange(npix)
diff --git a/src/toast/tests/ops_sim_tod_conviqt.py b/src/toast/tests/ops_sim_tod_conviqt.py
index f78ceff2a..4891990bd 100644
--- a/src/toast/tests/ops_sim_tod_conviqt.py
+++ b/src/toast/tests/ops_sim_tod_conviqt.py
@@ -42,53 +42,61 @@ def setUp(self):
         self.fname_sky = os.path.join(self.outdir, "sky_alm.fits")
         self.fname_beam = os.path.join(self.outdir, "beam_alm.fits")
 
-        # Synthetic sky and beam (a_lm expansions)
-        self.slm = create_fake_sky_alm(self.lmax, self.fwhm_sky)
-        self.slm[1:] = 0  # No polarization
-        hp.write_alm(self.fname_sky, self.slm, lmax=self.lmax, overwrite=True)
-
-        self.blm = create_fake_beam_alm(
-            self.lmax,
-            self.mmax,
-            fwhm_x=self.fwhm_beam,
-            fwhm_y=self.fwhm_beam,
-        )
-        hp.write_alm(
-            self.fname_beam,
-            self.blm,
-            lmax=self.lmax,
-            mmax_in=self.mmax,
-            overwrite=True,
-        )
+        self.rank = 0
+        if self.comm is not None:
+            self.rank = self.comm.rank
+
+        if self.rank == 0:
+            # Synthetic sky and beam (a_lm expansions)
+            self.slm = create_fake_sky_alm(self.lmax, self.fwhm_sky)
+            self.slm[1:] = 0  # No polarization
+            hp.write_alm(self.fname_sky, self.slm, lmax=self.lmax, overwrite=True)
+
+            self.blm = create_fake_beam_alm(
+                self.lmax,
+                self.mmax,
+                fwhm_x=self.fwhm_beam,
+                fwhm_y=self.fwhm_beam,
+            )
+            hp.write_alm(
+                self.fname_beam,
+                self.blm,
+                lmax=self.lmax,
+                mmax_in=self.mmax,
+                overwrite=True,
+            )
 
-        blm_I, blm_Q, blm_U = create_fake_beam_alm(
-            self.lmax,
-            self.mmax,
-            fwhm_x=self.fwhm_beam,
-            fwhm_y=self.fwhm_beam,
-            separate_IQU=True,
-        )
-        hp.write_alm(
-            self.fname_beam.replace(".fits", "_I000.fits"),
-            blm_I,
-            lmax=self.lmax,
-            mmax_in=self.mmax,
-            overwrite=True,
-        )
-        hp.write_alm(
-            self.fname_beam.replace(".fits", "_0I00.fits"),
-            blm_Q,
-            lmax=self.lmax,
-            mmax_in=self.mmax,
-            overwrite=True,
-        )
-        hp.write_alm(
-            self.fname_beam.replace(".fits", "_00I0.fits"),
-            blm_U,
-            lmax=self.lmax,
-            mmax_in=self.mmax,
-            overwrite=True,
-        )
+            blm_I, blm_Q, blm_U = create_fake_beam_alm(
+                self.lmax,
+                self.mmax,
+                fwhm_x=self.fwhm_beam,
+                fwhm_y=self.fwhm_beam,
+                separate_IQU=True,
+            )
+            hp.write_alm(
+                self.fname_beam.replace(".fits", "_I000.fits"),
+                blm_I,
+                lmax=self.lmax,
+                mmax_in=self.mmax,
+                overwrite=True,
+            )
+            hp.write_alm(
+                self.fname_beam.replace(".fits", "_0I00.fits"),
+                blm_Q,
+                lmax=self.lmax,
+                mmax_in=self.mmax,
+                overwrite=True,
+            )
+            hp.write_alm(
+                self.fname_beam.replace(".fits", "_00I0.fits"),
+                blm_U,
+                lmax=self.lmax,
+                mmax_in=self.mmax,
+                overwrite=True,
+            )
+
+        if self.comm is not None:
+            self.comm.Barrier()
 
         return
 
@@ -154,13 +162,9 @@ def test_sim(self):
         toast_bin_path = os.path.join(self.outdir, "toast_bin.fits")
         write_healpix_fits(data[binner.binned], toast_bin_path, nest=False)
 
-        rank = 0
-        if self.comm is not None:
-            rank = self.comm.rank
-
         fail = False
 
-        if rank == 0:
+        if self.rank == 0:
             import matplotlib.pyplot as plt
 
             mapfile = os.path.join(self.outdir, "toast_bin.fits")
@@ -194,38 +198,45 @@ def test_sim(self):
                 pol=False,
             )
 
-            sky = hp.alm2map(self.slm[0], self.nside, lmax=self.lmax, verbose=False)
-            beam = hp.alm2map(
-                self.blm[0], self.nside, lmax=self.lmax, mmax=self.mmax, verbose=False
-            )
-
-            fig = plt.figure(figsize=[12, 8])
-            nrow, ncol = 2, 3
-            hp.mollview(sky, title="input sky", sub=[nrow, ncol, 1])
-            hp.mollview(mdata, title="output sky", sub=[nrow, ncol, 2])
-            hp.mollview(smoothed, title="smoothed sky", sub=[nrow, ncol, 3])
-            hp.mollview(beam, title="beam", sub=[nrow, ncol, 4], rot=[0, 90])
-
-            ell = np.arange(self.lmax + 1)
-            ax = fig.add_subplot(nrow, ncol, 5)
-            ax.plot(ell[1:], cl_in[1:], label="input")
-            ax.plot(ell[1:], cl_smoothed[1:], label="smoothed")
-            ax.plot(ell[1:], blsq[1:], label="beam")
-            ax.plot(ell[1:], gauss_blsq[1:], label="gauss beam")
-            ax.plot(ell[1:], 1 / deconv[1:] ** 2, label="1 / deconv")
-            ax.plot(
-                ell[1:],
-                cl_in[1:] * blsq[1:] * deconv[1:] ** 2,
-                label="input x beam x deconv",
-            )
-            ax.plot(ell[1:], cl_out[1:], label="output")
-            ax.legend(loc="best")
-            ax.set_xscale("log")
-            ax.set_yscale("log")
-            ax.set_ylim([1e-20, 1e1])
-
-            outfile = os.path.join(self.outdir, "cl_comparison.png")
-            fig.savefig(outfile)
+            if self.comm is None or self.comm.size == 1:
+                sky = hp.alm2map(self.slm[0], self.nside, lmax=self.lmax, verbose=False)
+                beam = hp.alm2map(
+                    self.blm[0],
+                    self.nside,
+                    lmax=self.lmax,
+                    mmax=self.mmax,
+                    verbose=False,
+                )
+
+                fig = plt.figure(figsize=[12, 8])
+                nrow, ncol = 2, 3
+                hp.mollview(sky, title="input sky", sub=[nrow, ncol, 1])
+                hp.mollview(mdata, title="output sky", sub=[nrow, ncol, 2])
+                hp.mollview(smoothed, title="smoothed sky", sub=[nrow, ncol, 3])
+                hp.mollview(beam, title="beam", sub=[nrow, ncol, 4], rot=[0, 90])
+
+                ell = np.arange(self.lmax + 1)
+                ax = fig.add_subplot(nrow, ncol, 5)
+                ax.plot(ell[1:], cl_in[1:], label="input")
+                ax.plot(ell[1:], cl_smoothed[1:], label="smoothed")
+                ax.plot(ell[1:], blsq[1:], label="beam")
+                ax.plot(ell[1:], gauss_blsq[1:], label="gauss beam")
+                ax.plot(ell[1:], 1 / deconv[1:] ** 2, label="1 / deconv")
+                ax.plot(
+                    ell[1:],
+                    cl_in[1:] * blsq[1:] * deconv[1:] ** 2,
+                    label="input x beam x deconv",
+                )
+                ax.plot(ell[1:], cl_out[1:], label="output")
+                ax.legend(loc="best")
+                ax.set_xscale("log")
+                ax.set_yscale("log")
+                ax.set_ylim([1e-20, 1e1])
+
+                # For some reason, matplotlib hangs with multiple tasks,
+                # even if only one writes.
+                outfile = os.path.join(self.outdir, "cl_comparison.png")
+                fig.savefig(outfile)
 
             compare = blsq > 1e-5
             ref = cl_in[compare] * blsq[compare] * deconv[compare] ** 2
@@ -238,13 +249,15 @@ def test_sim(self):
             ):
                 fail = True
 
-        if data.comm.comm_world is not None:
-            fail = data.comm.comm_world.bcast(fail, root=0)
+        if self.comm is not None:
+            fail = self.comm.bcast(fail, root=0)
 
         self.assertFalse(fail)
 
         return
 
+    """
+
     def test_sim_hwp(self):
         if not ops.conviqt.available():
             print("libconviqt not available, skipping tests")
@@ -268,3 +281,5 @@ def test_sim_hwp(self):
         sim_conviqt.exec(data)
 
         return
+
+    """
diff --git a/src/toast/tests/ops_sim_tod_totalconvolve.py b/src/toast/tests/ops_sim_tod_totalconvolve.py
new file mode 100644
index 000000000..d3f7d4ee9
--- /dev/null
+++ b/src/toast/tests/ops_sim_tod_totalconvolve.py
@@ -0,0 +1,467 @@
+# Copyright (c) 2015-2020 by the parties listed in the AUTHORS file.
+# All rights reserved.  Use of this source code is governed by
+# a BSD-style license that can be found in the LICENSE file.
+
+import os
+
+import numpy as np
+
+from astropy import units as u
+
+import healpy as hp
+
+from .mpi import MPITestCase
+
+from ..vis import set_matplotlib_backend
+
+from .. import qarray as qa
+
+from .. import ops as ops
+
+from ..pixels_io import write_healpix_fits
+
+from ._helpers import (
+    create_outdir,
+    create_healpix_ring_satellite,
+    create_satellite_data,
+    create_fake_sky_alm,
+    create_fake_beam_alm,
+)
+
+
+class SimTotalconvolveTest(MPITestCase):
+    def setUp(self):
+        fixture_name = os.path.splitext(os.path.basename(__file__))[0]
+        self.outdir = create_outdir(self.comm, fixture_name)
+
+        self.nside = 128
+        self.lmax = 128
+        self.fwhm_sky = 10 * u.degree
+        self.fwhm_beam = 30 * u.degree
+        self.mmax = self.lmax
+        self.fname_sky = os.path.join(self.outdir, "sky_alm.fits")
+        self.fname_beam = os.path.join(self.outdir, "beam_alm.fits")
+        # Asymmetric beam for Conviqt comparison
+        self.fname_beam_asym = os.path.join(self.outdir, "beam_alm.asym.fits")
+        # Point source sky for Conviqt comparison
+        self.fname_sky_ps = os.path.join(self.outdir, "sky_alm.ps.fits")
+        myrank = 0 if self.comm is None else self.comm.rank
+
+        # Synthetic sky and beam (a_lm expansions)
+        self.slm = create_fake_sky_alm(self.lmax, self.fwhm_sky)
+        self.slm[1:] = 0  # No polarization
+        if myrank == 0:
+            hp.write_alm(self.fname_sky, self.slm, lmax=self.lmax, overwrite=True)
+
+        self.slm_ps = create_fake_sky_alm(self.lmax, self.fwhm_sky, pointsources=True)
+        self.slm_ps[1:] = 0  # No polarization
+        if myrank == 0:
+            hp.write_alm(self.fname_sky_ps, self.slm_ps, lmax=self.lmax, overwrite=True)
+
+        self.blm = create_fake_beam_alm(
+            self.lmax,
+            self.mmax,
+            fwhm_x=self.fwhm_beam,
+            fwhm_y=self.fwhm_beam,
+        )
+        if myrank == 0:
+            hp.write_alm(
+                self.fname_beam,
+                self.blm,
+                lmax=self.lmax,
+                mmax_in=self.mmax,
+                overwrite=True,
+            )
+
+        self.blm_asym = create_fake_beam_alm(
+            self.lmax,
+            self.mmax,
+            fwhm_x=self.fwhm_beam * 0.5,
+            fwhm_y=self.fwhm_beam,
+        )
+        if myrank == 0:
+            hp.write_alm(
+                self.fname_beam_asym,
+                self.blm_asym,
+                lmax=self.lmax,
+                mmax_in=self.mmax,
+                overwrite=True,
+            )
+
+        blm_I, blm_Q, blm_U = create_fake_beam_alm(
+            self.lmax,
+            self.mmax,
+            fwhm_x=self.fwhm_beam,
+            fwhm_y=self.fwhm_beam,
+            separate_IQU=True,
+        )
+        if myrank == 0:
+            hp.write_alm(
+                self.fname_beam.replace(".fits", "_I000.fits"),
+                blm_I,
+                lmax=self.lmax,
+                mmax_in=self.mmax,
+                overwrite=True,
+            )
+            hp.write_alm(
+                self.fname_beam.replace(".fits", "_0I00.fits"),
+                blm_Q,
+                lmax=self.lmax,
+                mmax_in=self.mmax,
+                overwrite=True,
+            )
+            hp.write_alm(
+                self.fname_beam.replace(".fits", "_00I0.fits"),
+                blm_U,
+                lmax=self.lmax,
+                mmax_in=self.mmax,
+                overwrite=True,
+            )
+
+        if self.comm is not None:
+            self.comm.Barrier()
+
+        return
+
+    def test_conviqt(self):
+        if not ops.totalconvolve.available():
+            print("ducc0.totalconvolve not available, skipping tests")
+            return
+
+        if not ops.conviqt.available():
+            print("libconviqt not available, skipping tests")
+            return
+
+        # Create a fake scan strategy that hits every pixel once.
+        data = create_healpix_ring_satellite(self.comm, nside=self.nside)
+
+        # Generate timestreams
+
+        detpointing = ops.PointingDetectorSimple()
+
+        totalconvolve_key = "totalconvolve_tod"
+        sim_totalconvolve = ops.SimTotalconvolve(
+            comm=self.comm,
+            detector_pointing=detpointing,
+            sky_file=self.fname_sky_ps,
+            beam_file=self.fname_beam_asym,
+            dxx=False,
+            det_data=totalconvolve_key,
+            normalize_beam=True,
+            fwhm=self.fwhm_sky,
+        )
+        sim_totalconvolve.exec(data)
+
+        conviqt_key = "conviqt_tod"
+        sim_conviqt = ops.SimConviqt(
+            comm=self.comm,
+            detector_pointing=detpointing,
+            sky_file=self.fname_sky_ps,
+            beam_file=self.fname_beam_asym,
+            dxx=False,
+            det_data=conviqt_key,
+            normalize_beam=True,
+            fwhm=self.fwhm_sky,
+        )
+        sim_conviqt.exec(data)
+
+        # Bin both signals into maps
+
+        pointing = ops.PointingHealpix(
+            nside=self.nside,
+            nest=False,
+            mode="I",
+            detector_pointing=detpointing,
+        )
+        pointing.apply(data)
+
+        default_model = ops.DefaultNoiseModel()
+        default_model.apply(data)
+
+        cov_and_hits = ops.CovarianceAndHits(
+            pixel_dist="pixel_dist",
+            pointing=pointing,
+            noise_model=default_model.noise_model,
+            rcond_threshold=1.0e-6,
+            sync_type="alltoallv",
+        )
+        cov_and_hits.apply(data)
+
+        binner = ops.BinMap(
+            pixel_dist="pixel_dist",
+            covariance=cov_and_hits.covariance,
+            det_data=totalconvolve_key,
+            pointing=pointing,
+            noise_model=default_model.noise_model,
+            sync_type="alltoallv",
+        )
+        binner.apply(data)
+        path_totalconvolve = os.path.join(self.outdir, "toast_bin.totalconvolve.fits")
+        write_healpix_fits(data[binner.binned], path_totalconvolve, nest=False)
+
+        binner = ops.BinMap(
+            pixel_dist="pixel_dist",
+            covariance=cov_and_hits.covariance,
+            det_data=conviqt_key,
+            pointing=pointing,
+            noise_model=default_model.noise_model,
+            sync_type="alltoallv",
+        )
+        binner.apply(data)
+        path_conviqt = os.path.join(self.outdir, "toast_bin.conviqt.fits")
+        write_healpix_fits(data[binner.binned], path_conviqt, nest=False)
+
+        rank = 0
+        if self.comm is not None:
+            rank = self.comm.rank
+
+        fail = False
+
+        if rank == 0:
+            import matplotlib.pyplot as plt
+
+            sky = hp.alm2map(self.slm_ps[0], self.nside, lmax=self.lmax, verbose=False)
+            beam = hp.alm2map(
+                self.blm_asym[0],
+                self.nside,
+                lmax=self.lmax,
+                mmax=self.mmax,
+                verbose=False,
+            )
+
+            map_totalconvolve = hp.read_map(path_totalconvolve)
+            map_conviqt = hp.read_map(path_conviqt)
+
+            # For some reason, matplotlib hangs with multiple tasks,
+            # even if only one writes.
+            if self.comm is None or self.comm.size == 1:
+                fig = plt.figure(figsize=[12, 8])
+                nrow, ncol = 2, 2
+                hp.mollview(sky, title="input sky", sub=[nrow, ncol, 1])
+                hp.mollview(beam, title="beam", sub=[nrow, ncol, 2], rot=[0, 90])
+                amp = np.amax(map_totalconvolve) / 4
+                hp.mollview(
+                    map_totalconvolve,
+                    min=-amp,
+                    max=amp,
+                    title="totalconvolve",
+                    sub=[nrow, ncol, 3],
+                )
+                hp.mollview(
+                    map_conviqt,
+                    min=-amp,
+                    max=amp,
+                    title="conviqt",
+                    sub=[nrow, ncol, 4],
+                )
+                outfile = os.path.join(self.outdir, "map_comparison.png")
+                fig.savefig(outfile)
+
+            for obs in data.obs:
+                for det in obs.local_detectors:
+                    tod_totalconvolve = obs.detdata[totalconvolve_key][det]
+                    tod_conviqt = obs.detdata[conviqt_key][det]
+                    if not np.allclose(
+                        tod_totalconvolve,
+                        tod_conviqt,
+                        rtol=1e-3,
+                        atol=1e-3,
+                    ):
+                        import matplotlib.pyplot as plt
+                        import pdb
+
+                        pdb.set_trace()
+                        fail = True
+                        break
+                if fail:
+                    break
+
+        if data.comm.comm_world is not None:
+            fail = data.comm.comm_world.bcast(fail, root=0)
+
+        self.assertFalse(fail)
+
+        return
+
+    def test_sim(self):
+        if not ops.totalconvolve.available():
+            print("ducc0.totalconvolve not available, skipping tests")
+            return
+
+        # Create a fake scan strategy that hits every pixel once.
+        data = create_healpix_ring_satellite(self.comm, nside=self.nside)
+
+        # Generate timestreams
+
+        detpointing = ops.PointingDetectorSimple()
+
+        key = "totalconvolve_tod"
+        sim_totalconvolve = ops.SimTotalconvolve(
+            comm=self.comm,
+            detector_pointing=detpointing,
+            sky_file=self.fname_sky,
+            beam_file=self.fname_beam,
+            dxx=False,
+            det_data=key,
+            normalize_beam=True,
+            fwhm=self.fwhm_sky,
+        )
+        sim_totalconvolve.exec(data)
+
+        # Bin a map to study
+
+        pointing = ops.PointingHealpix(
+            nside=self.nside,
+            nest=False,
+            mode="I",
+            detector_pointing=detpointing,
+        )
+        pointing.apply(data)
+
+        default_model = ops.DefaultNoiseModel()
+        default_model.apply(data)
+
+        cov_and_hits = ops.CovarianceAndHits(
+            pixel_dist="pixel_dist",
+            pointing=pointing,
+            noise_model=default_model.noise_model,
+            rcond_threshold=1.0e-6,
+            sync_type="alltoallv",
+        )
+        cov_and_hits.apply(data)
+
+        binner = ops.BinMap(
+            pixel_dist="pixel_dist",
+            covariance=cov_and_hits.covariance,
+            det_data=key,
+            pointing=pointing,
+            noise_model=default_model.noise_model,
+            sync_type="alltoallv",
+        )
+        binner.apply(data)
+
+        # Study the map on the root process
+
+        toast_bin_path = os.path.join(self.outdir, "toast_bin.fits")
+        write_healpix_fits(data[binner.binned], toast_bin_path, nest=False)
+
+        rank = 0
+        if self.comm is not None:
+            rank = self.comm.rank
+
+        fail = False
+
+        if rank == 0:
+            import matplotlib.pyplot as plt
+
+            mapfile = os.path.join(self.outdir, "toast_bin.fits")
+            mdata = hp.read_map(mapfile, nest=False)
+
+            deconv = 1 / hp.gauss_beam(
+                self.fwhm_sky.to_value(u.radian),
+                lmax=self.lmax,
+                pol=False,
+            )
+
+            smoothed = hp.alm2map(
+                hp.almxfl(self.slm[0], deconv),
+                self.nside,
+                lmax=self.lmax,
+                fwhm=self.fwhm_beam.to_value(u.radian),
+                verbose=False,
+                pixwin=False,
+            )
+
+            cl_out = hp.anafast(mdata, lmax=self.lmax)
+            cl_smoothed = hp.anafast(smoothed, lmax=self.lmax)
+
+            cl_in = hp.alm2cl(self.slm[0], lmax=self.lmax)
+            blsq = hp.alm2cl(self.blm[0], lmax=self.lmax, mmax=self.mmax)
+            blsq /= blsq[1]
+
+            gauss_blsq = hp.gauss_beam(
+                self.fwhm_beam.to_value(u.radian),
+                lmax=self.lmax,
+                pol=False,
+            )
+
+            sky = hp.alm2map(self.slm[0], self.nside, lmax=self.lmax, verbose=False)
+            beam = hp.alm2map(
+                self.blm[0], self.nside, lmax=self.lmax, mmax=self.mmax, verbose=False
+            )
+
+            # For some reason, matplotlib hangs with multiple tasks,
+            # even if only one writes.
+            if self.comm is None or self.comm.size == 1:
+                fig = plt.figure(figsize=[12, 8])
+                nrow, ncol = 2, 3
+                hp.mollview(sky, title="input sky", sub=[nrow, ncol, 1])
+                hp.mollview(mdata, title="output sky", sub=[nrow, ncol, 2])
+                hp.mollview(smoothed, title="smoothed sky", sub=[nrow, ncol, 3])
+                hp.mollview(beam, title="beam", sub=[nrow, ncol, 4], rot=[0, 90])
+
+                ell = np.arange(self.lmax + 1)
+                ax = fig.add_subplot(nrow, ncol, 5)
+                ax.plot(ell[1:], cl_in[1:], label="input")
+                ax.plot(ell[1:], cl_smoothed[1:], label="smoothed")
+                ax.plot(ell[1:], blsq[1:], label="beam")
+                ax.plot(ell[1:], gauss_blsq[1:], label="gauss beam")
+                ax.plot(ell[1:], 1 / deconv[1:] ** 2, label="1 / deconv")
+                ax.plot(
+                    ell[1:],
+                    cl_in[1:] * blsq[1:] * deconv[1:] ** 2,
+                    label="input x beam x deconv",
+                )
+                ax.plot(ell[1:], cl_out[1:], label="output")
+                ax.legend(loc="best")
+                ax.set_xscale("log")
+                ax.set_yscale("log")
+                ax.set_ylim([1e-20, 1e1])
+
+                if self.comm is None or self.comm.size == 1:
+                    # For some reason, matplotlib hangs with multiple tasks,
+                    # even if only one writes.
+                    outfile = os.path.join(self.outdir, "cl_comparison.png")
+                    fig.savefig(outfile)
+
+            compare = blsq > 1e-5
+            ref = cl_in[compare] * blsq[compare] * deconv[compare] ** 2
+            norm = np.mean(cl_out[compare] / ref)
+            if not np.allclose(
+                norm * ref,
+                cl_out[compare],
+                rtol=1e-5,
+                atol=1e-5,
+            ):
+                fail = True
+
+        if data.comm.comm_world is not None:
+            fail = data.comm.comm_world.bcast(fail, root=0)
+
+        self.assertFalse(fail)
+
+        return
+
+    def test_sim_hwp(self):
+        if not ops.totalconvolve.available():
+            print("ducc0.totalconvolve not available, skipping tests")
+            return
+
+        # Create a fake scan strategy that hits every pixel once.
+        data = create_satellite_data(self.comm)
+
+        # make a simple pointing matrix
+        detpointing = ops.PointingDetectorSimple()
+
+        # Generate timestreams
+        sim_totalconvolve = ops.SimTotalconvolve(
+            comm=self.comm,
+            detector_pointing=detpointing,
+            sky_file=self.fname_sky,
+            beam_file=self.fname_beam,
+            dxx=False,
+            hwp_angle="hwp_angle",
+        )
+        sim_totalconvolve.exec(data)
+
+        return
diff --git a/src/toast/tests/runner.py b/src/toast/tests/runner.py
index 30df86f0b..9d90c2262 100644
--- a/src/toast/tests/runner.py
+++ b/src/toast/tests/runner.py
@@ -39,6 +39,7 @@
 from . import ops_sim_tod_noise as test_ops_sim_tod_noise
 from . import ops_sim_tod_dipole as test_ops_sim_tod_dipole
 from . import ops_sim_tod_conviqt as test_ops_sim_tod_conviqt
+from . import ops_sim_tod_totalconvolve as test_ops_sim_tod_totalconvolve
 from . import ops_mapmaker_utils as test_ops_mapmaker_utils
 from . import ops_mapmaker_binning as test_ops_mapmaker_binning
 from . import ops_mapmaker_solve as test_ops_mapmaker_solve
@@ -154,6 +155,7 @@ def test(name=None, verbosity=2):
         suite.addTest(loader.loadTestsFromModule(test_ops_sim_tod_noise))
         suite.addTest(loader.loadTestsFromModule(test_ops_sim_tod_dipole))
         suite.addTest(loader.loadTestsFromModule(test_ops_sim_tod_conviqt))
+        suite.addTest(loader.loadTestsFromModule(test_ops_sim_tod_totalconvolve))
         suite.addTest(loader.loadTestsFromModule(test_ops_mapmaker_utils))
         suite.addTest(loader.loadTestsFromModule(test_ops_mapmaker_binning))
         suite.addTest(loader.loadTestsFromModule(test_ops_mapmaker_solve))