laplace_operator.py 4.98 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
# 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/>.
#
Martin Reinecke's avatar
Martin Reinecke committed
14
# Copyright(C) 2013-2018 Max-Planck-Society
15
16
17
#
# 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
from ..field import Field
Martin Reinecke's avatar
Martin Reinecke committed
21
from ..domains.power_space import PowerSpace
Martin Reinecke's avatar
Martin Reinecke committed
22
from .endomorphic_operator import EndomorphicOperator
23
from ..utilities import infer_space
Martin Reinecke's avatar
Martin Reinecke committed
24
from ..domain_tuple 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

    Parameters
    ----------
Martin Reinecke's avatar
Martin Reinecke committed
38
    logarithmic : bool, optional
Jakob Knollmueller's avatar
Jakob Knollmueller committed
39
40
        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
        pos = self.domain[self._space].k_lengths.copy()
Martin Reinecke's avatar
Martin Reinecke committed
54
        if logarithmic:
Martin Reinecke's avatar
Martin Reinecke committed
55
56
            pos[1:] = np.log(pos[1:])
            pos[0] = pos[1]-1.
57

58
59
60
61
62
        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
63
64
65
        self._dposc[-1] = 0.
        self._dposc[1:] += self._dpos
        self._dposc *= 0.5
Jakob Knollmueller's avatar
Jakob Knollmueller committed
66
67
68
69
70
71

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

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

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

    def _adjoint_times(self, x):
        axes = x.domain.axes[self._space]
103
        axis = axes[0]
104
        nval = len(self._dposc)
105
        prefix = (slice(None),) * axis
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
106
107
        sl_l = prefix + (slice(None, -1),)  # "left" slice
        sl_r = prefix + (slice(1, None),)  # "right" slice
108
109
        dpos = self._dpos.reshape((1,)*axis + (nval-1,))
        dposc = self._dposc.reshape((1,)*axis + (nval,))
110
        yf = x.val
Martin Reinecke's avatar
Martin Reinecke committed
111
112
        if axis == dobj.distaxis(yf):
            yf = dobj.redistribute(yf, nodist=(axis,))
Martin Reinecke's avatar
Martin Reinecke committed
113
        y = dobj.local_data(yf)
Martin Reinecke's avatar
Martin Reinecke committed
114
        y = y/dposc
115
116
117
        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
118
        ret = np.empty_like(y)
119
120
121
        ret[sl_l] = deriv
        ret[prefix + (-1,)] = 0.
        ret[sl_r] -= deriv
Martin Reinecke's avatar
Martin Reinecke committed
122
        ret = dobj.from_local_data(x.shape, ret, dobj.distaxis(yf))
Martin Reinecke's avatar
Martin Reinecke committed
123
        if dobj.distaxis(yf) != dobj.distaxis(x.val):
Martin Reinecke's avatar
Martin Reinecke committed
124
            ret = dobj.redistribute(ret, dist=dobj.distaxis(x.val))
125
        return Field(self.domain, val=ret)
Martin Reinecke's avatar
Martin Reinecke committed
126
127
128
129
130
131

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