Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

# -*- coding: utf-8 -*- 

# This file is part of the pyMOR project (http://www.pymor.org). 

# Copyright Holders: Rene Milk, Stephan Rave, Felix Schindler 

# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) 

 

from __future__ import absolute_import, division, print_function 

 

from itertools import chain 

 

from pymor.algorithms.timestepping import TimeStepperInterface 

from pymor.discretizations.interfaces import DiscretizationInterface 

from pymor.la import induced_norm, VectorArrayInterface 

from pymor.operators import OperatorInterface 

from pymor.operators.constructions import VectorOperator 

from pymor.parameters import Parameter 

from pymor.tools import method_arguments, FrozenDict 

 

 

class DiscretizationBase(DiscretizationInterface): 

    '''Base class for |Discretizations| providing some common functionality.''' 

 

    def __init__(self, operators, functionals, vector_operators, products=None, estimator=None, visualizer=None, 

                 cache_region='disk', name=None): 

        self.operators = FrozenDict(operators) 

        self.functionals = FrozenDict(functionals) 

        self.vector_operators = FrozenDict(vector_operators) 

        self.linear = all(op is None or op.linear for op in chain(operators.itervalues(), functionals.itervalues())) 

        self.products = products 

        self.estimator = estimator 

        self.visualizer = visualizer 

        self.cache_region = cache_region 

        self.name = name 

 

        if products: 

            for k, v in products.iteritems(): 

                setattr(self, '{}_product'.format(k), v) 

                setattr(self, '{}_norm'.format(k), induced_norm(v)) 

        if estimator is not None: 

            self.estimate = self.__estimate 

        if visualizer is not None: 

            self.visualize = self.__visualize 

 

    def __visualize(self, U, *args, **kwargs): 

        self.visualizer.visualize(U, self, *args, **kwargs) 

 

    def __estimate(self, U, mu=None): 

        return self.estimator.estimate(U, mu=mu, discretization=self) 

 

 

class StationaryDiscretization(DiscretizationBase): 

    '''Generic class for discretizations of stationary problems. 

 

    This class describes discrete problems given by the equation:: 

 

        L(u(μ), μ) = F(μ) 

 

    with a linear functional F and a (possibly) non-linear operator L. 

 

    Parameters 

    ---------- 

    operator 

        The |Operator| L. 

    rhs 

        The |Functional| F. 

    products 

        A dict of inner product |Operators| defined on the discrete space the 

        problem is posed on. For each product a corresponding norm 

        is added as a method of the discretization. 

    parameter_space 

        The |ParameterSpace| for which the discrete problem is posed. 

    estimator 

        An error estimator for the problem. This can be any object with 

        an `estimate(U, mu, discretization)` method. If `estimator` is 

        not `None` an `estimate(U, mu)` method is added to the 

        discretization. 

    visualizer 

        A visualizer for the problem. This can be any object with 

        a `visualize(U, discretization, ...)` method. If `visualizer` 

        is not `None` a `visualize(U, *args, **kwargs)` method is added 

        to the discretization, which forwards its arguments to the 

        visualizer's `visualize` method. 

    cache_region 

        `None` or name of the cache region to use. See 

        :mod:`pymor.core.cache`. 

    name 

        Name of the discretization. 

 

    Attributes 

    ---------- 

    operator 

        The |Operator| L. Synonymous for `operators['operator']`. 

    rhs 

        The |Functional| F. Synonymous for `functionals['rhs']`. 

    ''' 

 

    sid_ignore = ('visualizer', 'cache_region', 'name') 

 

    def __init__(self, operator, rhs, products=None, parameter_space=None, estimator=None, visualizer=None, 

                 cache_region='disk', name=None): 

        assert isinstance(operator, OperatorInterface) 

        assert isinstance(rhs, OperatorInterface) and rhs.linear 

        assert operator.dim_source == operator.dim_range == rhs.dim_source 

        assert rhs.dim_range == 1 

 

        operators = {'operator': operator} 

        functionals = {'rhs': rhs} 

        super(StationaryDiscretization, self).__init__(operators=operators, functionals=functionals, 

                                                       vector_operators={}, products=products, estimator=estimator, 

                                                       visualizer=visualizer, cache_region=cache_region, name=name) 

        self.dim_solution = operator.dim_source 

        self.type_solution = operator.type_source 

        self.operator = operator 

        self.rhs = rhs 

        self.operators = operators 

        self.build_parameter_type(inherits=(operator, rhs)) 

        self.parameter_space = parameter_space 

 

    with_arguments = frozenset(method_arguments(__init__)).union({'operators', 'functionals', 'vector_operators'}) 

 

    def with_(self, **kwargs): 

        assert set(kwargs.keys()) <= self.with_arguments 

        assert 'operators' not in kwargs or 'operator' not in kwargs 

        assert 'operators' not in kwargs or 'rhs' not in kwargs 

        assert 'operators' not in kwargs or kwargs['operators'].keys() == ['operator'] 

        assert 'functionals' not in kwargs or kwargs['functionals'].keys() == ['rhs'] 

        assert 'vector_operators' not in kwargs or not kwargs['vector_operators'].keys() 

 

        if 'operators' in kwargs: 

            kwargs.update(kwargs.pop('operators')) 

        if 'functionals' in kwargs: 

            kwargs.update(kwargs.pop('functionals')) 

        if 'vector_operators' in kwargs: 

            del kwargs['vector_operators'] 

 

        return self._with_via_init(kwargs) 

 

    def _solve(self, mu=None): 

        mu = self.parse_parameter(mu) 

 

        # explicitly checking if logging is disabled saves the str(mu) call 

        if not self.logging_disabled: 

            self.logger.info('Solving {} for {} ...'.format(self.name, mu)) 

 

        return self.operator.apply_inverse(self.rhs.as_vector(mu), mu=mu) 

 

 

class InstationaryDiscretization(DiscretizationBase): 

    '''Generic class for discretizations of stationary problems. 

 

    This class describes instationary problems given by the equations:: 

 

        M * ∂_t u(t, μ) + L(u(μ), t, μ) = F(t, μ) 

                                u(0, μ) = u_0(μ) 

 

    for t in [0,T], where L is a (possibly) non-linear time-dependent 

    |Operator|, F is a time-dependent linear |Functional|, and u_0 the 

    initial data. The mass operator is assumed to be linear, 

    time-independent and |Parameter|-independent. 

 

    Parameters 

    ---------- 

    T 

        The end-time T. 

    initial_data 

        The initial data u_0. Either a |VectorArray| of length 1 or 

        (for the |Parameter|-dependent case) a vector-like |Operator| 

        (i.e. a linear |Operator| with `dim_source == 1`). 

    operator 

        The |Operator| L. 

    rhs 

        The |Functional| F. 

    mass 

        The mass |Operator| `M`. If `None` the identity is assumed. 

    time_stepper 

        T time-stepper to be used by :meth:`solve`. Has to satisfy 

        the :class:`~pymor.algorithms.timestepping.TimeStepperInterface`. 

    num_values 

        The number of returned vectors of the solution trajectory. If `None`, each 

        intermediate vector that is calculated is returned. 

    products 

        A dict of product |Operators| defined on the discrete space the 

        problem is posed on. For each product a corresponding norm 

        is added as a method of the discretization. 

    parameter_space 

        The |ParameterSpace| for which the discrete problem is posed. 

    estimator 

        An error estimator for the problem. This can be any object with 

        an `estimate(U, mu, discretization)` method. If `estimator` is 

        not `None` an `estimate(U, mu)` method is added to the 

        discretization. 

    visualizer 

        A visualizer for the problem. This can be any object with 

        a `visualize(U, discretization, ...)` method. If `visualizer` 

        is not `None` a `visualize(U, *args, **kwargs)` method is added 

        to the discretization, which forwards its arguments to the 

        visualizer's `visualize` method. 

    cache_region 

        `None` or name of the cache region to use. See 

        :mod:`pymor.core.cache`. 

    name 

        Name of the discretization. 

 

    Attributes 

    ---------- 

    T 

        The end-time. 

    initial_data 

        The intial data u_0 given by a vector-like |Operator|. Synonymous 

        for `vector_operators['initial_data']`. 

    operator 

        The |Operator| L. Synonymous for `operators['operator']`. 

    rhs 

        The |Functional| F. Synonymous for `functionals['rhs']`. 

    mass 

        The mass operator M. Synonymous for `operators['mass']`. 

    time_stepper 

        The provided time-stepper. 

    ''' 

 

    sid_ignore = ('visualizer', 'cache_region', 'name') 

 

    def __init__(self, T, initial_data, operator, rhs=None, mass=None, time_stepper=None, num_values=None, 

                 products=None, parameter_space=None, estimator=None, visualizer=None, cache_region='disk', 

                 name=None): 

        assert isinstance(initial_data, (VectorArrayInterface, OperatorInterface)) 

        assert not isinstance(initial_data, OperatorInterface) or initial_data.dim_source == 1 

        assert isinstance(operator, OperatorInterface) 

        assert rhs is None or isinstance(rhs, OperatorInterface) and rhs.linear 

        assert mass is None or isinstance(mass, OperatorInterface) and mass.linear 

        if isinstance(initial_data, VectorArrayInterface): 

            initial_data = VectorOperator(initial_data, name='initial_data') 

        assert isinstance(time_stepper, TimeStepperInterface) 

        assert operator.dim_source == operator.dim_range == initial_data.dim_range 

        assert rhs is None or rhs.dim_source == operator.dim_source and rhs.dim_range == 1 

        assert mass is None or mass.dim_source == mass.dim_range == operator.dim_source 

 

        operators = {'operator': operator, 'mass': mass} 

        functionals = {'rhs': rhs} 

        vector_operators = {'initial_data': initial_data} 

        super(InstationaryDiscretization, self).__init__(operators=operators, functionals=functionals, 

                                                         vector_operators=vector_operators, 

                                                         products=products, estimator=estimator, 

                                                         visualizer=visualizer, cache_region=cache_region, name=name) 

        self.T = T 

        self.dim_solution = operator.dim_source 

        self.type_solution = operator.type_source 

        self.initial_data = initial_data 

        self.operator = operator 

        self.rhs = rhs 

        self.mass = mass 

        self.time_stepper = time_stepper 

        self.num_values = num_values 

        self.build_parameter_type(inherits=(initial_data, operator, rhs, mass), provides={'_t': 0}) 

        self.parameter_space = parameter_space 

 

        if hasattr(time_stepper, 'nt'): 

            self.with_arguments = self.with_arguments.union({'time_stepper_nt'}) 

 

    with_arguments = frozenset(method_arguments(__init__)).union({'operators', 'functionals', 'vector_operators'}) 

 

    def with_(self, **kwargs): 

        assert set(kwargs.keys()) <= self.with_arguments 

        assert 'operators' not in kwargs or kwargs['operators'].viewkeys() <= {'operator', 'mass'} 

        assert 'functionals' not in kwargs or kwargs['functionals'].viewkeys() <= {'rhs'} 

        assert 'vector_operators' not in kwargs or kwargs['vector_operators'].viewkeys() <= {'initial_data'} 

        assert 'operators' not in kwargs or not set(kwargs['operators']).intersection(kwargs.viewkeys()) 

        assert 'functionals' not in kwargs or not set(kwargs['functionals']).intersection(kwargs.viewkeys()) 

        assert 'vector_operators' not in kwargs or not set(kwargs['vector_operators']).intersection(kwargs.viewkeys()) 

        assert 'time_stepper_nt' not in kwargs or 'time_stepper' not in kwargs 

        if 'operators' in kwargs: 

            kwargs.update(kwargs.pop('operators')) 

        if 'functionals' in kwargs: 

            kwargs.update(kwargs.pop('functionals')) 

        if 'vector_operators' in kwargs: 

            kwargs.update(kwargs.pop('vector_operators')) 

        if 'time_stepper_nt' in kwargs: 

            kwargs['time_stepper'] = self.time_stepper.with_(nt=kwargs.pop('time_stepper_nt')) 

 

        return self._with_via_init(kwargs) 

 

    def _solve(self, mu=None): 

        mu = self.parse_parameter(mu).copy() if mu is not None else Parameter({}) 

 

        # explicitly checking if logging is disabled saves the expensive str(mu) call 

        if not self.logging_disabled: 

            self.logger.info('Solving {} for {} ...'.format(self.name, mu)) 

 

        mu['_t'] = 0 

        U0 = self.initial_data.as_vector(mu) 

        return self.time_stepper.solve(operator=self.operator, rhs=self.rhs, initial_data=U0, mass=self.mass, 

                                       initial_time=0, end_time=self.T, mu=mu, num_values=self.num_values)