Coverage for /usr/local/lib/python3.11/dist-packages/grond/problems/dynamic_rupture/problem.py: 30%
351 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-26 16:25 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-26 16:25 +0000
1import numpy as num
2import logging
4from pyrocko import gf, util
5from pyrocko.guts import String, Float, Dict, Int, Bool, StringChoice
6from pyrocko.guts_array import Array
7from pyrocko.gf.seismosizer import map_anchor
9from grond.meta import expand_template, Parameter, has_get_plot_classes, \
10 GrondError
12from ..base import Problem, ProblemConfig
13from .. import CMTProblem
15guts_prefix = 'grond'
17logger = logging.getLogger('grond.problems.dynamic_rupture.problem')
18km = 1e3
19d2r = num.pi / 180.
20as_km = dict(scale_factor=km, scale_unit='km')
21as_gpa = dict(scale_factor=1e9, scale_unit='GPa')
24class DynamicRuptureProblemConfig(ProblemConfig):
26 ranges = Dict.T(String.T(), gf.Range.T())
28 anchor = StringChoice.T(
29 choices=tuple(map_anchor.keys()),
30 default='top')
32 decimation_factor = Int.T(
33 default=1,
34 help='Decimation factor of the discretized sub-faults')
35 distance_min = Float.T(default=0.)
36 adaptive_resolution = StringChoice.T(
37 choices=('off', 'linear', 'uniqueness'),
38 default='off')
39 adaptive_start = Int.T(
40 default=0)
42 pure_shear = Bool.T(
43 default=True)
44 smooth_rupture = Bool.T(
45 default=False)
47 patch_mask = Array.T(
48 optional=True,
49 shape=(None,),
50 dtype=bool,
51 serialize_as='list')
53 point_source_target_balancing = Bool.T(
54 default=False,
55 help='If ``True``, target balancing (if used) is performed on a '
56 'moment tensor point source at the events location. It increases '
57 'the speed, but might lead to not fully optimized target weights.'
58 )
60 def get_problem(self, event, target_groups, targets):
61 if self.decimation_factor != 1:
62 logger.warn(
63 'Decimation factor for dynamic rupture source set to %i.'
64 ' Results may be inaccurate.' % self.decimation_factor)
66 base_source = gf.PseudoDynamicRupture.from_pyrocko_event(
67 event,
68 anchor=self.anchor,
69 nx=int(self.ranges['nx'].start),
70 ny=int(self.ranges['ny'].start),
71 decimation_factor=self.decimation_factor,
72 pure_shear=self.pure_shear,
73 smooth_rupture=self.smooth_rupture,
74 patch_mask=self.patch_mask,
75 stf=None)
77 subs = dict(
78 event_name=event.name,
79 event_time=util.time_to_str(event.time))
81 cmt_problem = None
82 if self.point_source_target_balancing:
83 base_source_cmt = gf.MTSource.from_pyrocko_event(event)
85 stf = gf.HalfSinusoidSTF()
86 stf.duration = event.duration or 0.0
88 base_source_cmt.stf = stf
90 ranges = dict(
91 time=self.ranges['time'],
92 north_shift=self.ranges['north_shift'],
93 east_shift=self.ranges['east_shift'],
94 depth=self.ranges['depth'],
95 magnitude=gf.Range(
96 start=event.magnitude - 1.,
97 stop=event.magnitude + 1.),
98 duration=gf.Range(start=0., stop=stf.duration * 2.),
99 rmnn=gf.Range(start=-1.41421, stop=1.41421),
100 rmee=gf.Range(start=-1.41421, stop=1.41421),
101 rmdd=gf.Range(start=-1.41421, stop=1.41421),
102 rmne=gf.Range(start=-1., stop=1.),
103 rmnd=gf.Range(start=-1., stop=1.),
104 rmed=gf.Range(start=-1., stop=1.))
106 cmt_problem = CMTProblem(
107 name=expand_template(self.name_template, subs),
108 base_source=base_source_cmt,
109 distance_min=self.distance_min,
110 target_groups=target_groups,
111 targets=targets,
112 ranges=ranges,
113 mt_type='dc',
114 stf_type='HalfSinusoidSTF',
115 norm_exponent=self.norm_exponent)
117 problem = DynamicRuptureProblem(
118 name=expand_template(self.name_template, subs),
119 base_source=base_source,
120 distance_min=self.distance_min,
121 target_groups=target_groups,
122 targets=targets,
123 ranges=self.ranges,
124 norm_exponent=self.norm_exponent,
125 adaptive_resolution=self.adaptive_resolution,
126 adaptive_start=self.adaptive_start,
127 cmt_problem=cmt_problem)
129 return problem
132@has_get_plot_classes
133class DynamicRuptureProblem(Problem):
135 problem_parameters = [
136 Parameter('east_shift', 'm', label='Easting', **as_km),
137 Parameter('north_shift', 'm', label='Northing', **as_km),
138 Parameter('depth', 'm', label='Depth', **as_km),
139 Parameter('length', 'm', label='Length', **as_km),
140 Parameter('width', 'm', label='Width', **as_km),
141 Parameter('slip', 'm', label='Peak slip', optional=True),
142 Parameter('strike', 'deg', label='Strike'),
143 Parameter('dip', 'deg', label='Dip'),
144 Parameter('rake', 'deg', label='Rake'),
145 Parameter('gamma', 'vr/vs', label=r'$\gamma$'),
146 Parameter('nx', 'patches', label='nx'),
147 Parameter('ny', 'patches', label='ny')
148 ]
150 problem_waveform_parameters = [
151 Parameter('nucleation_x', 'offset', label='Nucleation X'),
152 Parameter('nucleation_y', 'offset', label='Nucleation Y'),
153 Parameter('time', 's', label='Time'),
154 ]
156 dependants = [
157 Parameter('traction', 'Pa', label='Maximum traction'),
158 Parameter('magnitude', label='Magnitude'),
159 ]
161 distance_min = Float.T(default=0.0)
162 adaptive_resolution = String.T()
163 adaptive_start = Int.T()
165 cmt_problem = Problem.T(optional=True)
167 def __init__(self, **kwargs):
168 Problem.__init__(self, **kwargs)
169 self.deps_cache = {}
170 self.minmax_cache = {}
172 def set_engine(self, engine):
173 self._engine = engine
175 if self.cmt_problem is not None:
176 self.cmt_problem.set_engine(engine)
178 def pack(self, source):
179 arr = self.get_parameter_array(source)
180 for ip, p in enumerate(self.parameters):
181 if p.name == 'time':
182 arr[ip] -= self.base_source.time
183 return arr
185 def get_source(self, x):
186 src = self.base_source
188 d = self.get_parameter_dict(x)
189 p = {k: float(self.ranges[k].make_relative(src[k], d[k]))
190 for k in src.keys() if k in d}
192 p['nx'] = int(p['nx'])
193 p['ny'] = int(p['ny'])
195 if p['nx'] != src.nx or p['ny'] != src.ny:
196 logger.info('refining patches to %dx%d', p['nx'], p['ny'])
198 src.nx = p['nx']
199 src.ny = p['ny']
201 return src.clone(**p)
203 def random_uniform(self, xbounds, rstate, fixed_magnitude=None):
204 if fixed_magnitude is not None:
205 raise GrondError(
206 'Setting fixed magnitude in random model generation not '
207 'supported for this type of problem.')
209 x = num.zeros(self.nparameters)
210 for i in range(self.nparameters):
211 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
213 return x
215 def preconstrain(self, x, optimiser=None):
216 nx = self.ranges['nx']
217 ny = self.ranges['ny']
219 if optimiser and self.adaptive_resolution == 'linear' and \
220 optimiser.iiter >= self.adaptive_start:
222 progress = (optimiser.iiter - self.adaptive_start - 1) / \
223 (optimiser.niterations - self.adaptive_start)
225 nx = num.floor(
226 nx.start + progress * (nx.stop - nx.start + 1))
227 ny = num.floor(
228 ny.start + progress * (ny.stop - ny.start + 1))
230 elif optimiser and self.adaptive_resolution == 'uniqueness' and \
231 optimiser.iiter >= self.adaptive_start:
232 raise NotImplementedError
234 else:
235 nx = nx.start
236 ny = ny.start
238 idx_nx = self.get_parameter_index('nx')
239 idx_ny = self.get_parameter_index('ny')
240 x[idx_nx] = int(round(nx))
241 x[idx_ny] = int(round(ny))
243 return x
245 def make_dependant(self, xs, pname, store=None):
246 cache = self.deps_cache
247 if xs.ndim == 1:
248 return self.make_dependant(xs[num.newaxis, :], pname)[0]
250 if pname not in self.dependant_names:
251 raise KeyError(pname)
253 if store is None:
254 store_ids = self.get_gf_store_ids()
255 engine = self.get_engine()
256 store = engine.get_store(store_ids[0])
258 y = num.zeros(xs.shape[0])
259 for i, x in enumerate(xs):
260 k = tuple(x.tolist())
261 if k not in cache:
262 source = self.get_source(x)
264 if num.isnan(source.slip):
265 source.slip = 0.
267 source.discretize_basesource(store)
269 tractions = source.get_scaled_tractions(store)
270 traction = max(num.linalg.norm(tractions, axis=1))
272 magnitude = source.get_magnitude(store=store)
274 cache[k] = dict(
275 traction=traction,
276 magnitude=magnitude)
278 res = cache[k]
280 y[i] = res[pname]
282 return y
284 def get_dependant_bounds(self):
285 cache = self.minmax_cache
287 x_max_traction = self.get_min_model()
289 for i, p in enumerate(self.problem_parameters):
290 if p.name != 'slip':
291 continue
293 x_max_traction[i] = self.ranges[p.name].stop
295 xs = (self.get_min_model(), self.get_max_model(), x_max_traction)
297 if not any(tuple(x.tolist()) in cache for x in xs):
298 store_ids = self.get_gf_store_ids()
299 engine = self.get_engine()
300 store = engine.get_store(store_ids[0])
302 for x in xs:
303 source = self.get_source(x)
304 source.discretize_basesource(store)
306 tractions = source.get_scaled_tractions(store)
307 traction = max(num.linalg.norm(tractions, axis=1))
309 magnitude = source.get_magnitude(store=store)
311 cache[tuple(x.tolist())] = dict(
312 traction=traction,
313 magnitude=magnitude
314 )
316 traction = [cache[tuple(x.tolist())]['traction'] for x in xs]
317 magnitude = [cache[tuple(x.tolist())]['magnitude'] for x in xs]
319 out = [
320 (min(traction), max(traction)),
321 (min(magnitude), max(magnitude))]
323 return out
325 @classmethod
326 def get_plot_classes(cls):
327 from . import plot
328 plots = super(DynamicRuptureProblem, cls).get_plot_classes()
329 plots.extend([
330 plot.DynamicRuptureSlipMap,
331 plot.DynamicRuptureSTF,
332 plot.DynamicRuptureMap])
333 return plots
336class DoublePDRProblemConfig(ProblemConfig):
338 ranges = Dict.T(String.T(), gf.Range.T())
340 decimation_factor = Int.T(
341 default=1,
342 help='Decimation factor of the discretized sub-faults')
343 distance_min = Float.T(default=0.)
344 adaptive_resolution = StringChoice.T(
345 choices=('off', 'linear', 'uniqueness'),
346 default='off')
347 adaptive_start = Int.T(
348 default=0)
350 pure_shear = Bool.T(
351 default=True)
352 smooth_rupture = Bool.T(
353 default=False)
355 point_source_target_balancing = Bool.T(
356 default=False,
357 help='If ``True``, target balancing (if used) is performed on a '
358 'moment tensor point source at the events location. It increases '
359 'the speed, but might lead to not fully optimized target weights.'
360 )
362 def get_problem(self, event, target_groups, targets):
363 if self.decimation_factor != 1:
364 logger.warn(
365 'Decimation factor for dynamic rupture source set to %i.'
366 ' Results may be inaccurate.' % self.decimation_factor)
368 base_source = gf.DoublePDR.from_pyrocko_event(
369 event,
370 nx1=int(self.ranges['nx1'].start),
371 ny1=int(self.ranges['ny1'].start),
372 nx2=int(self.ranges['nx2'].start),
373 ny2=int(self.ranges['ny2'].start),
374 decimation_factor=self.decimation_factor,
375 pure_shear=self.pure_shear,
376 smooth_rupture=self.smooth_rupture,
377 stf=None)
379 subs = dict(
380 event_name=event.name,
381 event_time=util.time_to_str(event.time))
383 cmt_problem = None
384 if self.point_source_target_balancing:
385 base_source_cmt = gf.MTSource.from_pyrocko_event(event)
387 stf = gf.HalfSinusoidSTF()
388 stf.duration = event.duration or 0.0
390 base_source_cmt.stf = stf
392 ranges = dict(
393 time=self.ranges['time'],
394 north_shift=self.ranges['north_shift'],
395 east_shift=self.ranges['east_shift'],
396 depth=self.ranges['depth'],
397 magnitude=gf.Range(
398 start=event.magnitude - 1.,
399 stop=event.magnitude + 1.),
400 duration=gf.Range(start=0., stop=stf.duration * 2.),
401 rmnn=gf.Range(start=-1.41421, stop=1.41421),
402 rmee=gf.Range(start=-1.41421, stop=1.41421),
403 rmdd=gf.Range(start=-1.41421, stop=1.41421),
404 rmne=gf.Range(start=-1., stop=1.),
405 rmnd=gf.Range(start=-1., stop=1.),
406 rmed=gf.Range(start=-1., stop=1.))
408 cmt_problem = CMTProblem(
409 name=expand_template(self.name_template, subs),
410 base_source=base_source_cmt,
411 distance_min=self.distance_min,
412 target_groups=target_groups,
413 targets=targets,
414 ranges=ranges,
415 mt_type='dc',
416 stf_type='HalfSinusoidSTF',
417 norm_exponent=self.norm_exponent)
419 problem = DoublePDRProblem(
420 name=expand_template(self.name_template, subs),
421 base_source=base_source,
422 distance_min=self.distance_min,
423 target_groups=target_groups,
424 targets=targets,
425 ranges=self.ranges,
426 norm_exponent=self.norm_exponent,
427 adaptive_resolution=self.adaptive_resolution,
428 adaptive_start=self.adaptive_start,
429 cmt_problem=cmt_problem)
431 return problem
434@has_get_plot_classes
435class DoublePDRProblem(Problem):
437 problem_parameters = [
438 Parameter('east_shift', 'm', label='Easting', **as_km),
439 Parameter('north_shift', 'm', label='Northing', **as_km),
440 Parameter('depth', 'm', label='Depth', **as_km),
441 Parameter('length1', 'm', label='Length 1', **as_km),
442 Parameter('width1', 'm', label='Width 1', **as_km),
443 Parameter('slip1', 'm', label='Peak slip 1', optional=True),
444 Parameter('strike1', 'deg', label='Strike 1'),
445 Parameter('dip1', 'deg', label='Dip 1'),
446 Parameter('rake1', 'deg', label='Rake 1'),
447 Parameter('gamma1', 'vr/vs', label=r'$\gamma$ 1'),
448 Parameter('nx1', 'patches', label='nx 1'),
449 Parameter('ny1', 'patches', label='ny 1'),
450 Parameter('length2', 'm', label='Length 2', **as_km),
451 Parameter('width2', 'm', label='Width 2', **as_km),
452 Parameter('slip2', 'm', label='Peak slip 2', optional=True),
453 Parameter('strike2', 'deg', label='Strike 2'),
454 Parameter('dip2', 'deg', label='Dip 2'),
455 Parameter('rake2', 'deg', label='Rake 2'),
456 Parameter('gamma2', 'vr/vs', label=r'$\gamma$ 2'),
457 Parameter('nx2', 'patches', label='nx 2'),
458 Parameter('ny2', 'patches', label='ny 2'),
459 Parameter('delta_time', 's', label='$\\Delta$ Time'),
460 Parameter('delta_depth', 'm', label='$\\Delta$ Depth'),
461 Parameter('azimuth', 'deg', label='Azimuth'),
462 Parameter('distance', 'm', label='Distance'),
463 Parameter('mix', label='Mix')
464 ]
466 problem_waveform_parameters = [
467 Parameter('nucleation_x1', 'offset', label='Nucleation X 1'),
468 Parameter('nucleation_y1', 'offset', label='Nucleation Y 1'),
469 Parameter('nucleation_x2', 'offset', label='Nucleation X 2'),
470 Parameter('nucleation_y2', 'offset', label='Nucleation Y 2'),
471 Parameter('time', 's', label='Time'),
472 ]
474 dependants = []
476 distance_min = Float.T(default=0.0)
477 adaptive_resolution = String.T()
478 adaptive_start = Int.T()
480 cmt_problem = Problem.T(optional=True)
482 def set_engine(self, engine):
483 self._engine = engine
485 if self.cmt_problem is not None:
486 self.cmt_problem.set_engine(engine)
488 def pack(self, source):
489 arr = self.get_parameter_array(source)
490 for ip, p in enumerate(self.parameters):
491 if p.name == 'time':
492 arr[ip] -= self.base_source.time
493 return arr
495 def get_source(self, x):
496 src = self.base_source
498 d = self.get_parameter_dict(x)
499 p = {k: float(self.ranges[k].make_relative(src[k], d[k]))
500 for k in src.keys() if k in d}
502 p['nx1'] = int(p['nx1'])
503 p['ny1'] = int(p['ny1'])
505 if p['nx1'] != src.nx1 or p['ny1'] != src.ny1:
506 logger.info('refining patches to %dx%d', p['nx1'], p['ny1'])
508 src.nx1 = p['nx1']
509 src.ny1 = p['ny1']
511 p['nx2'] = int(p['nx2'])
512 p['ny2'] = int(p['ny2'])
514 if p['nx2'] != src.nx2 or p['ny2'] != src.ny2:
515 logger.info('refining patches to %dx%d', p['nx2'], p['ny2'])
517 src.nx2 = p['nx2']
518 src.ny2 = p['ny2']
520 return src.clone(**p)
522 def random_uniform(self, xbounds, rstate, fixed_magnitude=None):
523 if fixed_magnitude is not None:
524 raise GrondError(
525 'Setting fixed magnitude in random model generation not '
526 'supported for this type of problem.')
528 x = num.zeros(self.nparameters)
529 for i in range(self.nparameters):
530 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
532 return x
534 def preconstrain(self, x, optimiser=None):
535 nx1 = self.ranges['nx1']
536 ny1 = self.ranges['ny1']
537 nx2 = self.ranges['nx2']
538 ny2 = self.ranges['ny2']
540 if optimiser and self.adaptive_resolution == 'linear' and \
541 optimiser.iiter >= self.adaptive_start:
543 progress = (optimiser.iiter - self.adaptive_start - 1) / \
544 (optimiser.niterations - self.adaptive_start)
546 nx1 = num.floor(
547 nx1.start + progress * (nx1.stop - nx1.start + 1))
548 ny1 = num.floor(
549 ny1.start + progress * (ny1.stop - ny1.start + 1))
550 nx2 = num.floor(
551 nx2.start + progress * (nx2.stop - nx2.start + 1))
552 ny2 = num.floor(
553 ny2.start + progress * (ny2.stop - ny2.start + 1))
555 elif optimiser and self.adaptive_resolution == 'uniqueness' and \
556 optimiser.iiter >= self.adaptive_start:
557 raise NotImplementedError
559 else:
560 nx1 = nx1.start
561 ny1 = ny1.start
562 nx2 = nx2.start
563 ny2 = ny2.start
565 idx_nx1 = self.get_parameter_index('nx1')
566 idx_ny1 = self.get_parameter_index('ny1')
567 idx_nx2 = self.get_parameter_index('nx2')
568 idx_ny2 = self.get_parameter_index('ny2')
569 x[idx_nx1] = int(round(nx1))
570 x[idx_ny1] = int(round(ny1))
571 x[idx_nx2] = int(round(nx2))
572 x[idx_ny2] = int(round(ny2))
574 return x
576 @classmethod
577 def get_plot_classes(cls):
578 plots = super(DoublePDRProblem, cls).get_plot_classes()
580 return plots
583class TriplePDRProblemConfig(ProblemConfig):
585 ranges = Dict.T(String.T(), gf.Range.T())
587 decimation_factor = Int.T(
588 default=1,
589 help='Decimation factor of the discretized sub-faults')
590 distance_min = Float.T(default=0.)
591 adaptive_resolution = StringChoice.T(
592 choices=('off', 'linear', 'uniqueness'),
593 default='off')
594 adaptive_start = Int.T(
595 default=0)
597 pure_shear = Bool.T(
598 default=True)
599 smooth_rupture = Bool.T(
600 default=False)
602 def get_problem(self, event, target_groups, targets):
603 if self.decimation_factor != 1:
604 logger.warn(
605 'Decimation factor for dynamic rupture source set to %i.'
606 ' Results may be inaccurate.' % self.decimation_factor)
608 base_source = gf.TriplePDR.from_pyrocko_event(
609 event,
610 nx1=int(self.ranges['nx1'].start),
611 ny1=int(self.ranges['ny1'].start),
612 nx2=int(self.ranges['nx2'].start),
613 ny2=int(self.ranges['ny2'].start),
614 nx3=int(self.ranges['nx3'].start),
615 ny3=int(self.ranges['ny3'].start),
616 decimation_factor=self.decimation_factor,
617 pure_shear=self.pure_shear,
618 smooth_rupture=self.smooth_rupture,
619 stf=None)
621 subs = dict(
622 event_name=event.name,
623 event_time=util.time_to_str(event.time))
625 cmt_problem = None
627 problem = TriplePDRProblem(
628 name=expand_template(self.name_template, subs),
629 base_source=base_source,
630 distance_min=self.distance_min,
631 target_groups=target_groups,
632 targets=targets,
633 ranges=self.ranges,
634 norm_exponent=self.norm_exponent,
635 adaptive_resolution=self.adaptive_resolution,
636 adaptive_start=self.adaptive_start,
637 cmt_problem=cmt_problem)
639 return problem
642@has_get_plot_classes
643class TriplePDRProblem(Problem):
645 problem_parameters = [
646 Parameter('delta_time1', 's', label=r'$\Delta$Time 1'),
647 Parameter('north_shift1', 'm', label='Northing 1', **as_km),
648 Parameter('east_shift1', 'm', label='Easting 1', **as_km),
649 Parameter('depth1', 'm', label='Depth 1', **as_km),
650 Parameter('length1', 'm', label='Length 1', **as_km),
651 Parameter('width1', 'm', label='Width 1', **as_km),
652 Parameter('slip1', 'm', label='Peak slip 1', optional=True),
653 Parameter('strike1', 'deg', label='Strike 1'),
654 Parameter('dip1', 'deg', label='Dip 1'),
655 Parameter('rake1', 'deg', label='Rake 1'),
656 Parameter('gamma1', 'vr/vs', label=r'$\gamma$ 1'),
657 Parameter('nx1', 'patches', label='nx 1'),
658 Parameter('ny1', 'patches', label='ny 1'),
659 Parameter('delta_time2', 's', label=r'$\Delta$Time 2'),
660 Parameter('north_shift2', 'm', label='Northing 2', **as_km),
661 Parameter('east_shift2', 'm', label='Easting 2', **as_km),
662 Parameter('depth2', 'm', label='Depth 2', **as_km),
663 Parameter('length2', 'm', label='Length 2', **as_km),
664 Parameter('width2', 'm', label='Width 2', **as_km),
665 Parameter('slip2', 'm', label='Peak slip 2', optional=True),
666 Parameter('strike2', 'deg', label='Strike 2'),
667 Parameter('dip2', 'deg', label='Dip 2'),
668 Parameter('rake2', 'deg', label='Rake 2'),
669 Parameter('gamma2', 'vr/vs', label=r'$\gamma$ 2'),
670 Parameter('nx2', 'patches', label='nx 2'),
671 Parameter('ny2', 'patches', label='ny 2'),
672 Parameter('delta_time3', 's', label=r'$\Delta$Time 3'),
673 Parameter('north_shift3', 'm', label='Northing 3', **as_km),
674 Parameter('east_shift3', 'm', label='Easting 3', **as_km),
675 Parameter('depth3', 'm', label='Depth 3', **as_km),
676 Parameter('length3', 'm', label='Length 3', **as_km),
677 Parameter('width3', 'm', label='Width 3', **as_km),
678 Parameter('slip3', 'm', label='Peak slip 3', optional=True),
679 Parameter('strike3', 'deg', label='Strike 3'),
680 Parameter('dip3', 'deg', label='Dip 3'),
681 Parameter('rake3', 'deg', label='Rake 3'),
682 Parameter('gamma3', 'vr/vs', label=r'$\gamma$ 3'),
683 Parameter('nx3', 'patches', label='nx 3'),
684 Parameter('ny3', 'patches', label='ny 3')
685 ]
687 problem_waveform_parameters = [
688 Parameter('nucleation_x1', 'offset', label='Nucleation X 1'),
689 Parameter('nucleation_y1', 'offset', label='Nucleation Y 1'),
690 Parameter('nucleation_x2', 'offset', label='Nucleation X 2'),
691 Parameter('nucleation_y2', 'offset', label='Nucleation Y 2'),
692 Parameter('nucleation_x3', 'offset', label='Nucleation X 3'),
693 Parameter('nucleation_y3', 'offset', label='Nucleation Y 3')
694 ]
696 dependants = []
698 distance_min = Float.T(default=0.0)
699 adaptive_resolution = String.T()
700 adaptive_start = Int.T()
702 cmt_problem = Problem.T(optional=True)
704 def set_engine(self, engine):
705 self._engine = engine
707 if self.cmt_problem is not None:
708 self.cmt_problem.set_engine(engine)
710 def pack(self, source):
711 return self.get_parameter_array(source)
712 # return arr
714 def get_source(self, x):
715 src = self.base_source
717 d = self.get_parameter_dict(x)
718 p = {k: float(self.ranges[k].make_relative(src[k], d[k]))
719 for k in src.keys() if k in d}
721 p['nx1'] = int(p['nx1'])
722 p['ny1'] = int(p['ny1'])
724 if p['nx1'] != src.nx1 or p['ny1'] != src.ny1:
725 logger.info('refining patches to %dx%d', p['nx1'], p['ny1'])
727 src.nx1 = p['nx1']
728 src.ny1 = p['ny1']
730 p['nx2'] = int(p['nx2'])
731 p['ny2'] = int(p['ny2'])
733 if p['nx2'] != src.nx2 or p['ny2'] != src.ny2:
734 logger.info('refining patches to %dx%d', p['nx2'], p['ny2'])
736 src.nx2 = p['nx2']
737 src.ny2 = p['ny2']
739 p['nx3'] = int(p['nx3'])
740 p['ny3'] = int(p['ny3'])
742 if p['nx3'] != src.nx3 or p['ny3'] != src.ny3:
743 logger.info('refining patches to %dx%d', p['nx3'], p['ny3'])
745 src.nx3 = p['nx3']
746 src.ny3 = p['ny3']
748 return src.clone(**p)
750 def random_uniform(self, xbounds, rstate, fixed_magnitude=None):
751 if fixed_magnitude is not None:
752 raise GrondError(
753 'Setting fixed magnitude in random model generation not '
754 'supported for this type of problem.')
756 x = num.zeros(self.nparameters)
757 for i in range(self.nparameters):
758 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
760 return x
762 def preconstrain(self, x, optimiser=None):
763 nx1 = self.ranges['nx1']
764 ny1 = self.ranges['ny1']
765 nx2 = self.ranges['nx2']
766 ny2 = self.ranges['ny2']
767 nx3 = self.ranges['nx3']
768 ny3 = self.ranges['ny3']
770 if optimiser and self.adaptive_resolution == 'linear' and \
771 optimiser.iiter >= self.adaptive_start:
773 progress = (optimiser.iiter - self.adaptive_start - 1) / \
774 (optimiser.niterations - self.adaptive_start)
776 nx1 = num.floor(
777 nx1.start + progress * (nx1.stop - nx1.start + 1))
778 ny1 = num.floor(
779 ny1.start + progress * (ny1.stop - ny1.start + 1))
780 nx2 = num.floor(
781 nx2.start + progress * (nx2.stop - nx2.start + 1))
782 ny2 = num.floor(
783 ny2.start + progress * (ny2.stop - ny2.start + 1))
784 nx3 = num.floor(
785 nx3.start + progress * (nx3.stop - nx3.start + 1))
786 ny3 = num.floor(
787 ny3.start + progress * (ny3.stop - ny3.start + 1))
789 elif optimiser and self.adaptive_resolution == 'uniqueness' and \
790 optimiser.iiter >= self.adaptive_start:
791 raise NotImplementedError
793 else:
794 nx1 = nx1.start
795 ny1 = ny1.start
796 nx2 = nx2.start
797 ny2 = ny2.start
798 nx3 = nx3.start
799 ny3 = ny3.start
801 idx_nx1 = self.get_parameter_index('nx1')
802 idx_ny1 = self.get_parameter_index('ny1')
803 idx_nx2 = self.get_parameter_index('nx2')
804 idx_ny2 = self.get_parameter_index('ny2')
805 idx_nx3 = self.get_parameter_index('nx3')
806 idx_ny3 = self.get_parameter_index('ny3')
807 x[idx_nx1] = int(round(nx1))
808 x[idx_ny1] = int(round(ny1))
809 x[idx_nx2] = int(round(nx2))
810 x[idx_ny2] = int(round(ny2))
811 x[idx_nx3] = int(round(nx3))
812 x[idx_ny3] = int(round(ny3))
814 return x
816 @classmethod
817 def get_plot_classes(cls):
818 plots = super(TriplePDRProblem, cls).get_plot_classes()
820 return plots
823__all__ = '''
824 DynamicRuptureProblem
825 DynamicRuptureProblemConfig
826 DoublePDRProblem
827 DoublePDRProblemConfig
828 TriplePDRProblem
829 TriplePDRProblemConfig
830'''.split()