laplace_operator.py 5.21 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 ..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
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)