|
41 | 41 | "cell_type": "markdown", |
42 | 42 | "metadata": {}, |
43 | 43 | "source": [ |
44 | | - "We will need to use subdomains to accomodate the modified equations in the PML regions. As `Grid` objects cannot have subdomains added retroactively, we must define our subdomains beforehand." |
| 44 | + "We will need to use subdomains to accommodate the modified equations in the PML regions." |
45 | 45 | ] |
46 | 46 | }, |
47 | 47 | { |
|
52 | 52 | "source": [ |
53 | 53 | "from devito import SubDomain\n", |
54 | 54 | "\n", |
55 | | - "s_o = 6 # Space order\n", |
| 55 | + "so = 6 # Space order\n", |
56 | 56 | "\n", |
57 | 57 | "class MainDomain(SubDomain): # Main section with no damping\n", |
58 | 58 | " name = 'main'\n", |
59 | | - " def __init__(self, PMLS, S_O):\n", |
60 | | - " super().__init__()\n", |
61 | | - " self.PMLS = PMLS\n", |
62 | | - " self.S_O = S_O\n", |
| 59 | + " def __init__(self, pmls, so, grid=None):\n", |
| 60 | + " # NOTE: These attributes are used in `define`, and thus must be\n", |
| 61 | + " # set up before `super().__init__` is called.\n", |
| 62 | + " self.pmls = pmls\n", |
| 63 | + " self.so = so\n", |
| 64 | + " super().__init__(grid=grid)\n", |
63 | 65 | "\n", |
64 | 66 | " def define(self, dimensions):\n", |
65 | 67 | " x, y = dimensions\n", |
66 | | - " return {x: ('middle', self.PMLS, self.PMLS),\n", |
67 | | - " y: ('middle', self.S_O//2, self.PMLS)}\n", |
| 68 | + " return {x: ('middle', self.pmls, self.pmls),\n", |
| 69 | + " y: ('middle', self.so//2, self.pmls)}\n", |
68 | 70 | "\n", |
69 | 71 | "\n", |
70 | 72 | "class Left(SubDomain): # Left PML region\n", |
71 | 73 | " name = 'left'\n", |
72 | | - " def __init__(self, PMLS):\n", |
73 | | - " super().__init__()\n", |
74 | | - " self.PMLS = PMLS\n", |
| 74 | + " def __init__(self, pmls, grid=None):\n", |
| 75 | + " self.pmls = pmls\n", |
| 76 | + " super().__init__(grid=grid)\n", |
75 | 77 | "\n", |
76 | 78 | " def define(self, dimensions):\n", |
77 | 79 | " x, y = dimensions\n", |
78 | | - " return {x: ('left', self.PMLS), y: y}\n", |
| 80 | + " return {x: ('left', self.pmls), y: y}\n", |
79 | 81 | "\n", |
80 | 82 | "\n", |
81 | 83 | "class Right(SubDomain): # Right PML region\n", |
82 | 84 | " name = 'right'\n", |
83 | | - " def __init__(self, PMLS):\n", |
84 | | - " super().__init__()\n", |
85 | | - " self.PMLS = PMLS\n", |
| 85 | + " def __init__(self, pmls, grid=None):\n", |
| 86 | + " self.pmls = pmls\n", |
| 87 | + " super().__init__(grid=grid)\n", |
86 | 88 | "\n", |
87 | 89 | " def define(self, dimensions):\n", |
88 | 90 | " x, y = dimensions\n", |
89 | | - " return {x: ('right', self.PMLS), y: y}\n", |
| 91 | + " return {x: ('right', self.pmls), y: y}\n", |
90 | 92 | "\n", |
91 | 93 | "class Base(SubDomain): # Base PML region\n", |
92 | 94 | " name = 'base'\n", |
93 | | - " def __init__(self, PMLS):\n", |
94 | | - " super().__init__()\n", |
95 | | - " self.PMLS = PMLS\n", |
| 95 | + " def __init__(self, pmls, grid=None):\n", |
| 96 | + " self.pmls = pmls\n", |
| 97 | + " super().__init__(grid=grid)\n", |
96 | 98 | "\n", |
97 | 99 | " def define(self, dimensions):\n", |
98 | 100 | " x, y = dimensions\n", |
99 | | - " return {x: ('middle', self.PMLS, self.PMLS), y: ('right', self.PMLS)}\n", |
| 101 | + " return {x: ('middle', self.pmls, self.pmls), y: ('right', self.pmls)}\n", |
100 | 102 | "\n", |
101 | 103 | "class FreeSurface(SubDomain): # Free surface region\n", |
102 | 104 | " name = 'freesurface'\n", |
103 | | - " def __init__(self, PMLS, S_O):\n", |
104 | | - " super().__init__()\n", |
105 | | - " self.PMLS = PMLS\n", |
106 | | - " self.S_O = S_O\n", |
| 105 | + " def __init__(self, pmls, so, grid=None):\n", |
| 106 | + " self.pmls = pmls\n", |
| 107 | + " self.so = so\n", |
| 108 | + " super().__init__(grid=grid)\n", |
107 | 109 | "\n", |
108 | 110 | " def define(self, dimensions):\n", |
109 | 111 | " x, y = dimensions\n", |
110 | | - " return {x: ('middle', self.PMLS, self.PMLS), y: ('left', self.S_O//2)}\n", |
111 | | - "\n", |
112 | | - "main_domain = MainDomain(nbpml, s_o)\n", |
113 | | - "left = Left(nbpml)\n", |
114 | | - "right = Right(nbpml)\n", |
115 | | - "base = Base(nbpml)\n", |
116 | | - "freesurface = FreeSurface(nbpml, s_o)" |
| 112 | + " return {x: ('middle', self.pmls, self.pmls), y: ('left', self.so//2)}" |
117 | 113 | ] |
118 | 114 | }, |
119 | 115 | { |
120 | 116 | "cell_type": "markdown", |
121 | 117 | "metadata": {}, |
122 | 118 | "source": [ |
123 | | - "And create the grid, attaching our subdomains:" |
| 119 | + "We create the grid and set up our subdomains:" |
124 | 120 | ] |
125 | 121 | }, |
126 | 122 | { |
|
131 | 127 | "source": [ |
132 | 128 | "from devito import Grid\n", |
133 | 129 | "\n", |
134 | | - "grid = Grid(shape=shape, extent=extent,\n", |
135 | | - " subdomains=(main_domain, left, right, base, freesurface))\n", |
| 130 | + "grid = Grid(shape=shape, extent=extent)\n", |
| 131 | + "\n", |
| 132 | + "main = MainDomain(nbpml, so, grid=grid)\n", |
| 133 | + "left = Left(nbpml, grid=grid)\n", |
| 134 | + "right = Right(nbpml, grid=grid)\n", |
| 135 | + "base = Base(nbpml, grid=grid)\n", |
| 136 | + "freesurface = FreeSurface(nbpml, so, grid=grid)\n", |
| 137 | + "\n", |
136 | 138 | "x, y = grid.dimensions" |
137 | 139 | ] |
138 | 140 | }, |
|
170 | 172 | "from devito import TimeFunction, VectorTimeFunction, NODE\n", |
171 | 173 | "\n", |
172 | 174 | "p = TimeFunction(name='p', grid=grid, time_order=1,\n", |
173 | | - " space_order=s_o, staggered=NODE)\n", |
| 175 | + " space_order=so, staggered=NODE)\n", |
174 | 176 | "v = VectorTimeFunction(name='v', grid=grid, time_order=1,\n", |
175 | | - " space_order=s_o)" |
| 177 | + " space_order=so)" |
176 | 178 | ] |
177 | 179 | }, |
178 | 180 | { |
|
323 | 325 | "from devito import Eq, grad, div\n", |
324 | 326 | "\n", |
325 | 327 | "eq_v = Eq(v.forward, v - dt*grad(p)/density,\n", |
326 | | - " subdomain=grid.subdomains['main'])\n", |
| 328 | + " subdomain=main)\n", |
327 | 329 | "\n", |
328 | 330 | "eq_p = Eq(p.forward, p - dt*velocity**2*density*div(v.forward),\n", |
329 | | - " subdomain=grid.subdomains['main'])" |
| 331 | + " subdomain=main)" |
330 | 332 | ] |
331 | 333 | }, |
332 | 334 | { |
|
361 | 363 | "# Left side\n", |
362 | 364 | "eq_v_damp_left = Eq(v.forward,\n", |
363 | 365 | " (1-d_l)*v - dt*grad(p)/density,\n", |
364 | | - " subdomain=grid.subdomains['left'])\n", |
| 366 | + " subdomain=left)\n", |
365 | 367 | "\n", |
366 | 368 | "eq_p_damp_left = Eq(p.forward,\n", |
367 | 369 | " (1-gamma*velocity**2*dt-d_l*dt)*p\n", |
368 | 370 | " - d_l*gamma*velocity**2*p_i\n", |
369 | 371 | " - dt*velocity**2*density*div(v.forward),\n", |
370 | | - " subdomain=grid.subdomains['left'])\n", |
| 372 | + " subdomain=left)\n", |
371 | 373 | "\n", |
372 | 374 | "# Right side\n", |
373 | 375 | "eq_v_damp_right = Eq(v.forward,\n", |
374 | 376 | " (1-d_r)*v - dt*grad(p)/density,\n", |
375 | | - " subdomain=grid.subdomains['right'])\n", |
| 377 | + " subdomain=right)\n", |
376 | 378 | "\n", |
377 | 379 | "eq_p_damp_right = Eq(p.forward,\n", |
378 | 380 | " (1-gamma*velocity**2*dt-d_r*dt)*p\n", |
379 | 381 | " - d_r*gamma*velocity**2*p_i\n", |
380 | 382 | " - dt*velocity**2*density*div(v.forward),\n", |
381 | | - " subdomain=grid.subdomains['right'])\n", |
| 383 | + " subdomain=right)\n", |
382 | 384 | "\n", |
383 | 385 | "# Base edge\n", |
384 | 386 | "eq_v_damp_base = Eq(v.forward,\n", |
385 | 387 | " (1-d_b)*v - dt*grad(p)/density,\n", |
386 | | - " subdomain=grid.subdomains['base'])\n", |
| 388 | + " subdomain=base)\n", |
387 | 389 | "\n", |
388 | 390 | "eq_p_damp_base = Eq(p.forward,\n", |
389 | 391 | " (1-gamma*velocity**2*dt-d_b*dt)*p\n", |
390 | 392 | " - d_b*gamma*velocity**2*p_i\n", |
391 | 393 | " - dt*velocity**2*density*div(v.forward),\n", |
392 | | - " subdomain=grid.subdomains['base'])" |
| 394 | + " subdomain=base)" |
393 | 395 | ] |
394 | 396 | }, |
395 | 397 | { |
|
456 | 458 | " mapper.update({f: f.subs({yind: INT(abs(yind))})})\n", |
457 | 459 | " return Eq(lhs, rhs.subs(mapper), subdomain=subdomain)\n", |
458 | 460 | " \n", |
459 | | - "fs_p = free_surface_top(eq_p, grid.subdomains['freesurface'], 'pressure')\n", |
460 | | - "fs_v = free_surface_top(eq_v, grid.subdomains['freesurface'], 'velocity')" |
| 461 | + "fs_p = free_surface_top(eq_p, freesurface, 'pressure')\n", |
| 462 | + "fs_v = free_surface_top(eq_v, freesurface, 'velocity')" |
461 | 463 | ] |
462 | 464 | }, |
463 | 465 | { |
|
616 | 618 | "name": "python", |
617 | 619 | "nbconvert_exporter": "python", |
618 | 620 | "pygments_lexer": "ipython3", |
619 | | - "version": "3.9.16" |
| 621 | + "version": "3.11.5" |
620 | 622 | } |
621 | 623 | }, |
622 | 624 | "nbformat": 4, |
|
0 commit comments