laplace_operator.py 5.19 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18 19

import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
20 21 22
from ..field import Field
from ..spaces.power_space import PowerSpace
from .endomorphic_operator import EndomorphicOperator
23
from ..utilities import infer_space
Martin Reinecke's avatar
Martin Reinecke committed
24
from .. import DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
25
from .. import dobj
Jakob Knollmueller's avatar
Jakob Knollmueller committed
26 27


28
class LaplaceOperator(EndomorphicOperator):
29
    """An irregular LaplaceOperator with free boundary and excluding monopole.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
30

Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
31
    This LaplaceOperator implements the second derivative of a Field in
Martin Reinecke's avatar
Martin Reinecke committed
32
    PowerSpace on logarithmic or linear scale with vanishing curvature at the
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
33 34
    boundary, starting at the second entry of the Field. The second derivative
    of the Field on the irregular grid is calculated using finite differences.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
35 36 37 38 39 40

    Parameters
    ----------
    logarithmic : boolean,
        Whether smoothness is calculated on a logarithmic scale or linear scale
        default : True
Martin Reinecke's avatar
Martin Reinecke committed
41 42
    space : int
        The index of the domain on which the operator acts
Jakob Knollmueller's avatar
Jakob Knollmueller committed
43 44
    """

45 46
    def __init__(self, domain, space=None, logarithmic=True):
        super(LaplaceOperator, self).__init__()
Martin Reinecke's avatar
Martin Reinecke committed
47
        self._domain = DomainTuple.make(domain)
48
        self._space = infer_space(self._domain, space)
49 50 51

        if not isinstance(self._domain[self._space], PowerSpace):
            raise ValueError("Operator must act on a PowerSpace.")
52

53 54
        self._logarithmic = bool(logarithmic)

55
        pos = self.domain[self._space].k_lengths.copy()
56
        if self.logarithmic:
Martin Reinecke's avatar
Martin Reinecke committed
57 58
            pos[1:] = np.log(pos[1:])
            pos[0] = pos[1]-1.
59

60 61 62 63 64
        self._dpos = pos[1:]-pos[:-1]  # defined between points
        # centered distances (also has entries for the first and last point
        # for convenience, but they will never affect the result)
        self._dposc = np.empty_like(pos)
        self._dposc[:-1] = self._dpos
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
65 66 67
        self._dposc[-1] = 0.
        self._dposc[1:] += self._dpos
        self._dposc *= 0.5
Jakob Knollmueller's avatar
Jakob Knollmueller committed
68 69 70 71 72 73

    @property
    def domain(self):
        return self._domain

    @property
Martin Reinecke's avatar
Martin Reinecke committed
74 75
    def capability(self):
        return self.TIMES | self.ADJOINT_TIMES
Jakob Knollmueller's avatar
Jakob Knollmueller committed
76

77 78 79 80
    @property
    def logarithmic(self):
        return self._logarithmic

81 82
    def _times(self, x):
        axes = x.domain.axes[self._space]
83
        axis = axes[0]
Martin Reinecke's avatar
Martin Reinecke committed
84 85 86 87
        locval = x.val
        if axis == dobj.distaxis(locval):
            locval = dobj.redistribute(locval, nodist=(axis,))
        val = dobj.local_data(locval)
88
        nval = len(self._dposc)
89
        prefix = (slice(None),) * axis
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
90 91
        sl_l = prefix + (slice(None, -1),)  # "left" slice
        sl_r = prefix + (slice(1, None),)  # "right" slice
92 93
        dpos = self._dpos.reshape((1,)*axis + (nval-1,))
        dposc = self._dposc.reshape((1,)*axis + (nval,))
Martin Reinecke's avatar
Martin Reinecke committed
94 95
        deriv = (val[sl_r]-val[sl_l])/dpos  # defined between points
        ret = np.empty_like(val)
96 97 98
        ret[sl_l] = deriv
        ret[prefix + (-1,)] = 0.
        ret[sl_r] -= deriv
Martin Reinecke's avatar
Martin Reinecke committed
99
        ret /= np.sqrt(dposc)
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
100
        ret[prefix + (slice(None, 2),)] = 0.
101
        ret[prefix + (-1,)] = 0.
Martin Reinecke's avatar
Martin Reinecke committed
102
        ret = dobj.from_local_data(locval.shape, ret, dobj.distaxis(locval))
Martin Reinecke's avatar
Martin Reinecke committed
103
        if dobj.distaxis(locval) != dobj.distaxis(x.val):
Martin Reinecke's avatar
Martin Reinecke committed
104 105
            ret = dobj.redistribute(ret, dist=dobj.distaxis(x.val))
        return Field(self.domain, val=ret).weight(-0.5, spaces=self._space)
106 107 108

    def _adjoint_times(self, x):
        axes = x.domain.axes[self._space]
109
        axis = axes[0]
110
        nval = len(self._dposc)
111
        prefix = (slice(None),) * axis
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
112 113
        sl_l = prefix + (slice(None, -1),)  # "left" slice
        sl_r = prefix + (slice(1, None),)  # "right" slice
114 115
        dpos = self._dpos.reshape((1,)*axis + (nval-1,))
        dposc = self._dposc.reshape((1,)*axis + (nval,))
Martin Reinecke's avatar
Martin Reinecke committed
116 117 118
        yf = x.weight(0.5, spaces=self._space).val
        if axis == dobj.distaxis(yf):
            yf = dobj.redistribute(yf, nodist=(axis,))
Martin Reinecke's avatar
Martin Reinecke committed
119
        y = dobj.local_data(yf)
Martin Reinecke's avatar
Martin Reinecke committed
120
        y /= np.sqrt(dposc)
121 122 123
        y[prefix + (slice(None, 2),)] = 0.
        y[prefix + (-1,)] = 0.
        deriv = (y[sl_r]-y[sl_l])/dpos  # defined between points
Martin Reinecke's avatar
Martin Reinecke committed
124
        ret = np.empty_like(y)
125 126 127
        ret[sl_l] = deriv
        ret[prefix + (-1,)] = 0.
        ret[sl_r] -= deriv
Martin Reinecke's avatar
Martin Reinecke committed
128
        ret = dobj.from_local_data(x.shape, ret, dobj.distaxis(yf))
Martin Reinecke's avatar
Martin Reinecke committed
129
        if dobj.distaxis(yf) != dobj.distaxis(x.val):
Martin Reinecke's avatar
Martin Reinecke committed
130 131
            ret = dobj.redistribute(ret, dist=dobj.distaxis(x.val))
        return Field(self.domain, val=ret).weight(-1, spaces=self._space)
Martin Reinecke's avatar
Martin Reinecke committed
132 133 134 135 136 137

    def apply(self, x, mode):
        self._check_input(x, mode)
        if mode == self.TIMES:
            return self._times(x)
        return self._adjoint_times(x)