laplace_operator.py 3.76 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/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
17

18
import numpy as np
19
20
21

from .. import dobj
from ..domain_tuple import DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
22
from ..domains.power_space import PowerSpace
23
from ..field import Field
24
from ..utilities import infer_space
25
from .endomorphic_operator import EndomorphicOperator
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
    def __init__(self, domain, space=None, logarithmic=True):
Martin Reinecke's avatar
Martin Reinecke committed
46
        self._domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
47
        self._capability = self.TIMES | self.ADJOINT_TIMES
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

Martin Reinecke's avatar
Martin Reinecke committed
67
68
    def apply(self, x, mode):
        self._check_input(x, mode)
69
        axes = x.domain.axes[self._space]
70
        axis = axes[0]
71
        nval = len(self._dposc)
72
        prefix = (slice(None),) * axis
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
73
74
        sl_l = prefix + (slice(None, -1),)  # "left" slice
        sl_r = prefix + (slice(1, None),)  # "right" slice
75
76
        dpos = self._dpos.reshape((1,)*axis + (nval-1,))
        dposc = self._dposc.reshape((1,)*axis + (nval,))
Martin Reinecke's avatar
Martin Reinecke committed
77
        v, val = dobj.ensure_not_distributed(x.val, (axis,))
Martin Reinecke's avatar
Martin Reinecke committed
78
        ret = np.empty_like(val)
Martin Reinecke's avatar
Martin Reinecke committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
        if mode == self.TIMES:
            deriv = (val[sl_r]-val[sl_l])/dpos  # defined between points
            ret[sl_l] = deriv
            ret[prefix + (-1,)] = 0.
            ret[sl_r] -= deriv
            ret /= dposc
            ret[prefix + (slice(None, 2),)] = 0.
            ret[prefix + (-1,)] = 0.
        else:
            val = val/dposc
            val[prefix + (slice(None, 2),)] = 0.
            val[prefix + (-1,)] = 0.
            deriv = (val[sl_r]-val[sl_l])/dpos  # defined between points
            ret[sl_l] = deriv
            ret[prefix + (-1,)] = 0.
            ret[sl_r] -= deriv
Martin Reinecke's avatar
Martin Reinecke committed
95
96
        ret = dobj.from_local_data(x.shape, ret, dobj.distaxis(v))
        return Field(self.domain, dobj.ensure_default_distributed(ret))