diff --git a/CHANGELOG.md b/CHANGELOG.md index ceb2320..3493f51 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -65,6 +65,28 @@ and this project adheres to local-window FWI and an x-strip domain decomposition (`SEMModelParallel`). User guide: `docs/user-guide/sem.md`; examples: the SEM section (notebooks 41-44). +- **`BiotVSSEM2DDeformed`: Biot poroelastic waves on the SEM.** The FD `Biot`'s + velocity-stress-pressure form with the Darcy flux on the Quad9 mesh, its eight + coefficients and its `source_type` / `receiver_type` names through `PropSEM`. The + pore pressure is continuous and the flux element-local, so the pore fluid is the + acoustic pressure SEM (with the flux continuous, as in SPECFEM2D, a point source + excites grid-scale modes along its grid lines); the viscous drag integrated + exactly per node, C-PML, an open-pore free surface (`p = 0`) that a `topography` + carries too. Checked against analytic and wavenumber-integration solutions, the + FD records and SPECFEM2D's `poroelastic_semi_infinite_homogeneous` example. 2-D, + eager, no boundary saving. The receiver gather now also reads element-local + fields (the flux and the stresses). With SPECFEM2D's other poroelastic + features: a point force `fx` / `fz` and flux sources `qx` / `qz`; an anisotropic + permeability (`BiotAnisoPermSEM2DDeformed`, the drag a tensor); the viscous + attenuation (`equation_options={"viscous_q": (Q0, f0)}`, a standard linear + solid); fluid (`rho = 0`) and dry elastic (`M = 0`) regions in one model, coupled + as in SPECFEM2D (open pores to a fluid, sealed to a dry solid), sources loading + only the fields of their medium. `sweep.equations.biot_params`: the drag tensor, + the fluid and dry-solid coefficients, and the parameter families of Morency, Luo & + Tromp (2009) with their kernels by the exact chain rule. Checked against + SPECFEM2D runs of each and its `poroelastic_acoustic` example. +- **`PropSEM(equation_options=...)`**: an equation's own `build` keywords; a name the + equation does not take is refused. - **`impl='c'` for `sweep.sem`.** The SEM CUDA kernels are 43 entries of the core (kind 3 of the C API, called through `sweep_sem_call`), their input structs generated from `cuda/equations/sem{2,3}d/*_types.h` like the FD ones. diff --git a/docs/api/index.md b/docs/api/index.md index d0acc79..c754bab 100644 --- a/docs/api/index.md +++ b/docs/api/index.md @@ -228,7 +228,7 @@
SEM equations
2D + 3D
-

Acoustic, elastic (displacement and velocity-stress, VTI) and coupled fluid-solid, on flat and curved meshes.

+

Acoustic, elastic (displacement and velocity-stress, VTI), coupled fluid-solid and poroelastic (Biot), on flat and curved meshes.

registry: equation_classes describe
 diff --git a/docs/api/sem/equations.md b/docs/api/sem/equations.md index 70e8545..06510e4 100644 --- a/docs/api/sem/equations.md +++ b/docs/api/sem/equations.md @@ -159,3 +159,25 @@ each; records `chi`, `u_x`, `u_z` and the derived `p`, `vx`, `vz`. members: [build] inherited_members: true show_if_no_docstring: true + +## Poroelastic (Biot) + +The FD `sweep.equations.Biot` on the spectral-element mesh: solid velocity +and pore pressure on the continuous nodes, Darcy flux and total stress in the +elements, with the same eight coefficients `[H, mu, C, M, rho, rho_f, m, b]` +and the same source and receiver names. See +[Poroelastic waves (Biot)](../../user-guide/sem.md#poroelastic-waves-biot). + +::: sweep.sem.BiotVSSEM2DDeformed + options: + heading_level: 3 + members: [build] + inherited_members: true + show_if_no_docstring: true + +::: sweep.sem.BiotAnisoPermSEM2DDeformed + options: + heading_level: 3 + members: [build] + inherited_members: true + show_if_no_docstring: true diff --git a/docs/user-guide/cli.md b/docs/user-guide/cli.md index aa280c1..efa9f26 100644 --- a/docs/user-guide/cli.md +++ b/docs/user-guide/cli.md @@ -75,22 +75,24 @@ Finite difference (sweep.equations, run with PropTorch): Spectral element (sweep.sem, run with PropSEM): - Equation Mesh Models Torch Binding Binding Ready - ----------------------- --------------------- --------------------------------------- ------------- ------------- - AcousticElasticSEM2D quad9 2-D, topography ['vp', 'vs', 'rho'] yes yes - AcousticPotentialSEM2D quad9 2-D, topography ['vp'] yes yes - AcousticSEM2D quad4 2-D ['vp'] no no - AcousticSEM2DDeformed quad9 2-D, topography ['vp'] no no - AcousticSEM3D hex8 3-D ['vp'] yes yes - AcousticSEM3DDeformed hex27 3-D, topography ['vp'] yes yes - ElasticSEM2D quad4 2-D ['vp', 'vs', 'rho'] no no - ElasticSEM2DDeformed quad9 2-D, topography ['vp', 'vs', 'rho'] yes yes - ElasticSEM3D hex8 3-D ['vp', 'vs', 'rho'] yes yes - ElasticSEM3DDeformed hex27 3-D, topography ['vp', 'vs', 'rho'] yes yes - ElasticVSSEM2DDeformed quad9 2-D, topography ['vp', 'vs', 'rho'] no no - ElasticVSSEM3D hex8 3-D ['vp', 'vs', 'rho'] yes yes - ElasticVSSEM3DDeformed hex27 3-D, topography ['vp', 'vs', 'rho'] no no - ElasticVTISEM2DDeformed quad9 2-D, topography ['vp', 'vs', 'rho', 'epsilon', 'delta'] no no + Equation Mesh Models Torch Binding Binding Ready + -------------------------- --------------------- ------------------------------------------------------------------ ------------- ------------- + AcousticElasticSEM2D quad9 2-D, topography ['vp', 'vs', 'rho'] yes yes + AcousticPotentialSEM2D quad9 2-D, topography ['vp'] yes yes + AcousticSEM2D quad4 2-D ['vp'] no no + AcousticSEM2DDeformed quad9 2-D, topography ['vp'] no no + AcousticSEM3D hex8 3-D ['vp'] yes yes + AcousticSEM3DDeformed hex27 3-D, topography ['vp'] yes yes + BiotAnisoPermSEM2DDeformed quad9 2-D, topography ['H', 'mu', 'C', 'M', 'rho', 'rho_f', 'm', 'b_xx', 'b_xz', 'b_zz'] no no + BiotVSSEM2DDeformed quad9 2-D, topography ['H', 'mu', 'C', 'M', 'rho', 'rho_f', 'm', 'b'] no no + ElasticSEM2D quad4 2-D ['vp', 'vs', 'rho'] no no + ElasticSEM2DDeformed quad9 2-D, topography ['vp', 'vs', 'rho'] yes yes + ElasticSEM3D hex8 3-D ['vp', 'vs', 'rho'] yes yes + ElasticSEM3DDeformed hex27 3-D, topography ['vp', 'vs', 'rho'] yes yes + ElasticVSSEM2DDeformed quad9 2-D, topography ['vp', 'vs', 'rho'] no no + ElasticVSSEM3D hex8 3-D ['vp', 'vs', 'rho'] yes yes + ElasticVSSEM3DDeformed hex27 3-D, topography ['vp', 'vs', 'rho'] no no + ElasticVTISEM2DDeformed quad9 2-D, topography ['vp', 'vs', 'rho', 'epsilon', 'delta'] no no ``` The unified facades `DAS` (formerly `DASModeler`) and `AcousticAniso` are @@ -172,8 +174,16 @@ sweep show AcousticPotentialSEM2D On the elastic SEM classes `source_type` also takes `PropTorch`'s names: `vx`, `vz` (and `vy`) for a force and, in 2-D, `sxx`, `szz`, `sxz` for a moment -tensor (see [SEM · Source types](sem.md#source-types)). A class without an -absorbing layer says so, as `ElasticSEM2D` does: +tensor (see [SEM · Source types](sem.md#source-types)). The Biot classes +name the `source_type` values `PropSEM` takes for them (a force, a moment, a +pore-pressure or flux injection), so `show` prints those instead of their +source fields: + +```text + Source types: ['vx', 'vz', 'qx', 'qz', 'v_x', 'v_z', 'q_x', 'q_z', 'fx', 'fz', 'p', 'sxx', 'szz', 'sxz'] +``` + +A class without an absorbing layer says so, as `ElasticSEM2D` does: ```text C-PML: no -- a reflecting box, abcn=0 only @@ -193,7 +203,7 @@ sweep show Acoustic3d ; echo "exit=$?" ```text No such wave equation: Acoustic3d Did you mean: Acoustic3D, Acoustic, AcousticVRZ3D? - `sweep list equations` names all 43 finite-difference and 14 spectral-element equations. + `sweep list equations` names all 43 finite-difference and 16 spectral-element equations. exit=1 ``` diff --git a/docs/user-guide/equations.md b/docs/user-guide/equations.md index 1a3ba86..9ec2ebb 100644 --- a/docs/user-guide/equations.md +++ b/docs/user-guide/equations.md @@ -86,8 +86,8 @@ The `impl="c"` column says whether `PropSEM(impl='c')` can route the class to the CUDA kernels; every class runs on `impl="eager"`. The `BoundarySaving()` column is the CUDA-side coverage of [boundary saving](sem.md#boundary-saving) -- on `impl="eager"` it works for -every class; ✅ = both backends, "eager" = eager only, and the C-PML is -required wherever the ring IS the C-PML band. +every class but the dissipative Biot one; ✅ = both backends, "eager" = eager +only, and the C-PML is required wherever the ring IS the C-PML band. | Equation | Models | Mesh | Dim | Notes | `impl="c"` | `BoundarySaving()` | | --- | --- | --- | --- | --- | --- | --- | @@ -105,6 +105,8 @@ required wherever the ring IS the C-PML band. | `ElasticVSSEM3DDeformed` | `['vp', 'vs', 'rho']` | hex27 | 3D | 3-D elastic with **both** topography and absorbing boundaries (the velocity–stress family is where 3-D C-PML exists) | ❌ | eager | | `ElasticVTISEM2DDeformed` | `['vp', 'vs', 'rho', 'epsilon', 'delta']` | quad9 | 2D | VTI anisotropic velocity–stress: topography + C-PML | ❌ | eager | | `AcousticElasticSEM2D` | `['vp', 'vs', 'rho']` | quad9 | 2D | Coupled fluid–solid (ocean over seabed): split acoustic + elastic C-PML, records `chi, u_x, u_z` (+ derived `p, vx, vz`) | ✅ | ✅ (C-PML; fp32; +anchors) | +| `BiotVSSEM2DDeformed` | `['H', 'mu', 'C', 'M', 'rho', 'rho_f', 'm', 'b']` | quad9 | 2D | [Biot poroelastic](sem.md#poroelastic-waves-biot) velocity–stress with Darcy flux, the FD `Biot`'s equations and names (pore pressure continuous, flux element-local): topography + C-PML, open-pore surface, fluid (`rho = 0`) and dry elastic (`M = 0`) regions, viscous attenuation (`viscous_q`), records `vx, vz, p, qx, qz, sxx, szz, sxz` | ❌ | ❌ (dissipative) | +| `BiotAnisoPermSEM2DDeformed` | `['H', 'mu', 'C', 'M', 'rho', 'rho_f', 'm', 'b_xx', 'b_xz', 'b_zz']` | quad9 | 2D | `BiotVSSEM2DDeformed` with an [anisotropic permeability](sem.md#anisotropic-permeability-and-viscous-attenuation) (drag tensor) | ❌ | ❌ (dissipative) | For the authoritative list of equations exported by your installation, with their models and compiled-binding support, use the CLI: diff --git a/docs/user-guide/sem.md b/docs/user-guide/sem.md index 2bb3f6c..ae99ef3 100644 --- a/docs/user-guide/sem.md +++ b/docs/user-guide/sem.md @@ -507,6 +507,160 @@ corner): the moment tensor differentiates the basis, whose gradient jumps across the edge, and stamped in one element the explosion came out 1-2% lopsided. A force needs none of this; its stamp has no derivative. +## Poroelastic waves (Biot) + +`BiotVSSEM2DDeformed` solves the equations of the FD `sweep.equations.Biot` — +solid velocity, Darcy flux, total stress and pore pressure, the low-frequency +Biot theory — on the spectral-element mesh, with the same eight coefficients +`[H, mu, C, M, rho, rho_f, m, b]` (build them with +`sweep.equations.biot.biot_from_rock`) and the same source and receiver names: + +```python +from sweep.equations.biot import biot_from_rock +from sweep.sem import BiotVSSEM2DDeformed, PropSEM + +models = biot_from_rock(k_s=36e9, rho_s=2650.0, k_f=2.25e9, rho_f=1000.0, eta=1e-3, + k_d=k_d, mu=mu, phi=phi, kappa=1e-13, tortuosity=2.0) +solver = PropSEM(BiotVSSEM2DDeformed, shape=(nz, nx), dh=2.0, dt=1e-4, + free_surface=True, f_dominant=30.0, + source_type=["sxx", "szz"], receiver_type=["vx", "vz", "p"]) +seis = solver(wavelet, sources, receivers, models=models) # (nshots, nt, nrec, 3) +``` + +| `source_type` | source | the wavelet is | +| --- | --- | --- | +| `["sxx", "szz"]` (default), `["sxz"]` | total stress, as on the elastics | the moment rate | +| `["fx"]`, `["fz"]` | a point force on the solid and the fluid alike (SPECFEM2D's poroelastic force) | the force | +| `["p"]` | pore pressure; the total stress does not move | its rate | +| `["vx"]`, `["vz"]` (`v_x`, `v_z`) | solid velocity; the flux does not move | the velocity source | +| `["qx"]`, `["qz"]` (`q_x`, `q_z`) | Darcy flux; the solid does not move | the flux source | + +As on `PropTorch`, a stress source is a moment tensor that `PropSEM` drives +with the wavelet's integral, while the others take the wavelet itself; the two +kinds do not mix in one call (run two and add the records). The force loads +the total momentum and the relative flow by the same amount, which moves the +solid and the fluid in proportion; the velocity and flux sources are injected +into their fields after the update, so they do not mix with a force either. +Receivers take `vx`, `vz`, `p`, `qx`, `qz`, `sxx`, `szz` and `sxz` (default +`["vx", "vz"]`); the flux and the stresses are read from the receiver's +element. + +The solid is discretised as the velocity-stress elastic equation is: +velocity on the continuous nodes, total stress in the elements, momentum in +weak form with a lumped mass. The pore fluid is the pressure form of the +acoustic SEM: `p` lives on the continuous nodes with its equation in weak +form, and the Darcy flux in the elements, driven by the gradient of `p`. On +its own (`rho_f = C = 0`) the fluid is exactly the acoustic pressure SEM. +The other placement -- the flux continuous and `p` in the elements, as in +SPECFEM2D's poroelastic solver -- leaves the relative flow, which has no +shear stiffness, on a continuous field. A point source then excites +grid-scale modes along the grid lines through it. On those lines the pore +pressure came out 4 to 14 times the analytic one and grew like `1 / dh` +under refinement, while the solid velocity was off by a few percent. + +Where an element node is shared, the solid momentum (one value) and the +flux copies of the elements there are coupled through the fluid's inertia. +Per node they leave `A dq/dt + W_b q = r` over the copies. This is +integrated exactly over the step (one small matrix exponential per node, +computed once per run), so any permeability is stable. The time step is +that of the inviscid fast P wave +(`sweep.equations.biot.biot_fast_p_velocity`), which is also the speed the +C-PML is designed for. A `free_surface` side is open-pore: traction-free +with `p = 0` on its nodes. A `topography` gives a curved one. Without an +absorbing layer the boundary is traction-free and closed to flow. When the +surface's drainage layer `sqrt(2 D / omega)` is thinner than the mesh (the +Darcy regime), its first nodes are fully drained. This is about 1% on +surface records on 8 m elements, which is the FD `Biot`'s figure on 2 m +cells. The undrained limit (`b` → ∞) agrees with `ElasticVSSEM2DDeformed` +(`vp = sqrt(H / rho)`, `vs = sqrt(mu / rho)`) to the discretisation: the +volumetric `alpha^2 M` part of the modulus there passes through the +continuous `p`. + +### Anisotropic permeability and viscous attenuation + +`BiotAnisoPermSEM2DDeformed` takes the drag as a symmetric tensor +`eta kappa^-1`: the models end in `b_xx, b_xz, b_zz` instead of `b`, which +`sweep.equations.biot_params.biot_drag_tensor(eta, k_xx, k_xz, k_zz)` gives +from a permeability tensor (x right, z down -- flip the sign of `k_xz` from a +z-up convention such as SPECFEM2D's). The flux copies at a node then advance +with both components coupled, still exactly. + +`equation_options={"viscous_q": (Q0, f0)}` turns on SPECFEM2D's viscous +attenuation (`ATTENUATION_PORO_FLUID_PART`): the drag becomes a standard +linear solid, `b (1 + i w tau_e) / (1 + i w tau_s)` with +`tau_e, tau_s = (sqrt(Q0^2 + 1) ± 1) / (2 pi f0 Q0)` -- the plain drag at low +frequency, `tau_e / tau_s` times it at high frequency, `1 / Q0` at `f0`. It +costs a memory variable per flux component. `equation_options` passes any of +an equation's own `build` keywords; a name the equation does not take is +refused. + +### Fluid and dry elastic regions + +An element without a solid frame (`rho = 0`) is a fluid, one without a pore +fluid (`M = 0`) a dry elastic solid, so a sea floor, a reservoir under a +cap rock or a porous layer on a basement is one model. +`sweep.equations.biot_params.biot_fluid(k_f, rho_f)` and +`biot_elastic(vp, vs, rho)` give their coefficients; the ones a medium does +not have are zero, which is checked, and get no gradient, so an inversion +keeps the media as they are. Each element must be wholly one medium, +so give these models on the GLL nodes (`(E, ngll, ngll)` arrays, e.g. from +the centre depth of each element in `solver.mesh.z_e` -- the mesh includes +the absorbing layer) with the interfaces on element edges. In a fluid, `q` is the fluid velocity (the `qx` / `qz` +receivers) and `p` the pressure; `v` and the stresses are zero there, as `p` +and `q` are in a dry solid. + +The interfaces are those of SPECFEM2D. Between a solid (porous or dry) and a +fluid, the fluid's pressure acts on the solid and the solid's normal +velocity enters the fluid's continuity equation, by two boundary integrals +that are exact transposes of each other; with the continuous `p` this is +the open-pore condition at a porous one (pressure, traction and the fluid's +volume flux `(v + q)·n` continuous). Between a porous and a dry solid +nothing is added: the velocity and the traction are continuous through the +shared nodes and the weak form, and the pores are sealed, `q·n = 0`, which +is the natural condition of the pressure equation. (SPECFEM2D sets the whole +relative displacement to zero on those nodes.) A source loads only the fields +of the element it is in, so one that would load nothing there -- a stress or +solid-velocity source in a fluid, a pore pressure or flux in a dry solid -- +is refused, and a force in a fluid moves the fluid alone; on an interface +each field takes its load from the side that has it. An +interface may cross the side or bottom absorbing bands but not a corner +region. + +`sweep.equations.biot_params` also maps the coefficients to the parameter +families of Morency, Luo & Tromp (2009) that SPECFEM2D writes its kernels +in -- `"primary"`, `"density"` (density-normalised) and `"wavespeed"` -- and +turns a coefficient gradient into a family's kernels by the exact chain rule +(`biot_family_kernels`). + +### Validation + +`test/sem/test_biot_sem.py` checks the analytic line source (solid velocity +and pore pressure) from inviscid to the Darcy limit, with the viscous +attenuation against the frequency-dependent drag, the open-pore half-space +against wavenumber integration, the pressure-release image of a pore fluid on +its own, the pore pressure on the source's grid line, the FD records, the +undrained limit, the drag tensor's symmetries, the parameter families and +the gradient (also through the drag tensor, the viscous memory and three +media). A dry solid is `ElasticVSSEM2DDeformed` and water on a dry solid is +`AcousticElasticSEM2D`, both to round-off; across water, porous rock and +shale the Green's functions are reciprocal to round-off. +`test/sem/specfem_compare/` runs SPECFEM2D (inputs from +`make_specfem_poro_cases.py`): `run_compare_poro.py` its +`poroelastic_semi_infinite_homogeneous` example against its reference +seismograms (4e-4); `run_compare_poro_features.py` the point force, the +anisotropic permeability and the viscous drag with and without attenuation +(relative L2 0.02-0.2%; the attenuation moves the records by 0.3-1%, and +ours by the same to 3-15% of that); `run_compare_poro_coupling.py` its +`poroelastic_acoustic` example and porous-elastic and water-porous-elastic +layers (0.01-0.08% in the source band; the example 0.07-0.85% on its own +mesh, where its inviscid slow wave has about one element per wavelength, +and 0.003-0.013% with both on elements half as large). SPECFEM2D computes +no pressure in poroelastic media, so these compare velocities. + +2-D and `impl='eager'` only. Boundary saving is refused because the viscous +drag is dissipative everywhere, so the field cannot be stepped back; use +`Full()` or `Ckpt(chunks=...)`. + ## Wavefield snapshots `return_wavefield=True` (eager) returns `(seis, snapshots)` with diff --git a/src/sweep/cli.py b/src/sweep/cli.py index 200c912..4ec9adf 100644 --- a/src/sweep/cli.py +++ b/src/sweep/cli.py @@ -225,7 +225,12 @@ def _show_sem_equation(class_name, cls): f"topography: {_format_yes_no(info['topography'])}") print(f" C-PML: {_SEM_CPML[info['cpml']]}") print(f" Wavefields: {info['fields']}") - print(f" Source fields: {info['sources']}") + # an equation that names its PropSEM sources (Biot) takes exactly those + names = list(getattr(cls, "propsem_source_names", None) or ()) + if names: + print(f" Source types: {names}") + else: + print(f" Source fields: {info['sources']}") print(f" Receiver types: {info['receivers'] + info['derived']}") print(f" Needed models: {propsem_models(cls)}") print(f" Torch binding support: {_format_yes_no(binding_supported)}") diff --git a/src/sweep/equations/biot_params.py b/src/sweep/equations/biot_params.py new file mode 100644 index 0000000..67f0f7d --- /dev/null +++ b/src/sweep/equations/biot_params.py @@ -0,0 +1,162 @@ +"""Parameter families of the Biot coefficients and their kernels. + +The solvers (FD :class:`sweep.equations.Biot`, SEM ``BiotVSSEM2DDeformed``) take the eight +coefficients ``[H, mu, C, M, rho, rho_f, m, b]`` (:func:`sweep.equations.biot.biot_from_rock`), +so autograd gives the gradient with respect to those. The three families of Morency, Luo & +Tromp (2009) -- the ones SPECFEM2D writes out -- are other coordinates of the same medium: + +* ``"primary"``: ``rho, rho_f, m, b, B, C, M, mu_fr`` with ``B = H - 4 mu / 3``; +* ``"density"`` (density-normalised): ``rho, rho_f, phi_c, B / rho, C / rho_f, M / m, + mu / rho, b`` with ``phi_c = rho_f / m = phi / T``; +* ``"wavespeed"``: ``rho, rho_f, phi_c, c_PI, c_PII, c_S, r_I, b`` -- the inviscid fast-P, + slow-P and S speeds and ``r_I``, the fluid-to-solid amplitude ratio of the fast P wave. + +:func:`biot_to_family` / :func:`biot_from_family` map between them exactly, and +:func:`biot_family_kernels` turns a coefficient gradient into a family's kernels by the +exact chain rule (a vector-Jacobian product through :func:`biot_from_family`). SPECFEM2D's +printed wave-speed kernels are not that chain rule. + +With an anisotropic permeability the drag is a tensor ``[[b_xx, b_xz], [b_xz, b_zz]]`` +(:func:`biot_drag_tensor`); the coefficient lists then end in its three components, which +every family carries unchanged in place of ``b``. + +:func:`biot_fluid` and :func:`biot_elastic` give the coefficients of a fluid and of a dry +elastic solid, the other two media of the SEM ``BiotVSSEM2DDeformed`` (a fluid has no solid +frame, ``rho = 0``; a dry solid no pore fluid, ``M = 0``). The families are those of a +porous medium only. +""" + +from __future__ import annotations + +import torch + +__all__ = ["BIOT_FAMILIES", "biot_drag_tensor", "biot_fluid", "biot_elastic", + "biot_to_family", "biot_from_family", "biot_family_kernels"] + +BIOT_FAMILIES = { + "primary": ["rho", "rho_f", "m", "b", "B", "C", "M", "mu_fr"], + "density": ["rho", "rho_f", "phi_c", "B_bar", "C_bar", "M_bar", "mu_bar", "b"], + "wavespeed": ["rho", "rho_f", "phi_c", "c_PI", "c_PII", "c_S", "r_I", "b"], +} + + +def biot_drag_tensor(eta, k_xx, k_xz, k_zz): + """The drag ``eta kappa^-1`` of a permeability tensor ``[[k_xx, k_xz], [k_xz, k_zz]]`` + (m^2, x right and z down -- flip the sign of ``k_xz`` from a z-up convention such as + SPECFEM2D's), as ``(b_xx, b_xz, b_zz)``. Tensors or floats; gradients flow.""" + det = k_xx * k_zz - k_xz * k_xz + return eta * k_zz / det, -eta * k_xz / det, eta * k_xx / det + + +def biot_fluid(k_f, rho_f, drag_tensor=False): + """The coefficients ``[H, mu, C, M, rho, rho_f, m, b]`` of a fluid of bulk modulus + ``k_f`` and density ``rho_f``: no solid frame (``rho = 0``), the fluid's own + compliance and inertia in ``M`` and ``m``, nothing else. ``drag_tensor=True`` ends the + list in ``(b_xx, b_xz, b_zz)`` instead of ``b``. Tensors or floats.""" + k_f, rho_f = torch.broadcast_tensors(_t(k_f), _t(rho_f)) + z = torch.zeros_like(k_f) + return [z, z, z, k_f, z, z, rho_f] + [z] * (3 if drag_tensor else 1) + + +def biot_elastic(vp, vs, rho, drag_tensor=False): + """The coefficients ``[H, mu, C, M, rho, rho_f, m, b]`` of a dry elastic solid of P and + S velocities ``vp``, ``vs`` and density ``rho``: ``H = rho vp^2``, ``mu = rho vs^2``, + no pore fluid (``M = 0``) and nothing else. ``drag_tensor=True`` ends the list in + ``(b_xx, b_xz, b_zz)`` instead of ``b``. Tensors or floats.""" + vp, vs, rho = torch.broadcast_tensors(_t(vp), _t(vs), _t(rho)) + z = torch.zeros_like(vp) + return [rho * vp * vp, rho * vs * vs, z, z, rho, z, z] + [z] * (3 if drag_tensor else 1) + + +def _t(x): + """A tensor; a Python number becomes a float64 one.""" + return x if isinstance(x, torch.Tensor) else torch.as_tensor(x, dtype=torch.float64) + + +def _split(coef): + coef = list(coef) + if len(coef) not in (8, 10): + raise ValueError(f"Biot coefficients are [H, mu, C, M, rho, rho_f, m, b] or end in " + f"(b_xx, b_xz, b_zz); got {len(coef)} entries") + return coef[:7], coef[7:] + + +def _fast_p_ratio(h, mu, c, bm, rho, rho_f, m): + g = (rho * c - rho_f * h) / m + s0 = h - rho * bm / m + return (s0 + torch.sqrt(s0 * s0 + 4 * (c - rho_f * bm / m) * g)) / (2 * g) + + +def biot_to_family(family, coef): + """``coef`` (``[H, mu, C, M, rho, rho_f, m, b]``, or the drag tensor in place of ``b``) + in the coordinates of ``family`` (see :data:`BIOT_FAMILIES`), as a list of tensors.""" + (h, mu, c, bm, rho, rho_f, m), drag = _split(_t(x) for x in coef) + big_b = h - 4.0 / 3.0 * mu + if family == "primary": + return [rho, rho_f, m, *drag, big_b, c, bm, mu] + if family == "density": + return [rho, rho_f, rho_f / m, big_b / rho, c / rho_f, bm / m, mu / rho, *drag] + if family == "wavespeed": + a = rho * m - rho_f ** 2 + bq = h * m + bm * rho - 2 * c * rho_f + disc = torch.sqrt(bq * bq - 4 * a * (h * bm - c * c)) + return [rho, rho_f, rho_f / m, torch.sqrt((bq + disc) / (2 * a)), + torch.sqrt((bq - disc) / (2 * a)), torch.sqrt(mu / (rho - rho_f ** 2 / m)), + _fast_p_ratio(h, mu, c, bm, rho, rho_f, m), *drag] + raise ValueError(f"unknown family {family!r}; choose from {sorted(BIOT_FAMILIES)}") + + +def biot_from_family(family, params): + """The inverse of :func:`biot_to_family`: ``[H, mu, C, M, rho, rho_f, m, b...]``.""" + params = [_t(x) for x in params] + n = len(params) + if n not in (8, 10): + raise ValueError(f"a {family!r} parameter list has 8 entries (10 with a drag " + f"tensor); got {n}") + if family == "primary": + rho, rho_f, m = params[:3] + drag = params[3:n - 4] + big_b, c, bm, mufr = params[n - 4:] + return [big_b + 4.0 / 3.0 * mufr, mufr, c, bm, rho, rho_f, m, *drag] + drag = params[7:] + if family == "density": + rho, rho_f, phic, bbar, cbar, mbar, mubar = params[:7] + m = rho_f / phic + return [rho * (bbar + 4.0 / 3.0 * mubar), rho * mubar, rho_f * cbar, m * mbar, + rho, rho_f, m, *drag] + if family == "wavespeed": + rho, rho_f, phic, cp1, cp2, cs, r = params[:7] + m = rho_f / phic + a1, a2 = cp1 * cp1, cp2 * cp2 + u1, u2 = rho * r + rho_f, rho_f * r + m + nn = rho * r * r + 2 * rho_f * r + m + return [a2 * rho + (a1 - a2) * u1 * u1 / nn, cs * cs * (rho - rho_f * rho_f / m), + a2 * rho_f + (a1 - a2) * u1 * u2 / nn, a2 * m + (a1 - a2) * u2 * u2 / nn, + rho, rho_f, m, *drag] + raise ValueError(f"unknown family {family!r}; choose from {sorted(BIOT_FAMILIES)}") + + +def biot_family_kernels(family, coef, grads, relative=True): + """``{name: kernel}`` of ``family`` from the gradient ``grads`` of a misfit with respect + to ``coef`` (same layout: one tensor per coefficient, e.g. what ``loss.backward()`` + leaves on the models given to the propagator). + + ``relative=True`` gives ``y dJ/dy`` (Morency et al.'s ``K_y``, the misfit change per + relative change of ``y``), ``False`` the plain ``dJ/dy``. Computed in float64 by the + exact chain rule through :func:`biot_from_family`.""" + cf = [_t(c).detach().double() for c in coef] + g = [_t(x).detach().double() for x in grads] + if len(g) != len(cf): + raise ValueError(f"{len(g)} gradients for {len(cf)} coefficients") + y = [t.detach().clone().requires_grad_() for t in biot_to_family(family, cf)] + back = biot_from_family(family, y) + jac = torch.autograd.grad(back, y, grad_outputs=g, allow_unused=True) + names = list(BIOT_FAMILIES[family]) + if len(cf) == 10: + i = names.index("b") + names[i:i + 1] = ["b_xx", "b_xz", "b_zz"] + out = {} + for name, yy, jj in zip(names, y, jac): + jj = torch.zeros_like(yy) if jj is None else jj + out[name] = (yy.detach() * jj) if relative else jj + return out diff --git a/src/sweep/sem/_assembly.py b/src/sweep/sem/_assembly.py index a24dd9e..7a2ec43 100644 --- a/src/sweep/sem/_assembly.py +++ b/src/sweep/sem/_assembly.py @@ -38,6 +38,18 @@ def assemble(r_e: torch.Tensor, ibool: torch.Tensor, like: torch.Tensor): return out.index_add(0, idx, r_e.reshape(-1, k)).unbind(-1) +def element_mean(mesh, f_e): + """An element-local field ``(E, ngll, ...)`` as a global one ``(n_global,)``: + at a node the elements share, the ``w |J|``-weighted mean of their values + (the lumped L2 projection onto the continuous space). For snapshots of the + element-local fields; a continuous field comes back unchanged.""" + idx = mesh.ibool.reshape(-1) + w = mesh.w_jac + num = f_e.new_zeros(mesh.n_global).index_add(0, idx, (w * f_e).reshape(-1)) + den = w.new_zeros(mesh.n_global).index_add(0, idx, w.reshape(-1)) + return num / den + + def assemble_each(fields, ibool, like): """Assemble several element-local fields to global nodes, one single-field ``index_add`` each — the unstacked counterpart of ``assemble(stack(fields))``. diff --git a/src/sweep/sem/equations/biot_vs_deformed.py b/src/sweep/sem/equations/biot_vs_deformed.py new file mode 100644 index 0000000..0eacb60 --- /dev/null +++ b/src/sweep/sem/equations/biot_vs_deformed.py @@ -0,0 +1,645 @@ +"""1st-order velocity-stress Biot poroelastic SEM on a Quad9 deformed mesh. + +The low-frequency Biot equations in the form of the FD :class:`sweep.equations.Biot` +(solid velocity ``v``, Darcy flux ``q = phi (v_f - v)``, total stress ``sigma``, pore +pressure ``p``; Biot 1962, Carcione *Wave Fields in Real Media* ch. 7):: + + rho dv/dt + rho_f dq/dt = div sigma + rho_f dv/dt + m dq/dt + B q = -grad p + d sigma/dt = (H - 2 mu) div v I + mu (grad v + grad v^T) + C div q I + d p/dt = -C div v - M div q + +with the coefficients ``[H, mu, C, M, rho, rho_f, m, b]`` of +:func:`sweep.equations.biot.biot_from_rock`; ``B = b I`` here, a tensor ``eta kappa^-1`` in +:class:`BiotAnisoPermSEM2DDeformed`. + +Discretisation +-------------- +The solid is :class:`ElasticVSSEM2DDeformed`'s: ``v_x, v_z`` continuous (CG) nodal +fields, the total stress element-local, the momentum equation in weak form with a lumped +mass. The pore fluid is the pressure form of the acoustic SEM: ``p`` is a CG nodal field +whose equation is taken in weak form, and the Darcy flux ``q_x, q_z`` is element-local, +driven by the strong gradient of ``p``:: + + M_p dp/dt = -A(w alpha div v) + D^T W q, M_p = A(w / M), alpha = C / M, + sigma = sigma' - alpha p I, d sigma'/dt = lambda_d div v I + mu (grad v + grad v^T) + +(``A`` assembles, ``D^T W`` is the weak divergence), so that a pore fluid on its own +(``rho_f = C = 0``) is exactly the acoustic pressure SEM. The other placement -- ``q`` +continuous and ``p`` element-local, as in SPECFEM2D's poroelastic solver -- leaves the +relative flow, which has no shear stiffness, on a continuous field: a point source then +excites grid-scale modes that run along the grid lines through it, and the pore pressure +there grows like ``1 / dh`` under refinement. The two couplings (``alpha p`` in the +momentum, ``alpha div v`` in the pressure) are transposes of each other, as are ``D^T W q`` +and the strong ``grad p``, so the semi-discrete system conserves the Biot energy up to the +viscous drag. + +The two momentum equations meet at every node: the solid's (one, CG) and the element +copies of the flux there (one per element sharing the node). Eliminating ``dv/dt`` leaves, +per node, ``A dq/dt + W_B q = r`` over the copies (and both components when ``B`` is a +tensor), with ``A = diag(w m) - w rho_f (w rho_f)^T / M_rho``; with the forces frozen over +the step this is integrated exactly, ``q <- E q + Psi r`` (``E`` and ``Psi`` from one +matrix exponential per node, computed in :meth:`init_state` for the run's ``dt``), and +then ``v <- v + (dt F - sum w rho_f (q_new - q)) / M_rho``. So any permeability is stable, +the time step is that of the inviscid fast P wave (see +:func:`sweep.equations.biot.biot_fast_p_velocity`), and ``b -> inf`` gives the undrained +elastic medium. + +Viscous attenuation (``viscous_q=(Q0, f0)``, SPECFEM2D's ``ATTENUATION_PORO_FLUID_PART``): +the drag becomes a standard linear solid, ``B (1 + i w tau_e) / (1 + i w tau_s)`` with +``tau_e, tau_s = (sqrt(Q0^2 + 1) +- 1) / (2 pi f0 Q0)`` -- ``B`` at low frequency, +``B tau_e / tau_s`` at high frequency, ``1 / Q0`` at ``f0``. The unrelaxed drag goes into +the exact step and a memory ``r``, ``dr/dt = -(r + (1 - tau_e / tau_s) B q) / tau_s``, acts +as a force on the flux. + +Fluid and elastic regions +------------------------- +An element without a solid frame (``rho = 0``) is a fluid, one without a pore fluid +(``M = 0``) a dry elastic solid -- see :func:`sweep.equations.biot_params.biot_fluid` and +:func:`~sweep.equations.biot_params.biot_elastic`; every other coefficient that does not +exist there must be zero, and an element is wholly one medium. In a fluid ``q`` is the +fluid velocity and ``m`` the fluid density (the equations above with ``rho = rho_f = H = +mu = C = b = 0``), so the fluid is the acoustic pressure SEM; a dry solid is +:class:`ElasticVSSEM2DDeformed`. ``v`` lives on the nodes of the solid elements, ``p`` on +those of the fluid-filled ones (masked inverse masses keep the others at zero) and the +flux copies of a dry element stay zero. The interfaces: + +* solid (porous or dry) -- fluid: ``-∮ p n φ`` on the solid momentum and ``+∮ ψ v·n`` on + the pressure (``n`` out of the solid), the two exact transposes of + :class:`sweep.sem.coupling.FluidSolidInterface`; with the continuous ``p`` that is the + open-pore condition at a porous one (pressure, traction and the fluid volume flux + ``(v + q)·n`` continuous), SPECFEM2D's; +* porous -- dry solid: nothing to add. ``v`` and the traction are continuous through the + shared nodes and the weak form, and the dropped boundary integral of ``D^T W q`` is the + sealed pores ``q·n = 0`` (SPECFEM2D sets the whole ``w = 0`` on those nodes). + +A source loads only the fields of the element it sits in -- a stress source in a fluid or +a pore pressure in a dry solid is refused, a force in a fluid moves the fluid alone -- and +on an interface it is shared among the elements that carry each field. + +Boundaries: an open-pore free surface (traction-free, ``p = 0``) on the ``dirichlet`` +sides -- ``p`` pinned to zero on their nodes, the traction natural; elsewhere the boundary +integrals are dropped (traction-free, no flow), behind the C-PML. C-PML: the profile of +:class:`ElasticVSSEM2DDeformed`, field memories for the weak forms (stress, flux) and +gradient memories for the strong ones (``grad v``, ``grad p``): 12 element fields. An +interface may cross the side or bottom bands, with the plain coupling there, but no corner +region (as in :class:`~sweep.sem.equations.acoustic_elastic2d.AcousticElasticSEM2D`). + +Sources +------- +``src`` is aligned with ``(v_x, v_z, p, q_x, q_z)``. ``src[2]`` injects pore pressure: +``p`` jumps by ``dt src / A(w)`` at the start of the step and the total stress does not move +(FD ``Biot``'s ``p``). The others follow ``source_mode``: + +* ``"load"`` (default): weak-form loads, ``src[0:2]`` on the total momentum and ``src[3:5]`` + on the relative flow (shared among the flux copies at a node by ``w``). A moment tensor + ``-M : grad(phi)`` on the first injects total stress (FD's ``sxx`` / ``szz`` / ``sxz``); the + same force on both is SPECFEM2D's poroelastic point force. +* ``"velocity"``: FD ``Biot``'s ``vx`` / ``vz`` / ``qx`` / ``qz``, after the momentum update: + ``v`` jumps by ``dt src / M_rho`` with the flux unchanged, or every flux copy by + ``dt src / A(w m)`` with ``v`` unchanged. + +Eager only, no boundary saving: the viscous drag is dissipative at every node, so the +reverse reconstruction does not exist (as on the FD ``Biot``). +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from types import SimpleNamespace +from typing import ClassVar, Optional + +import torch + +from .._assembly import assemble, phys_grads +from ..coupling import FluidSolidInterface, assert_no_corner_pml_on_interface +from ..cpml.elastic_vs import SimpleElasticCPMLProfile, build_simple_elastic_cpml +from ..fields import FieldSpec, check_element_arrays, init_state_from_specs +from ..mesh.boundaries import free_surface_node_idx +from ..mesh.domains import DomainMask, masked_inv_mass +from ..mesh.quad9 import Quad9Mesh + +_SOURCE_MODES = ("load", "velocity") + + +def _node_slots(ibool): + """``(node, slot, kmax)`` of every element-local node in flattened ``ibool`` order: its + global node and its rank among the element copies of that node.""" + g = ibool.reshape(-1) + order = torch.argsort(g, stable=True) + gs = g[order] + pos = torch.arange(gs.numel(), device=g.device) + start = torch.ones_like(gs, dtype=torch.bool) + start[1:] = gs[1:] != gs[:-1] + run = torch.cummax(torch.where(start, pos, torch.zeros_like(pos)), 0).values + slot = torch.empty_like(pos) + slot[order] = pos - run + return g, slot, int(slot.max()) + 1 + + +def _media(H_e, mu_e, C_e, M_e, rho_e, rho_f_e, m_e, b, tensor): + """``(solid, fluid)``, ``(E,)`` bool: the elements with a solid frame (``rho > 0``) and + those with a pore fluid (``M > 0``); both is a porous element, one alone a dry solid or + a fluid. The coefficients a medium does not have must be zero.""" + drag = ("b_xx", "b_xz", "b_zz") if tensor else ("b", "b_xz", "b") + names = ("H", "mu", "C", "M", "rho", "rho_f", "m") + drag + vals = dict(zip(names, (H_e, mu_e, C_e, M_e, rho_e, rho_f_e, m_e, *b))) + flags = [] + for key, what in (("rho", "a solid frame (rho > 0)"), ("M", "a pore fluid (M > 0)")): + x = vals[key].flatten(1) + if bool((x < 0).any()): + raise ValueError(f"{key} must be >= 0") + has, some = (x > 0).all(1), (x > 0).any(1) + if bool((some & ~has).any()): + bad = int((some & ~has).nonzero()[0]) + raise ValueError(f"element {bad} has {what} on some of its nodes only; an " + f"element is wholly porous, wholly fluid or wholly dry") + flags.append(has) + solid, fluid = flags + if bool((~solid & ~fluid).any()): + bad = int((~solid & ~fluid).nonzero()[0]) + raise ValueError(f"element {bad} has neither a solid frame (rho > 0) nor a pore " + f"fluid (M > 0)") + for where, zero, what in ((solid & ~fluid, ("C", "rho_f", "m") + drag, + "a dry elastic element (M = 0)"), + (fluid & ~solid, ("H", "mu", "C", "rho_f") + drag, + "a fluid element (rho = 0)")): + if bool(where.any()): + bad = [k for k in dict.fromkeys(zero) if bool((vals[k][where] != 0).any())] + if bad: + raise ValueError(f"{what} has nonzero {bad}: they do not exist there " + f"and must be 0") + if bool((m_e[fluid] <= 0).any()): + raise ValueError("m must be > 0 where there is a fluid (M > 0)") + return solid, fluid + + +def viscous_relaxation_times(q0: float, f0: float): + """``(tau_e, tau_s)`` of SPECFEM2D's viscous standard linear solid, quality ``q0`` at + ``f0`` Hz.""" + if not (q0 > 0 and f0 > 0): + raise ValueError(f"viscous_q needs Q0 > 0 and f0 > 0, got ({q0}, {f0})") + s = math.sqrt(q0 * q0 + 1.0) + return (s + 1.0) / (2.0 * math.pi * f0 * q0), (s - 1.0) / (2.0 * math.pi * f0 * q0) + + +@dataclass +class BiotVSSEM2DDeformed: + """1st-order velocity-stress Biot poroelastic SEM on a deformed Quad9 mesh.""" + + # ---- Field specification: the positional state layout (see fields.py) ---- # + field_specs: ClassVar[list] = [ + FieldSpec("v_x", supports_source=True, supports_receiver=True), + FieldSpec("v_z", supports_source=True, supports_receiver=True), + FieldSpec("p", supports_source=True, supports_receiver=True), + FieldSpec("q_x", supports_source=True, supports_receiver=True, layout="element"), + FieldSpec("q_z", supports_source=True, supports_receiver=True, layout="element"), + FieldSpec("sigma_xx", supports_receiver=True, layout="element"), + FieldSpec("sigma_zz", supports_receiver=True, layout="element"), + FieldSpec("sigma_xz", supports_receiver=True, layout="element"), + FieldSpec("r_x", internal=True, layout="element"), # viscous memory (viscous_q) + FieldSpec("r_z", internal=True, layout="element"), + FieldSpec("m_sxx_x", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_sxz_z", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_sxz_x", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_szz_z", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_qx_x", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_qz_z", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_px", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_pz", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_vxx", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_vzz", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_vxz", internal=True, boundary_related=True, layout="element"), + FieldSpec("m_vzx", internal=True, boundary_related=True, layout="element"), + ] + + # FD Biot's receiver names + derived_receiver_types: ClassVar[dict] = { + "vx": ("v_x", None), "vz": ("v_z", None), "qx": ("q_x", None), "qz": ("q_z", None), + "sxx": ("sigma_xx", None), "szz": ("sigma_zz", None), "sxz": ("sigma_xz", None)} + + # PropSEM's source names: kind and payload over (v_x, v_z, p, q_x, q_z), or the + # (Mxx, Mzz, Mxz) of a moment. "velocity": FD Biot's injections (source_mode + # "velocity"); "force": weak-form loads (source_mode "load"), fx / fz being SPECFEM2D's + # poroelastic point force; "injection": the pore pressure, either mode; "moment": the + # total-stress moment tensor, driven by the wavelet's integral. + propsem_source_names: ClassVar[dict] = { + "vx": ("velocity", (1.0, 0.0, 0.0, 0.0, 0.0)), + "vz": ("velocity", (0.0, 1.0, 0.0, 0.0, 0.0)), + "qx": ("velocity", (0.0, 0.0, 0.0, 1.0, 0.0)), + "qz": ("velocity", (0.0, 0.0, 0.0, 0.0, 1.0)), + "v_x": ("velocity", (1.0, 0.0, 0.0, 0.0, 0.0)), + "v_z": ("velocity", (0.0, 1.0, 0.0, 0.0, 0.0)), + "q_x": ("velocity", (0.0, 0.0, 0.0, 1.0, 0.0)), + "q_z": ("velocity", (0.0, 0.0, 0.0, 0.0, 1.0)), + "fx": ("force", (1.0, 0.0, 0.0, 1.0, 0.0)), + "fz": ("force", (0.0, 1.0, 0.0, 0.0, 1.0)), + "p": ("injection", (0.0, 0.0, 1.0, 0.0, 0.0)), + "sxx": ("moment", (1.0, 0.0, 0.0)), "szz": ("moment", (0.0, 1.0, 0.0)), + "sxz": ("moment", (0.0, 0.0, 1.0)), + } + default_source_type: ClassVar[tuple] = ("sxx", "szz") + default_receiver_type: ClassVar[tuple] = ("v_x", "v_z") + tensor_drag: ClassVar[bool] = False + + # PropSEM refuses memory=BoundarySaving() with this reason + bs_refusal: ClassVar[str] = ( + "Biot has no boundary saving: the viscous drag is dissipative at every node, so " + "stepping back multiplies round-off by exp(gamma dt) per step and the reverse " + "reconstruction does not exist. Use memory=Full() or Ckpt(chunks=...).") + + mesh: Quad9Mesh + H_e: torch.Tensor # (E, ngll, ngll) undrained P-wave modulus + mu_e: torch.Tensor + C_e: torch.Tensor + M_e: torch.Tensor + rho_e: torch.Tensor + rho_f_e: torch.Tensor + m_e: torch.Tensor + b_e: tuple # (b_xx, b_xz, b_zz) element-local, relaxed + dirichlet: tuple # the open-pore sides: p = 0 on their nodes + dirichlet_idx: torch.Tensor + viscous_q: Optional[tuple] # (Q0, f0) of the viscous attenuation, or None + + mass_rho: torch.Tensor # (n_global,) A(w rho) + _inv_M: torch.Tensor # 1 / mass_rho, 0 off the solid + _inv_Mm: torch.Tensor # 1 / A(w m): a flux injected at a point + _inv_Mp: torch.Tensor # 1 / A(w / M), the pressure's lumped compliance + _inv_M1: torch.Tensor # 1 / A(w) over the fluid: a pressure injected + _lam_d: torch.Tensor # drained lambda, H - 2 mu - C^2 / M + _alpha: torch.Tensor # C / M (0 in a dry solid) + # the media (see "Fluid and elastic regions"): (E,) bool, the flux mask (E, 1, 1), + # the solid-fluid faces + _solid_e: torch.Tensor + _fluid_e: torch.Tensor + _qmask: torch.Tensor + ifc: FluidSolidInterface + _tau: tuple # (tau_e, tau_s), equal without viscous_q + # the element copies of the flux at every node: flattened ibool -> node * kmax + slot + _slot_idx: torch.Tensor # (E ngll^2,) + _kmax: int + _wm: torch.Tensor # (n_global, kmax) w m per copy (1 on empty slots) + _wr: torch.Tensor # w rho_f (0 on empty slots) + _wb: tuple # w (b_xx, b_xz, b_zz) tau_e / tau_s: unrelaxed + _share: torch.Tensor # w / A(w): a copy's share of a nodal flux load + + cpml: Optional[SimpleElasticCPMLProfile] + + _D: torch.Tensor # (ngll, ngll) + _w_J: torch.Tensor + _ixx: torch.Tensor + _ixz: torch.Tensor + _iex: torch.Tensor + _iez: torch.Tensor + + source_mode: str = "load" # see the module docstring + _drag: tuple = None # (dt, E, Psi), set by init_state + + @classmethod + def build( + cls, + mesh: Quad9Mesh, + H_e: torch.Tensor, + mu_e: torch.Tensor, + C_e: torch.Tensor, + M_e: torch.Tensor, + rho_e: torch.Tensor, + rho_f_e: torch.Tensor, + m_e: torch.Tensor, + b_e: torch.Tensor, + *, + dirichlet: tuple = (), + cpml: Optional[SimpleElasticCPMLProfile] = None, + viscous_q: Optional[tuple] = None, + ) -> "BiotVSSEM2DDeformed": + """Build the operators from the eight Biot coefficients at the GLL nodes; + ``dirichlet`` names the open-pore sides, ``viscous_q = (Q0, f0)`` turns on the + viscous attenuation.""" + check_element_arrays(mesh, b_e=b_e) + return cls._build(mesh, H_e, mu_e, C_e, M_e, rho_e, rho_f_e, m_e, + (b_e, torch.zeros_like(b_e), b_e), dirichlet, cpml, viscous_q) + + @classmethod + def _build(cls, mesh, H_e, mu_e, C_e, M_e, rho_e, rho_f_e, m_e, b, dirichlet, cpml, + viscous_q): + check_element_arrays(mesh, H_e=H_e, mu_e=mu_e, C_e=C_e, M_e=M_e, rho_e=rho_e, + rho_f_e=rho_f_e, m_e=m_e) + b_xx, b_xz, b_zz = b + if bool(((b_xx < 0) | (b_zz < 0) | (b_xx * b_zz < b_xz * b_xz)).any()): + raise ValueError("the drag eta / kappa must be positive semi-definite") + solid, fluid = _media(H_e, mu_e, C_e, M_e, rho_e, rho_f_e, m_e, b, cls.tensor_drag) + tau = (1.0, 1.0) if viscous_q is None else viscous_relaxation_times(*viscous_q) + device = H_e.device + w = mesh.w_jac + ibool = mesh.ibool + like = H_e.new_zeros(mesh.n_global) + + def lump(c): + return assemble(w * c, ibool, like) + + def alive(elements): + out = torch.zeros(mesh.n_global, dtype=torch.bool, device=device) + return out.index_fill(0, ibool[elements].reshape(-1), True) + + node, slot, kmax = _node_slots(ibool) + slot_idx = node * kmax + slot + + def slots(x, fill): + out = torch.full((mesh.n_global * kmax,), fill, dtype=x.dtype, device=device) + return out.index_copy(0, slot_idx, x.reshape(-1)).view(mesh.n_global, kmax) + + v_alive, p_alive = alive(solid), alive(fluid) + fl = fluid.view(-1, 1, 1) + qmask = fl.to(H_e.dtype) + # the coefficients a medium does not have are zero (checked above); through the + # masks their gradient is zero too, so an update keeps them so -- a fluid's rho + # would otherwise feel the solid's mass on the nodes it shares with a solid + sm = solid.view(-1, 1, 1).to(H_e.dtype) + pm = sm * qmask + H_e, mu_e, rho_e, C_e, rho_f_e = H_e * sm, mu_e * sm, rho_e * sm, C_e * pm, rho_f_e * pm + m_e = m_e * qmask + b_xx, b_xz, b_zz = b_xx * pm, b_xz * pm, b_zz * pm + inv_M_e = torch.where(fl, 1.0 / torch.where(fl, M_e, torch.ones_like(M_e)), 0.0) + mass_rho = lump(rho_e) + mass_1 = lump(qmask.expand_as(M_e)) + inv_M1 = masked_inv_mass(mass_1, p_alive) + alpha = C_e * inv_M_e + unrelax = tau[0] / tau[1] + ifc = FluidSolidInterface.detect(mesh, DomainMask.from_element_flags(mesh, ~solid)) + if cpml is not None and ifc.n_face: + node_profile = {k: like.index_copy(0, ibool.reshape(-1), getattr(cpml, k).reshape(-1)) + for k in ("b_x", "b_z")} + assert_no_corner_pml_on_interface(ifc, SimpleNamespace(**node_profile)) + return cls( + mesh=mesh, H_e=H_e, mu_e=mu_e, C_e=C_e, M_e=M_e, rho_e=rho_e, + rho_f_e=rho_f_e, m_e=m_e, b_e=(b_xx, b_xz, b_zz), + dirichlet=tuple(dirichlet), + dirichlet_idx=free_surface_node_idx(mesh, dirichlet, device), + viscous_q=None if viscous_q is None else tuple(float(v) for v in viscous_q), + mass_rho=mass_rho, _inv_M=masked_inv_mass(mass_rho, v_alive), + _inv_Mm=masked_inv_mass(lump(m_e), p_alive), + _inv_Mp=masked_inv_mass(lump(inv_M_e), p_alive), _inv_M1=inv_M1, + _lam_d=H_e - 2.0 * mu_e - alpha * C_e, _alpha=alpha, + _solid_e=solid, _fluid_e=fluid, _qmask=qmask, ifc=ifc, _tau=tau, + _slot_idx=slot_idx, _kmax=kmax, + _wm=slots(torch.where(fl, w * m_e, 1.0), 1.0), _wr=slots(w * rho_f_e, 0.0), + _wb=tuple(slots(w * bb * unrelax, 0.0) for bb in (b_xx, b_xz, b_zz)), + _share=slots(w * qmask, 0.0) * inv_M1[:, None], + cpml=cpml, + _D=mesh.deriv, _w_J=w, + _ixx=mesh.inv_J[..., 0, 0], _ixz=mesh.inv_J[..., 0, 1], + _iex=mesh.inv_J[..., 1, 0], _iez=mesh.inv_J[..., 1, 1], + ) + + # The C-PML builder ``cpml=`` expects, and the speed PropSEM designs it for. + cpml_builder: ClassVar[object] = staticmethod(build_simple_elastic_cpml) + + @staticmethod + def pml_speed(models) -> float: + """The largest P velocity of the coefficients (the drag does not enter): the + inviscid fast P wave in a porous medium, ``sqrt(H / rho)`` in a dry solid, + ``sqrt(M / m)`` in a fluid.""" + from ...equations.biot import biot_wave_velocities + h, mu, c, bm, rho, rho_f, m = torch.broadcast_tensors( + *(torch.as_tensor(x).detach().double() for x in models[:7])) + solid, fluid = rho > 0, bm > 0 + speeds = [] + if bool((solid & fluid).any()): + sel = solid & fluid + speeds.append(biot_wave_velocities([x[sel] for x in (h, mu, c, bm, rho, rho_f, m)])[0]) + if bool((solid & ~fluid).any()): + speeds.append(torch.sqrt(h[solid & ~fluid] / rho[solid & ~fluid])) + if bool((fluid & ~solid).any()): + speeds.append(torch.sqrt(bm[fluid & ~solid] / m[fluid & ~solid])) + return float(torch.cat([v.reshape(-1) for v in speeds]).max()) + + @property + def source_field_elements(self): + """The elements each source field ``(v_x, v_z, p, q_x, q_z)`` lives in: a source + loads a field only through those (see :func:`sweep.sem.propagator._stamp_loads`).""" + s, f = self._solid_e, self._fluid_e + return (s, s, f, f, f) + + # ------------------------------------------------------------------ # + # The drag: one matrix exponential per node, for the run's dt + # ------------------------------------------------------------------ # + def drag_operators(self, dt: float, chunk: int = 1 << 16): + """``(E, Psi)`` advancing the flux copies at a node by ``q <- E q + Psi r`` over a + step of frozen forces ``r``: from ``A dq/dt + W_B q = r``, + ``exp([[-A^-1 W_B, A^-1], [0, 0]] dt) = [[E, Psi], [0, I]]``, in float64. Each is + ``(n_global, kmax, kmax)`` (both components alike), or ``(n_global, 2 kmax, + 2 kmax)`` over ``(q_x copies, q_z copies)`` with a drag tensor. ``chunk`` nodes at a + time: the exponential's work arrays are ~10 times its input.""" + wm, wr = self._wm.double(), self._wr.double() + wb = [t.double() for t in self._wb] + inv_rho = masked_inv_mass(self.mass_rho.double(), self._inv_M != 0) + dtype = self.mass_rho.dtype + Es, Psis = [], [] + for i in range(0, wm.shape[0], chunk): + n = slice(i, i + chunk) + A = torch.diag_embed(wm[n]) - wr[n, :, None] * wr[n, None, :] * inv_rho[n, None, None] + if self.tensor_drag: + z = torch.zeros_like(A) + A = torch.cat([torch.cat([A, z], 2), torch.cat([z, A], 2)], 1) + wxx, wxz, wzz = (torch.diag_embed(t[n]) for t in wb) + WB = torch.cat([torch.cat([wxx, wxz], 2), torch.cat([wxz, wzz], 2)], 1) + else: + WB = torch.diag_embed(wb[0][n]) + k = A.shape[-1] + A_inv = torch.linalg.inv(A) + top = torch.cat([-A_inv @ WB, A_inv], 2) + X = torch.linalg.matrix_exp(torch.cat([top, torch.zeros_like(top)], 1) * dt) + Es.append(X[:, :k, :k].to(dtype)) + Psis.append(X[:, :k, k:].to(dtype)) + return torch.cat(Es), torch.cat(Psis) + + def init_state(self, dt: float, device, dtype) -> tuple: + """Zero fields, and the drag operators for this ``dt`` -- computed here, once per + run, outside any checkpointed chunk of the time loop.""" + if self._drag is None or self._drag[0] != dt: + self._drag = (dt, *self.drag_operators(dt)) + return init_state_from_specs(self, device, dtype) + + def _to_slots(self, x): + """Element-local ``(E, ngll, ngll)`` -> ``(n_global, kmax)``, zero on empty slots.""" + k = self._kmax + out = x.new_zeros(self.mesh.n_global * k).index_add(0, self._slot_idx, x.reshape(-1)) + return out.view(self.mesh.n_global, k) + + def _from_slots(self, s, like): + return s.reshape(-1)[self._slot_idx].view_as(like) + + def _weak_div(self, *pairs): + """``sum_g w (d(phi)/dx X + d(phi)/dz Z)`` at the global nodes for every element-local + pair ``(X, Z)`` -- :class:`ElasticVSSEM2DDeformed`'s stress divergence for the + momentum, ``D^T W q`` for the pressure.""" + w_J, ixx, ixz, iex, iez = self._w_J, self._ixx, self._ixz, self._iex, self._iez + F_beta = torch.stack([w_J * (X * ixx + Z * ixz) for X, Z in pairs], dim=-1) + F_alpha = torch.stack([w_J * (X * iex + Z * iez) for X, Z in pairs], dim=-1) + R_e = (torch.einsum("db,eadk->eabk", self._D, F_beta) + + torch.einsum("ga,egbk->eabk", self._D, F_alpha)) + like = pairs[0][0].new_zeros(self.mesh.n_global) + if len(pairs) == 1: + return assemble(R_e[..., 0], self.mesh.ibool, like) + return assemble(R_e, self.mesh.ibool, like) + + def _pin(self, p): + """``p = 0`` on the open-pore surface.""" + idx = self.dirichlet_idx + return p if idx.numel() == 0 else p.index_fill(0, idx, 0.0) + + def _advance_flux(self, q_x, q_z, r_x, r_z): + """The flux copies after the step, ``E q + Psi r``, slot layout ``(n_global, k)``.""" + _, E, Psi = self._drag + if self.tensor_drag: + k = self._kmax + qs = torch.cat([q_x, q_z], 1) + new = torch.einsum("nij,nj->ni", E, qs) + torch.einsum("nij,nj->ni", Psi, + torch.cat([r_x, r_z], 1)) + return new[:, :k], new[:, k:] + return tuple(torch.einsum("nij,nj->ni", E, q) + torch.einsum("nij,nj->ni", Psi, r) + for q, r in ((q_x, r_x), (q_z, r_z))) + + # ================================================================== # + # SEMPropagator interface: init_state + step (staggered leapfrog). + # ================================================================== # + def step(self, fields: tuple, src, dt: float) -> tuple: + """One staggered-leapfrog step (v, q at n+1/2; sigma, p, r at n+1).""" + (v_x, v_z, p, q_x, q_z, s_xx, s_zz, s_xz, r_x, r_z, + m_sxx_x, m_sxz_z, m_sxz_x, m_szz_z, m_qx_x, m_qz_z, m_px, m_pz, + m_vxx, m_vzz, m_vxz, m_vzx) = fields + mode = self.source_mode + if mode not in _SOURCE_MODES: + raise ValueError(f"source_mode must be one of {_SOURCE_MODES}, got {mode!r}") + if self._drag is None or self._drag[0] != dt: + raise RuntimeError(f"step(dt={dt}) needs init_state(dt) first: the drag " + f"operators are computed there, once per run") + cp = self.cpml + ibool = self.mesh.ibool + w = self._w_J + + # ----- pore pressure injected at t^n; the total stress does not move ----- + p = self._pin(p + dt * self._inv_M1 * src[2]) + + # ----- forces at t^n: total-stress divergence and pressure gradient ----- + d_p_dx, d_p_dz = phys_grads(self._D, ibool, self._ixx, self._ixz, self._iex, + self._iez, p) + if cp is not None: + m_sxx_x = cp.a_x * m_sxx_x + cp.b_x * s_xx + m_sxz_z = cp.a_z * m_sxz_z + cp.b_z * s_xz + m_sxz_x = cp.a_x * m_sxz_x + cp.b_x * s_xz + m_szz_z = cp.a_z * m_szz_z + cp.b_z * s_zz + m_px = cp.a_x * m_px + cp.b_x * d_p_dx + m_pz = cp.a_z * m_pz + cp.b_z * d_p_dz + d_p_dx, d_p_dz = d_p_dx + m_px, d_p_dz + m_pz + R_x, R_z = self._weak_div((s_xx + m_sxx_x, s_xz + m_sxz_z), + (s_xz + m_sxz_x, s_zz + m_szz_z)) + else: + R_x, R_z = self._weak_div((s_xx, s_xz), (s_xz, s_zz)) + f_x, f_z = -R_x, -R_z + ifc = self.ifc + if ifc.n_face: + # a fluid's pressure on the solid: -∮ p n φ, n out of the solid + s_x, s_z = ifc.apply_el(p, p) + f_x, f_z = f_x + s_x, f_z + s_z + + # ----- v and the flux copies at every node, the drag integrated exactly ----- + wr, inv_m, qmask = self._wr, self._inv_M, self._qmask + # w (-grad p + r) per copy; nothing moves the copies of a dry element + g_x, g_z = w * (r_x - d_p_dx) * qmask, w * (r_z - d_p_dz) * qmask + if mode == "load": + f_x, f_z = f_x + src[0], f_z + src[1] + r_xs = self._to_slots(g_x) - wr * (f_x * inv_m)[:, None] + r_zs = self._to_slots(g_z) - wr * (f_z * inv_m)[:, None] + if mode == "load": + r_xs = r_xs + self._share * src[3][:, None] + r_zs = r_zs + self._share * src[4][:, None] + qs_x, qs_z = self._to_slots(q_x), self._to_slots(q_z) + new_x, new_z = self._advance_flux(qs_x, qs_z, r_xs, r_zs) + v_x = v_x + (dt * f_x - (wr * (new_x - qs_x)).sum(1)) * inv_m + v_z = v_z + (dt * f_z - (wr * (new_z - qs_z)).sum(1)) * inv_m + q_x, q_z = self._from_slots(new_x, q_x), self._from_slots(new_z, q_z) + if mode == "velocity": + v_x = v_x + dt * inv_m * src[0] + v_z = v_z + dt * inv_m * src[1] + q_x = q_x + (dt * self._inv_Mm * src[3])[ibool] * qmask + q_z = q_z + (dt * self._inv_Mm * src[4])[ibool] * qmask + + # ----- stress, pressure and the viscous memory from the new v, q ----- + d_vx_dx, d_vx_dz, d_vz_dx, d_vz_dz = phys_grads( + self._D, ibool, self._ixx, self._ixz, self._iex, self._iez, v_x, v_z) + if cp is not None: + m_vxx = cp.a_x * m_vxx + cp.b_x * d_vx_dx + m_vzz = cp.a_z * m_vzz + cp.b_z * d_vz_dz + m_vxz = cp.a_z * m_vxz + cp.b_z * d_vx_dz + m_vzx = cp.a_x * m_vzx + cp.b_x * d_vz_dx + d_vx_dx = d_vx_dx + m_vxx + d_vz_dz = d_vz_dz + m_vzz + d_vx_dz = d_vx_dz + m_vxz + d_vz_dx = d_vz_dx + m_vzx + m_qx_x = cp.a_x * m_qx_x + cp.b_x * q_x + m_qz_z = cp.a_z * m_qz_z + cp.b_z * q_z + div_q = self._weak_div((q_x + m_qx_x, q_z + m_qz_z)) + else: + div_q = self._weak_div((q_x, q_z)) + div_v = d_vx_dx + d_vz_dz + alpha = self._alpha + rhs = div_q - assemble(w * alpha * div_v, ibool, p.new_zeros(self.mesh.n_global)) + if ifc.n_face: + rhs = rhs - ifc.apply_ac(v_x, v_z, p) # +∮ ψ v·n, n out of the solid + dp_free = self._pin(dt * self._inv_Mp * rhs) + p = p + dp_free + lam_d, mu = self._lam_d, self.mu_e + dsig_p = alpha * dp_free[ibool] # sigma = sigma' - alpha p + s_xx = s_xx + dt * (lam_d * div_v + 2.0 * mu * d_vx_dx) - dsig_p + s_zz = s_zz + dt * (lam_d * div_v + 2.0 * mu * d_vz_dz) - dsig_p + s_xz = s_xz + dt * mu * (d_vx_dz + d_vz_dx) + if self.viscous_q is not None: + # dr/dt = -(r - (tau_e / tau_s - 1) B q) / tau_s, exact over the step at q^{n+1/2} + tau_e, tau_s = self._tau + decay = math.exp(-dt / tau_s) + gain = (1.0 - decay) * (tau_e / tau_s - 1.0) + b_xx, b_xz, b_zz = self.b_e + r_x = decay * r_x + gain * (b_xx * q_x + b_xz * q_z) + r_z = decay * r_z + gain * (b_xz * q_x + b_zz * q_z) + + return (v_x, v_z, p, q_x, q_z, s_xx, s_zz, s_xz, r_x, r_z, + m_sxx_x, m_sxz_z, m_sxz_x, m_szz_z, m_qx_x, m_qz_z, m_px, m_pz, + m_vxx, m_vzz, m_vxz, m_vzx) + + +@dataclass +class BiotAnisoPermSEM2DDeformed(BiotVSSEM2DDeformed): + """:class:`BiotVSSEM2DDeformed` with an anisotropic permeability: the drag + ``eta kappa^-1`` is the symmetric tensor ``[[b_xx, b_xz], [b_xz, b_zz]]`` (x right, z + down), SPECFEM2D's ``kxx, kxz, kzz``; build it with + :func:`sweep.equations.biot_params.biot_drag_tensor`. The flux copies at a node then + advance with both components coupled through the drag.""" + + tensor_drag: ClassVar[bool] = True + + @classmethod + def build( + cls, + mesh: Quad9Mesh, + H_e: torch.Tensor, + mu_e: torch.Tensor, + C_e: torch.Tensor, + M_e: torch.Tensor, + rho_e: torch.Tensor, + rho_f_e: torch.Tensor, + m_e: torch.Tensor, + b_xx_e: torch.Tensor, + b_xz_e: torch.Tensor, + b_zz_e: torch.Tensor, + *, + dirichlet: tuple = (), + cpml: Optional[SimpleElasticCPMLProfile] = None, + viscous_q: Optional[tuple] = None, + ) -> "BiotAnisoPermSEM2DDeformed": + """Build from the seven Biot coefficients and the three drag components at the + GLL nodes.""" + check_element_arrays(mesh, b_xx_e=b_xx_e, b_xz_e=b_xz_e, b_zz_e=b_zz_e) + return cls._build(mesh, H_e, mu_e, C_e, M_e, rho_e, rho_f_e, m_e, + (b_xx_e, b_xz_e, b_zz_e), dirichlet, cpml, viscous_q) diff --git a/src/sweep/sem/fields.py b/src/sweep/sem/fields.py index 43ded13..7f26bc5 100644 --- a/src/sweep/sem/fields.py +++ b/src/sweep/sem/fields.py @@ -172,15 +172,16 @@ def resolve_field_types(equation, kinds, role): ========================================== ================== ========== Neither exposes a stress receiver. In the velocity-stress family the - stresses DO exist as fields, but they are element-local (``E, ngll, ngll``, - DG-style, so material jumps can live on element faces) while the point - receiver gather works on the global CG nodes -- hence - ``supports_receiver=False`` on ``sigma_*``, a storage-layout limit rather - than a naming one. In the displacement family there is no stress state at - all; it is formed inside the stiffness apply and never stored. - - So FD's ``vx``/``sxx`` names transfer only partially: ``v_x`` exists in the - velocity-stress family, ``sxx`` in neither. + stresses DO exist as fields, element-local (``E, ngll, ngll``, DG-style, + so material jumps can live on element faces), and the receiver gather + reads an element-local field from the receiver's own element (Biot's flux + and stresses are recorded that way); the elastics just do not flag theirs + ``supports_receiver``. In the displacement family there is no + stress state at all; it is formed inside the stiffness apply and never + stored. + + So FD's ``vx``/``sxx`` names transfer only partially on the elastics: + ``v_x`` exists in the velocity-stress family, ``sxx`` in neither. """ attr = "supports_source" if role == "source" else "supports_receiver" specs = list(equation.field_specs) @@ -200,7 +201,10 @@ def resolve_field_types(equation, kinds, role): # equation may override (the fluid-solid form defaults to chi: a # scalar source in the water column). if role == "receiver": - return [s.name for s in avail] + # Biot records its solid velocity by default, as the FD Biot does, + # not all eight of its receiver fields + declared = getattr(equation, "default_receiver_type", None) + return list(declared) if declared else [s.name for s in avail] declared = getattr(equation, "default_source_type", None) return list(declared) if declared else [avail[-1].name] name = getattr(equation, "__name__", type(equation).__name__) diff --git a/src/sweep/sem/prop_base.py b/src/sweep/sem/prop_base.py index 49fb647..27bc7f8 100644 --- a/src/sweep/sem/prop_base.py +++ b/src/sweep/sem/prop_base.py @@ -85,7 +85,7 @@ from .equations.acoustic_potential2d import AcousticPotentialSEM2D from .fields import (MESH_NDIM, format_field_specs, mesh_kind, model_names, resolve_field_types) -from .propagator import SEMPropagator, receiver_gather +from .propagator import SEMPropagator, field_layouts, receiver_gather @@ -332,6 +332,65 @@ def _resolve_vector_source(equation, source_type, ndim): return out, frozenset(axes), None +def _resolve_named_source(equation, source_type): + """``source_type`` on an equation that declares its PropSEM source names + (``propsem_source_names``; Biot) -> ``(names, payload, moment, mode)``. + + A name is ``(kind, values)``: ``"moment"`` adds ``(Mxx, Mzz, Mxz)`` to the + moment tensor; the other kinds add a payload over the source fields -- + ``"velocity"`` injects into the field (the equation's ``source_mode + "velocity"``), ``"force"`` loads its equations (``source_mode "load"``), + ``"injection"`` goes with either. A moment takes the wavelet's integral and + the rest the wavelet itself, so they do not mix in one call, nor do an + injection and a load into the same fields. A name the equation lists in + ``propsem_source_refused`` raises with its reason. + """ + vocab = equation.propsem_source_names + refused = getattr(equation, "propsem_source_refused", {}) + name = equation.__name__ + fields = [s.name for s in equation.field_specs if s.supports_source] + out, payload, moment, kinds, owner = [], None, None, set(), {} + for want in (source_type or equation.default_source_type): + if want in out: + raise ValueError(f"duplicate source_type entry '{want}'") + if want in refused: + raise NotImplementedError( + f"source_type '{want}' on {name}: {refused[want]}") + if want not in vocab: + raise ValueError( + f"Unknown source_type entry '{want}' for {name}. Available " + f"source_type values: {', '.join(vocab)}") + kind, val = vocab[want] + if kind == "moment": + moment = val if moment is None else tuple( + a + b for a, b in zip(moment, val)) + else: + for i, a in enumerate(val): + if a: + if i in owner: + raise ValueError( + f"source_type entries '{owner[i]}' and '{want}' " + f"both drive {fields[i]} on {name}") + owner[i] = want + payload = val if payload is None else tuple( + a + b for a, b in zip(payload, val)) + kinds.add(kind) + out.append(want) + if moment is not None and payload is not None: + raise ValueError( + f"source_type {out} mixes a stress source with an injected field or " + f"a force: the stress source takes the wavelet's integral (a " + f"moment-tensor load), the others the wavelet itself, so they need " + f"two time functions. Run them as two calls and add the records.") + if {"velocity", "force"} <= kinds: + raise ValueError( + f"source_type {out} mixes an injected velocity with a force: one " + f"moves its field after the update, the other loads the equations. " + f"Run them as two calls and add the records.") + mode = "load" if (moment is not None or "force" in kinds) else "velocity" + return out, payload, moment, mode + + def _resolve_snapshot_steps(nt, return_wavefield, snapshot_times, snapshot_interval): """The steps to copy the wavefield at -- ``PropTorch``'s rule: the named @@ -552,7 +611,9 @@ class PropSEM: ``f_dominant`` (the C-PML's frequency), ``density`` (the fixed density of :class:`AcousticPotentialSEM2D`), ``bs_band_mode`` and ``bs_store_dtype_elem`` (what boundary saving keeps of the C-PML band, - and at what precision its element histories on ``impl='c'``). + and at what precision its element histories on ``impl='c'``), and + ``equation_options``, the equation's own ``build`` keywords (e.g. Biot's + ``viscous_q``). """ def __init__( @@ -578,6 +639,7 @@ def __init__( bs_store_dtype_elem=None, source_type=None, receiver_type=None, + equation_options=None, **kwargs, ): if kwargs: @@ -601,7 +663,7 @@ def __init__( # the C-PML builders and the DD strips count the layer in elements self.abcn_elements = self.abcn // self.n_order self._resolve_equation(equation, shape, self.abcn_elements, - source_type, receiver_type) + source_type, receiver_type, equation_options) if any((n - 1) % self.n_order for n in shape): raise ValueError( @@ -648,6 +710,9 @@ def _resolve_memory(self, memory, bs_store_dtype_elem): m = as_memory_strategy(memory) self.memory = Full() if m is None else m self._memory = self._bs_store = self._ckpt_chunk = None + refusal = getattr(self.equation, "bs_refusal", None) + if isinstance(m, BoundarySaving) and refusal: + raise NotImplementedError(refusal) def unsupported(obj, keep): return [f.name for f in dataclasses.fields(obj) @@ -688,7 +753,8 @@ def unsupported(obj, keep): "impl='c' and memory=BoundarySaving().") self._bs_store_elem = bs_store_dtype_elem - def _resolve_equation(self, equation, shape, abcn_elements, source_type, receiver_type): + def _resolve_equation(self, equation, shape, abcn_elements, source_type, receiver_type, + equation_options=None): """Derive and validate everything PropSEM needs from the equation: mesh family vs shape, the model list, what build()'s signature can carry, the calling convention, and the field-role selections. @@ -717,12 +783,21 @@ def _resolve_equation(self, equation, shape, abcn_elements, source_type, receive # the signature HERE, so it is a clear refusal at construction rather # than a TypeError from inside the first forward. import inspect as _inspect - _kw = {q.name for q in - _inspect.signature(equation.build).parameters.values() - if q.kind is q.KEYWORD_ONLY} + _params = _inspect.signature(equation.build).parameters + _kw = {q.name for q in _params.values() if q.kind is q.KEYWORD_ONLY} _split = set(getattr(equation, "cpml_builders", None) or ()) - _unsupported = (_kw - {"cpml", "dirichlet", "domain", "free_top"} - - _split) + # the rest of build()'s keywords are the equation's own options, which + # ``equation_options`` sets (Biot's viscous_q); a required one PropSEM + # cannot supply is refused + _own = _kw - {"cpml", "dirichlet", "domain", "free_top"} - _split + self._eq_options = dict(equation_options or {}) + _bad = sorted(set(self._eq_options) - _own) + if _bad: + raise TypeError( + f"equation_options {_bad} are not options of " + f"{equation.__name__}.build; it takes {sorted(_own) or 'none'}") + _unsupported = {k for k in _own - set(self._eq_options) + if _params[k].default is _params[k].empty} if _unsupported: raise NotImplementedError( f"{equation.__name__}.build takes keyword argument(s) " @@ -749,7 +824,12 @@ def _resolve_equation(self, equation, shape, abcn_elements, source_type, receive # the vector equations (a force with one component per axis) also # take the FD names; see _resolve_vector_source self._src_force_axes, self._src_moment = None, None - if not self.is_acoustic and _stamp is None and n_src == self.ndim: + self._src_payload, self._src_mode = None, None + self._n_src = n_src + if getattr(equation, "propsem_source_names", None): + (self.source_type, self._src_payload, self._src_moment, + self._src_mode) = _resolve_named_source(equation, source_type) + elif not self.is_acoustic and _stamp is None and n_src == self.ndim: (self.source_type, self._src_force_axes, self._src_moment) = _resolve_vector_source( equation, source_type, self.ndim) @@ -1041,7 +1121,7 @@ def _resolve_model_batch(self, models): raise ValueError( f"per-shot models disagree on the batch size: {sorted(sizes)}") B = sizes.pop() - c_max = float(torch.as_tensor(models[0]).max().detach()) # batch-max vp + c_max = self._pml_speed(models) # batch-max vp eqs = [self._build_eq([m[b] if ps(m) else m for m in models], c_max=c_max) for b in range(B)] if self._use_c: @@ -1055,8 +1135,21 @@ def _resolve_model_batch(self, models): def _build_eq(self, models, c_max=None): eq = self._build_eq_raw(models, c_max=c_max) eq.bs_band_mode = self.bs_band_mode + if self._src_mode is not None: + # injected fields, or loads (a force, a moment tensor); see + # _resolve_named_source + eq.source_mode = self._src_mode return eq + def _pml_speed(self, models): + """The speed the C-PML is designed for: the largest ``vp`` (the first + model), or what the equation declares (``pml_speed``; Biot's fast P + wave of its eight coefficients).""" + speed = getattr(self.equation, "pml_speed", None) + if speed is None: + return float(torch.as_tensor(models[0]).max().detach()) + return float(speed(models)) + def _cpml_cached(self, tag, make, *key): """Memoise one C-PML profile on an EXACT key. @@ -1088,7 +1181,7 @@ def _build_eq_raw(self, models, c_max=None): grids = [self._interp_to_gll(m.to(device=self.dev, dtype=self.dtype)) for m in models] if c_max is None: - c_max = float(grids[0].max().detach()) # vp drives the PML profile + c_max = self._pml_speed(grids) # vp drives the PML profile sides = self._absorbing_sides() if self.is_acoustic: dirichlet = self.free_surface @@ -1105,9 +1198,11 @@ def _build_eq_raw(self, models, c_max=None): if self.equation is AcousticPotentialSEM2D: # vp + constant density rho = torch.full_like(grids[0], self._c_density) return AcousticPotentialSEM2D.build( - self.mesh, grids[0], rho, dirichlet=dirichlet, cpml=cpml) + self.mesh, grids[0], rho, dirichlet=dirichlet, cpml=cpml, + **self._eq_options) return self.equation.build(self.mesh, grids[0], - dirichlet=dirichlet, cpml=cpml) + dirichlet=dirichlet, cpml=cpml, + **self._eq_options) # elastic: vp/vs/rho and a NATURAL free surface -- traction-free is what # the weak form does on a side with no absorbing layer, so there is no # dirichlet set to build. ``free_surface`` therefore acts only by @@ -1142,7 +1237,7 @@ def _build_eq_raw(self, models, c_max=None): pml_thickness_elements=self.abcn_elements, c_max=c_max, f_dominant=self.f_dominant), sides, self.abcn_elements, c_max, self.f_dominant) - return self.equation.build(self.mesh, *grids, **kw) + return self.equation.build(self.mesh, *grids, **kw, **self._eq_options) # receiver field indices (1 for acoustic, 3 components for elastic) def _recv_idxs(self, eq): @@ -1195,7 +1290,7 @@ def _default_source(self): """Is the source the equation's default -- the unit vertical force on the elastics, the one source field elsewhere? (What a bare node id means to the DD driver.)""" - if self._src_moment is not None: + if self._src_moment is not None or self._src_payload is not None: return False if self._src_force_axes is None: # scalar / equation-declared return True @@ -1221,7 +1316,16 @@ def _source_vector(self, override=None): f"tensor) source; give a force source_type to use it.") return None if override is not None: - return tuple(float(c) for c in override) + vec = tuple(float(c) for c in override) + if self._src_payload is not None and len(vec) != self._n_src: + own = [s.name for s in self.equation.field_specs + if s.supports_source] + raise ValueError( + f"source_vector on {self.equation.__name__} has one entry " + f"per source field ({', '.join(own)}); got {len(vec)}") + return vec + if self._src_payload is not None: # an equation's named sources + return tuple(self._src_payload) if self._src_force_axes is not None: return tuple(1.0 if a in self._src_force_axes else 0.0 for a in range(nd)) @@ -1627,19 +1731,19 @@ def _forward_eager(self, eq, pershot_eqs, props, acq, *, memory=None, # same two things every time -- the global node ids behind a stamp's element # and the stamp's Lagrange weights). # ------------------------------------------------------------------ # - def _gather_rows(self, mesh, prop): + def _gather_rows(self, mesh, prop, layout="global"): """One shot's ``(rec_G, rec_W)`` -- see ``propagator.receiver_gather``, which this and the C consts builder share.""" - return receiver_gather(mesh, prop.receivers) + return receiver_gather(mesh, prop.receivers, layout) def _on_dev(self, G, W): return G.to(self.dev), W.to(device=self.dev, dtype=self.dtype) - def _gather_batch(self, mesh, props): + def _gather_batch(self, mesh, props, layout="global"): """The same, stacked over shots: ``(B, nrec, ngll^d)``.""" - G, W = zip(*(self._gather_rows(mesh, p) for p in props)) + G, W = zip(*(self._gather_rows(mesh, p, layout) for p in props)) return self._on_dev(torch.stack(G), torch.stack(W)) # ------------------------------------------------------------------ # @@ -1707,6 +1811,11 @@ def _run_batched(self, eq, props, wav_per_shot, recv_idxs, *, patterns = [torch.stack([pl[f] for pl in pat_list]) for f in range(len(pat_list[0]))] # (B, n_global) rec_G, rec_W = self._gather_batch(mesh, props) # (B, nrec, ngll^d) + # an element-local field (Biot's flux, stresses) is gathered from its + # own layout + rec_GI = {ri: rec_G if lay == "global" else + self._gather_batch(mesh, props, lay)[0] + for ri, lay in zip(recv_idxs, field_layouts(eq, recv_idxs))} wav_B = torch.stack(wav_per_shot) # (B, nt) st0 = eq.init_state(dt, dev, dtype) # tuple @@ -1725,7 +1834,8 @@ def advance(n0, k, dyn_in): sval = wav_B[:, n0 + j] src = [pat * sval[:, None] for pat in patterns] d = vstep(d, src) - comps = [(d[ri][ar, rec_G] * rec_W).sum(-1) for ri in recv_idxs] + comps = [(d[ri].reshape(B, -1)[ar, rec_GI[ri]] * rec_W).sum(-1) + for ri in recv_idxs] chunk.append(torch.stack(comps, -1)) # (B, nrec, ncomp) return d, torch.stack(chunk, dim=1) # (B, k, nrec, ncomp) diff --git a/src/sweep/sem/propagator.py b/src/sweep/sem/propagator.py index 67520c7..5f8ee66 100644 --- a/src/sweep/sem/propagator.py +++ b/src/sweep/sem/propagator.py @@ -41,17 +41,19 @@ from .mesh.hex8 import Hex8Mesh from .mesh.hex27 import Hex27Mesh from .basis.quad9 import q2_basis_at, q2_basis_deriv_at -from ._assembly import scatter_point_sources +from ._assembly import element_mean, scatter_point_sources from .fields import resolve_field_types # --------------------------------------------------------------------------- # # Lagrange basis helpers (scalar evaluation at a point in [-1, 1]). # --------------------------------------------------------------------------- # -def receiver_gather(mesh, receivers): +def receiver_gather(mesh, receivers, layout="global"): """``(rec_G, rec_W)`` of one shot: the global node ids behind every receiver's element and that receiver's Lagrange weights, both - ``(nrec, ngll**dim)``. + ``(nrec, ngll**dim)``. With ``layout="element"`` the ids index an + element-local field flattened to ``(E * ngll**dim,)`` instead -- the + receiver's own element's nodes (Biot's flux and stresses). Every receiver samples ``sum_ab w_eta[a] w_xi[b] . field[ibool[e, a, b]]``, so ONE gather plus a weighted sum answers all of them at once -- the @@ -64,8 +66,21 @@ def receiver_gather(mesh, receivers): callers wrote these two lines out (the eager forward, PropSEM's batch gather, SEMCudaPropagator's consts). """ - return (torch.stack([mesh.ibool[r.elem_id].reshape(-1) for r in receivers]), - torch.stack([r.gather_weights() for r in receivers])) + if layout == "element": + m = mesh.ibool[0].numel() + nodes = torch.arange(m, device=mesh.ibool.device) + G = torch.stack([r.elem_id * m + nodes for r in receivers]) + elif layout == "global": + G = torch.stack([mesh.ibool[r.elem_id].reshape(-1) for r in receivers]) + else: + raise ValueError(f"receivers sample a 'global' or an 'element' field, " + f"not {layout!r}") + return G, torch.stack([r.gather_weights() for r in receivers]) + + +def field_layouts(eq, idxs): + """The storage layout of each state field in ``idxs`` (``field_specs``).""" + return [eq.field_specs[i].layout for i in idxs] def lagrange_basis_at(nodes: torch.Tensor, x: float) -> torch.Tensor: @@ -566,13 +581,24 @@ def _stamp_components(stamp, kind, payload): Returns a list aligned with the equation's source fields: ``[scalar]`` for an acoustic body source, ``[F_x, F_z(, F_y)]`` for an elastic body force, the two/three moment-tensor components for ``moment_tensor`` / ``explosion``. + + A ``force`` payload has one entry per source field (Biot's three: ``v_x, + v_z, p``), a 2-D ``moment_tensor`` one ``(Mxx, Mzz, Mxz)`` block per pair + of them; fields past the payload take no load. """ if kind == "scalar": return [stamp.scalar_stamp(payload[0])] if kind == "force": - return list(stamp.force_stamp(*payload)) + # one Lagrange stamp per entry: force_stamp's ``F * outer`` product, for + # any number of source fields + return [stamp.scalar_stamp(f) for f in payload] if kind == "moment_tensor": - return list(stamp.moment_tensor_stamp(payload[0], payload[1], payload[2])) + if len(payload) % 3: + raise ValueError( + f"a moment_tensor payload is (Mxx, Mzz, Mxz) blocks; got " + f"{len(payload)} entries") + return [c for k in range(0, len(payload), 3) + for c in stamp.moment_tensor_stamp(*payload[k:k + 3])] if kind == "explosion": # isotropic Mxx = Mzz = M0, Mxz = 0 m0 = payload[0] return list(stamp.moment_tensor_stamp(m0, m0, 0.0)) @@ -583,15 +609,42 @@ def _stamp_components(stamp, kind, payload): _DERIVATIVE_KINDS = ("moment_tensor", "explosion") -def _stamp_loads(stamp, kind, payload): +def _stamp_loads(stamp, kind, payload, alive=None, names=None): """``[(elem_id, components), ...]`` of one point source: its element for a :class:`PointStamp`, every element of an :class:`EdgeStamp`, each - weighted by its share.""" + weighted by its share. + + ``alive`` (an equation's ``source_field_elements``: per source field, the + ``(E,)`` bool of the elements it lives in) keeps each field's load to + those elements, the shares renormalised over them -- a source on an + interface between two media loads each field from its own side. A source + left loading nothing is refused (``names``: the source fields).""" parts = getattr(stamp, "parts", None) - if parts is None: - return [(stamp.elem_id, _stamp_components(stamp, kind, payload))] - return [(p.elem_id, [c * w for c in _stamp_components(p, kind, payload)]) - for p, w in zip(parts, stamp.weights)] + weights = [1.0] if parts is None else list(stamp.weights) + parts = [stamp] if parts is None else parts + comps = [_stamp_components(p, kind, payload) for p in parts] + if alive is None: + if len(parts) == 1: + return [(parts[0].elem_id, comps[0])] + return [(p.elem_id, [c * w for c in cs]) + for p, w, cs in zip(parts, weights, comps)] + out = [[] for _ in parts] + loaded = False + for f in range(len(comps[0])): + keep = [bool(alive[f][p.elem_id]) if f < len(alive) else True for p in parts] + total = sum(w for w, k in zip(weights, keep) if k) + for i, k in enumerate(keep): + c = comps[i][f] * (weights[i] / total) if k else comps[i][f] * 0.0 + loaded = loaded or bool(k and c.abs().max() > 0) + out[i].append(c) + if not loaded and any(bool(c.abs().max() > 0) for cs in comps for c in cs): + fields = [n for f, n in enumerate(names or ()) + if any(bool(cs[f].abs().max() > 0) for cs in comps if f < len(cs))] + raise ValueError( + f"a {kind!r} source in element(s) {[p.elem_id for p in parts]} loads " + f"{fields or 'its fields'}, none of which lives there (a fluid element " + f"carries no solid velocity or stress, a dry one no pore pressure or flux)") + return [(p.elem_id, cs) for p, cs in zip(parts, out)] # --------------------------------------------------------------------------- # @@ -732,8 +785,12 @@ def build( default = (1.0,) + (0.0,) * (n_src_fields - 1) source_payloads = [default] * len(source_xz) # a moment tensor differentiates the basis, so on an element edge it - # is stamped in every element there (EdgeStamp); the rest locate once - sources = [(build_moment_stamp if kind in _DERIVATIVE_KINDS + # is stamped in every element there (EdgeStamp); the rest locate once -- + # unless the fields live in different media (source_field_elements): + # then every source meets all the elements at its point, to take each + # field's load from the ones that carry it + edge_all = getattr(eq, "source_field_elements", None) is not None + sources = [(build_moment_stamp if edge_all or kind in _DERIVATIVE_KINDS else build_point_stamp)(eq.mesh, *xz) for xz, kind in zip(source_xz, source_kinds)] receivers = [build_point_stamp(eq.mesh, *xz) for xz in receiver_xz] @@ -742,6 +799,14 @@ def build( receiver_type=resolve_field_types( type(eq), receiver_type, "receiver")) + def _source_media(self): + """``(alive, names)`` for :func:`_stamp_loads`: the equation's + ``source_field_elements`` (``None`` when every field lives everywhere) + and its source field names.""" + alive = getattr(self.eq, "source_field_elements", None) + names = [s.name for s in self.eq.field_specs if s.supports_source] + return alive, names + # ---- source spatial pattern (assembled ONCE, scaled per step) ------- # def _build_source_patterns(self, n_global, device, dtype): """Assemble the constant per-field source spatial pattern. @@ -762,9 +827,16 @@ def _build_source_patterns(self, n_global, device, dtype): # the positional order is already the equation's own field order. patterns = [torch.zeros(n_global, dtype=dtype, device=device) for _ in range(n_src)] + alive, names = self._source_media() for stamp, kind, payload in zip(self.sources, self.source_kinds, self.source_payloads): - for e, comps in _stamp_loads(stamp, kind, payload): + for e, comps in _stamp_loads(stamp, kind, payload, alive, names): + if len(comps) > n_src: + # zip would drop the extra components without a word + raise ValueError( + f"a {kind!r} source with payload {tuple(payload)} has " + f"{len(comps)} components; {type(eq).__name__} has " + f"{n_src} source fields") g = mesh.ibool[e].reshape(-1) for pat, comp in zip(patterns, comps): pat.index_add_(0, g, comp.reshape(-1)) @@ -787,12 +859,17 @@ def _build_per_source_patterns(self, device, dtype): mesh = self.eq.mesh n_src = sum(s.supports_source for s in self.eq.field_specs) G_rows, W_rows = [], [[] for _ in range(n_src)] + alive, names = self._source_media() for stamp, kind, payload in zip(self.sources, self.source_kinds, self.source_payloads): - loads = _stamp_loads(stamp, kind, payload) + loads = _stamp_loads(stamp, kind, payload, alive, names) G_rows.append(torch.cat([mesh.ibool[e].reshape(-1) for e, _ in loads])) for f, rows in enumerate(W_rows): - rows.append(torch.cat([c[f].reshape(-1) for _, c in loads])) + # a payload may cover the leading source fields only (a moment + # tensor on Biot's v rows, none on p): the rest take no load + rows.append(torch.cat([ + (c[f] if f < len(c) else torch.zeros_like(c[0])).reshape(-1) + for _, c in loads])) width = max(g.numel() for g in G_rows) if any(g.numel() != width for g in G_rows): # an edge source spans several elements: pad the others with zero @@ -856,6 +933,7 @@ def forward( device, dtype = self._dev_dtype() recv = self._recv_components() rec_G, rec_W = self._receiver_gather(device, dtype) + layouts = field_layouts(self.eq, recv) source_spec, like = self._source_term(per_source, device, dtype) step_fn = self._step_fn(compile_step) @@ -865,6 +943,11 @@ def forward( snap_at is not None, compile_step, bs_store_dtype) n_rec = len(self.receivers) + # an element-local field (Biot's flux, stresses) is gathered from its + # own layout + rec_GI = {i: rec_G if lay == "global" else + receiver_gather(self.eq.mesh, self.receivers, lay)[0].to(device) + for i, lay in zip(recv, layouts)} if n_rec else {} state = self.eq.init_state(self.dt, device, dtype) # positional tuple seis = {i: torch.zeros(n_rec, self.nt, dtype=dtype, device=device) for i in recv} @@ -879,20 +962,22 @@ def forward( if n_rec > 0: for i in recv: - seis[i][:, n] = (state[i][rec_G] * rec_W).sum(dim=1) + seis[i][:, n] = (state[i].reshape(-1)[rec_GI[i]] + * rec_W).sum(dim=1) if snap_at is not None and n in snap_at: - for i in recv: - snaps[i].append(state[i].clone()) + for i, lay in zip(recv, layouts): + snaps[i].append(state[i].clone() if lay == "global" else + element_mean(self.eq.mesh, state[i])) # (nrec, nt) per component -> (nt, nrec, ncomp), the FD layout. out = torch.stack([seis[i] for i in recv], -1).transpose(0, 1) if snap_at is None: return out if not snap_at: # snapshot_times=[]: none asked for - return out, state[recv[0]].new_zeros( - (0, len(recv)) + tuple(state[recv[0]].shape)) - # (nsnap, ncomp, n_global), at the steps of snap_at in ascending order + return out, state[recv[0]].new_zeros((0, len(recv), self.eq.mesh.n_global)) + # (nsnap, ncomp, n_global), at the steps of snap_at in ascending order; + # an element-local field as the |J|-weighted mean at each shared node return out, torch.stack([torch.stack(snaps[i], 0) for i in recv], 1) def _snapshot_steps(self, every, times): diff --git a/test/sem/specfem_compare/make_specfem_poro_cases.py b/test/sem/specfem_compare/make_specfem_poro_cases.py new file mode 100644 index 0000000..6eceb2f --- /dev/null +++ b/test/sem/specfem_compare/make_specfem_poro_cases.py @@ -0,0 +1,149 @@ +"""SPECFEM2D input directories for ``run_compare_poro_features.py`` and +``run_compare_poro_coupling.py``, built from SPECFEM2D's own examples: + +* on a 2000 m x 1600 m box of 5 m Quad9 elements (free top, Stacey elsewhere, echoes after + the 0.25 s record), from ``poroelastic_semi_infinite_homogeneous``: ``force_z`` (a vertical + point force), ``aniso_k`` (an anisotropic permeability), ``visc`` and ``visc_atten`` (the + viscous drag, then with ``ATTENUATION_PORO_FLUID_PART``) -- the source 600 m deep, two of + the four receivers on its grid lines --, ``pe_layers`` (porous rock over a dry elastic one) + and ``ape_layers`` (water, porous rock, elastic rock); +* ``pa_example``: the ``poroelastic_acoustic`` example as it is, serial and without images, + and ``pa_example_fine`` on elements half as large. + +Run ``xmeshfem2D`` and then ``xspecfem2D`` in each ``/``. + + python make_specfem_poro_cases.py /EXAMPLES +""" +import re +import shutil +import sys +from pathlib import Path + +ex, out = Path(sys.argv[1]), Path(sys.argv[2]) + +LX, LZ, H = 2000.0, 1600.0, 5.0 +XS, ZS = 1000.0, 1000.0 # z up from the bottom: 600 m deep +RECS = [(1100.0, 1100.0), (1100.0, 1000.0), (1000.0, 1100.0), (1150.0, 950.0)] +ROCK = dict(rho_s=2650.0, rho_f=880.0, phi=0.1, tort=2.0, k_s=12.2e9, k_f=1.985e9, + k_fr=9.6e9, mu_fr=5.1e9) +SHALE = "1 2500.d0 4000.d0 2300.d0 0 0 9999 9999 0 0 0 0 0 0" # rho vp vs +WATER = "1 1020.d0 1500.d0 0.d0 0 0 9999 9999 0 0 0 0 0 0" + +# name: (source_type, angle, eta, (kxx, kxz, kzz) in SPECFEM's z-up axes, (Q0, f0) or None) +FEATURES = { + "force_z": (1, 0.0, 0.0, (1e-9, 0.0, 1e-9), None), + "aniso_k": (2, None, 1e-3, (2e-9, 5e-10, 5e-10), None), + "visc": (2, None, 1e-3, (1e-9, 0.0, 1e-9), None), + "visc_atten": (2, None, 1e-3, (1e-9, 0.0, 1e-9), (5.0, 30.0)), +} +# name: (layers top-down as (material line, thickness in m), source type, source (x, z up), +# receivers (x, z up), steps of 1e-4 s); the porous rock with eta = 1e-3 +LAYERS = { + "pe_layers": (["rock", 800.0, SHALE, 800.0], 2, (1000.0, 1000.0), + [(1100.0, 1100.0), (1100.0, 850.0), (1050.0, 750.0), (1150.0, 650.0)], 2500), + "ape_layers": ([WATER, 300.0, "rock", 500.0, SHALE, 800.0], 1, (1000.0, 1450.0), + [(1100.0, 1400.0), (1100.0, 1250.0), (1050.0, 1000.0), (1150.0, 700.0)], + 5000), +} + + +def setp(text, key, value): + pat = re.compile(rf"^({re.escape(key)}\s*=\s*)(\S+)", re.M) + assert pat.search(text), key + return pat.sub(lambda m: m.group(1) + str(value), text, count=1) + + +def quiet(par): + for k, v in (("NPROC", 1), ("output_color_image", ".false."), + ("NTSTEP_BETWEEN_OUTPUT_IMAGES", 100000), ("seismotype", 2)): + par = setp(par, k, v) + return par + + +def rock(eta, k=(1e-9, 0.0, 1e-9)): + r = ROCK + line = (f"3 {r['rho_s']}d0 {r['rho_f']}d0 {r['phi']}d0 {r['tort']} {k[0]:.6e} {k[1]:.6e} " + f"{k[2]:.6e} {r['k_s']:.6e} {r['k_f']:.6e} {r['k_fr']:.6e} {eta:.6e} {r['mu_fr']:.6e} " + f"9999999999") + return line.replace("e+", "d+").replace("e-", "d-") + + +def write(name, par, src): + d = out / name / "DATA" + d.mkdir(parents=True, exist_ok=True) + (d / "Par_file").write_text(par) + (d / "SOURCE").write_text(src) + (d / "interfaces.dat").write_text( + f"2\n2\n 0 0\n {LX:.0f} 0\n2\n 0 {LZ:.0f}\n {LX:.0f} {LZ:.0f}\n{int(LZ / H)}\n") + print("wrote", d) + + +def box(par, recs): + """The 2000 m x 1600 m box and one receiver set per receiver.""" + par = setp(par, "xmax", f"{LX}d0") + par = setp(par, "nx", int(LX / H)) + block = "".join(f"nrec = 1\nxdeb = {x}\nzdeb = {z}\nxfin = {x}\nzfin = {z}\n" + f"record_at_surface_same_vertical = .false.\n" for x, z in recs) + par = setp(par, "nreceiversets", len(recs)) + return re.sub(r"^nrec\s*=.*\n^xdeb\s*=.*\n^zdeb\s*=.*\n^xfin\s*=.*\n^zfin\s*=.*\n" + r"^record_at_surface_same_vertical\s*=.*\n", lambda m: block, par, count=1, + flags=re.M) + + +base = ex / "poroelastic_semi_infinite_homogeneous" / "DATA" +par0, src0 = quiet((base / "Par_file").read_text()), (base / "SOURCE").read_text() +nx, nz = int(LX / H), int(LZ / H) + +for name, (stype, ang, eta, k, atten) in FEATURES.items(): + par = box(par0, RECS) + par = setp(par, "ATTENUATION_PORO_FLUID_PART", ".true." if atten else ".false.") + if atten: + par = setp(par, "Q0_poroelastic", atten[0]) + par = setp(par, "freq0_poroelastic", atten[1]) + par = re.sub(r"^1 3 .*$", lambda m: "1 " + rock(eta, k), par, count=1, flags=re.M) + par = re.sub(r"^1 200 1 200 1$", f"1 {nx} 1 {nz} 1", par, count=1, flags=re.M) + src = setp(setp(setp(src0, "xs", XS), "zs", ZS), "source_type", stype) + if stype == 1: + src = setp(src, "anglesource", ang) + write(name, par, src) + +for name, (layers, stype, (xs, zs), recs, nstep) in LAYERS.items(): + mats = [rock(1e-3) if m == "rock" else m for m in layers[::2]] + par = setp(box(par0, recs), "NSTEP", nstep) + par = setp(par, "ATTENUATION_PORO_FLUID_PART", ".false.") + par = setp(par, "nbmodels", len(mats)) + par = re.sub(r"^1 3 .*$", lambda m: "\n".join(f"{i + 1} {line}" for i, line in + enumerate(mats)), par, count=1, flags=re.M) + regions, top = [], nz # element rows count up from the bottom + for i, thick in enumerate(layers[1::2]): + n = int(round(thick / H)) + regions.append(f"1 {nx} {top - n + 1} {top} {i + 1}") + top -= n + assert top == 0, name + par = setp(par, "nbregions", len(mats)) + par = re.sub(r"^1 200 1 200 1$", lambda m: "\n".join(regions), par, count=1, flags=re.M) + src = setp(setp(setp(src0, "xs", xs), "zs", zs), "source_type", stype) + if stype == 1: + src = setp(src, "anglesource", 0.0) + write(name, par, src) + +# the poroelastic_acoustic example, and on elements half as large +pa = ex / "poroelastic_acoustic" / "DATA" +for name, scale in (("pa_example", 1), ("pa_example_fine", 2)): + d = out / name / "DATA" + d.mkdir(parents=True, exist_ok=True) + par = quiet((pa / "Par_file").read_text()) + itf = (pa / "interfaces.dat").read_text() + if scale > 1: + par = setp(par, "NSTEP", 2500 * scale) + par = setp(par, "DT", f"{4.2 / scale}d-4") + par = setp(par, "nx", 220 * scale) + par = re.sub(r"^1 220 1 110 1$", f"1 {220 * scale} 1 {110 * scale} 1", par, count=1, + flags=re.M) + par = re.sub(r"^1 220 111 220 2$", f"1 {220 * scale} {110 * scale + 1} {220 * scale} 2", + par, count=1, flags=re.M) + itf = re.sub(r"^\s*110\s*$", str(110 * scale), itf, flags=re.M) + (d / "Par_file").write_text(par) + (d / "interfaces.dat").write_text(itf) + shutil.copy(pa / "SOURCE", d / "SOURCE") + print("wrote", d) diff --git a/test/sem/specfem_compare/run_compare_poro.py b/test/sem/specfem_compare/run_compare_poro.py new file mode 100644 index 0000000..7d2e2ff --- /dev/null +++ b/test/sem/specfem_compare/run_compare_poro.py @@ -0,0 +1,116 @@ +"""BiotVSSEM2DDeformed against SPECFEM2D's ``poroelastic_semi_infinite_homogeneous`` example. + +The example (Morency & Tromp 2008, Fig. 4): a homogeneous inviscid poroelastic square, +1000 m x 1000 m in 5 m Quad9 elements, free top and Stacey elsewhere, an explosion at +(500, 300) and a receiver at (600, 400) (z up), a Gaussian moment history of f0 = 30 Hz, +dt = 1e-4 and 2500 steps. Its reference solid velocities ``REF_SEIS/AA.S0001.BX[XZ].semv`` +ship with SPECFEM2D. The record ends before any boundary echo reaches the receiver. + +SPECFEM2D loads a poroelastic moment tensor as ``+M : grad(phi)`` on both the solid and the +fluid equation of Morency & Tromp's (u_s, w) form, i.e. the same vector on the +total-momentum and the relative-flow equation. In PropSEM's names that is a total-stress +source ``["sxx", "szz"]`` of history ``-M(t)`` plus a pore-pressure source ``["p"]`` of +``+M(t)`` (both the injected quantity, FD's sign); PropSEM takes their rates. + + PYTHONPATH=src python test/sem/specfem_compare/run_compare_poro.py \\ + --ref /EXAMPLES/poroelastic_semi_infinite_homogeneous/REF_SEIS --out +""" + +from __future__ import annotations + +import argparse +import math +from pathlib import Path + +import numpy as np +import torch + +from sweep.equations.biot import biot_from_rock, biot_wave_velocities +from sweep.sem import BiotVSSEM2DDeformed, PropSEM + +# Par_file, model 1: rho_s rho_f phi tort kxx kxz kzz K_s K_f K_fr eta_f mu_fr +ROCK = dict(k_s=12.2e9, rho_s=2650.0, k_f=1.985e9, rho_f=880.0, eta=0.0, k_d=9.6e9, + mu=5.1e9, phi=0.1, kappa=1e-9, tortuosity=2.0) +L, H_EL, DT, NT, F0, FACTOR = 1000.0, 5.0, 1e-4, 2500, 30.0, 1e10 +SRC, REC = (500.0, 300.0), (600.0, 400.0) # SPECFEM (x, z), z up from the bottom + + +def specfem_stf(t): + """SPECFEM2D's time_function_type 3: ``-factor * (-exp(-a t^2) / (2 a))``, a = (pi f0)^2, + with its t0 = 1.2 / f0 already taken off ``t``.""" + a = (math.pi * F0) ** 2 + return FACTOR * np.exp(-a * t**2) / (2.0 * a) + + +def run(device, dtype, n_order=4): + dh = H_EL / n_order + n = int(round(L / dh)) + 1 + models = biot_from_rock(**ROCK) + t0 = 1.2 / F0 + t = np.arange(NT) * DT + m_hist = specfem_stf(t - t0) + rate = torch.from_numpy(np.gradient(m_hist, t)).to(dtype) + # grid indices [ix, iz], z down from the free top + src = np.array([[SRC[0] / dh, (L - SRC[1]) / dh]]) + rec = np.array([[[REC[0] / dh, (L - REC[1]) / dh]]]) + coef = [torch.full((n, n), float(c), dtype=dtype) for c in models] + out = 0.0 + for names, sign in ((["sxx", "szz"], -1.0), (["p"], 1.0)): + prop = PropSEM(BiotVSSEM2DDeformed, shape=(n, n), dh=dh, dt=DT, n_order=n_order, + free_surface=True, f_dominant=F0, dev=device, impl="eager", + dtype=dtype, source_type=names, receiver_type=["vx", "vz"]) + with torch.no_grad(): + out = out + prop(sign * rate, src, rec, models=coef)[0, :, 0].double().cpu().numpy() + vx, vz_down = out[:, 0], out[:, 1] + # a velocity sample n is at (n + 1) dt for an integrated source (see test_biot_sem.py); + # SPECFEM sample m is at m dt - t0 + t_sem = (np.arange(NT) + 1) * DT - t0 + return t_sem, vx, -vz_down, models + + +def read_semv(path): + d = np.loadtxt(path) + return d[:, 0], d[:, 1] + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--ref", type=Path, required=True) + ap.add_argument("--out", type=Path, required=True) + ap.add_argument("--device", default="cuda" if torch.cuda.is_available() else "cpu") + ap.add_argument("--float64", action="store_true") + args = ap.parse_args() + args.out.mkdir(parents=True, exist_ok=True) + dtype = torch.float64 if args.float64 else torch.float32 + + t_sem, vx, vz, models = run(args.device, dtype) + vf, vsl, vs = biot_wave_velocities([float(m) for m in models]) + print(f"fast P {vf:.1f} slow P {vsl:.1f} S {vs:.1f} m/s") + np.savez(args.out / "sem_biot_specfem_example.npz", t=t_sem, vx=vx, vz=vz) + + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + fig, axes = plt.subplots(2, 1, figsize=(8, 6), sharex=True) + for ax, comp, ours in zip(axes, ("BXX", "BXZ"), (vx, vz)): + t_ref, ref = read_semv(args.ref / f"AA.S0001.{comp}.semv") + on_ref = np.interp(t_ref, t_sem, ours) + err = np.linalg.norm(on_ref - ref) / np.linalg.norm(ref) + corr = float(on_ref @ ref) / float(np.linalg.norm(on_ref) * np.linalg.norm(ref)) + print(f"{comp}: rel L2 {err:.4f} corr {corr:.5f} " + f"amp {np.linalg.norm(on_ref) / np.linalg.norm(ref):.4f}") + ax.plot(t_ref, ref, "k", lw=1.6, label="SPECFEM2D") + ax.plot(t_ref, on_ref, "r--", lw=1.0, label="SEM Biot") + ax.plot(t_ref, 5 * (on_ref - ref), color="0.6", lw=0.8, label="5x diff") + ax.set_ylabel(f"{'vx' if comp == 'BXX' else 'vz (up)'} (m/s)") + ax.text(0.01, 0.92, f"({'ab'[comp == 'BXZ']}) rel. L2 {err:.3f}", transform=ax.transAxes) + axes[0].legend(loc="upper right", fontsize=8) + axes[1].set_xlabel("time (s)") + fig.tight_layout() + fig.savefig(args.out / "sem_biot_vs_specfem2d.png", dpi=220) + print("wrote", args.out / "sem_biot_vs_specfem2d.png") + + +if __name__ == "__main__": + main() diff --git a/test/sem/specfem_compare/run_compare_poro_coupling.py b/test/sem/specfem_compare/run_compare_poro_coupling.py new file mode 100644 index 0000000..643b811 --- /dev/null +++ b/test/sem/specfem_compare/run_compare_poro_coupling.py @@ -0,0 +1,158 @@ +"""BiotVSSEM2DDeformed with fluid and dry elastic media against SPECFEM2D: + +* ``pa_example``: SPECFEM2D's ``poroelastic_acoustic`` example (Morency & Tromp 2008, Fig. 8) + -- water over an inviscid porous half-space, a pressure source in the water --, and + ``pa_example_fine``, the same on elements half as large; +* ``pe_layers``: a viscous porous rock over a dry elastic one, an explosion in the rock; +* ``ape_layers``: water, the porous rock, the elastic rock; a pressure source in the water, + +the last two on the box of ``run_compare_poro_features.py`` (cases from +``make_specfem_poro_cases.py``). Velocities only (SPECFEM2D's ``seismotype = 2``): in the water +the fluid's, our ``q``; elsewhere the solid's, our ``v``. The source time function is +SPECFEM2D's own (``plot_source_time_function.txt``). + +SPECFEM2D's z points up (its BXZ is -v_z here). A force in an acoustic element is its +pressure source, ``chi'' = kappa div(grad(chi) / rho) + stf delta`` with ``p = -chi''``: +our pore-pressure injection of the wavelet ``-d stf / dt``. Its moment tensor +``+M : grad(phi)`` in a porous element is our total stress ``-M`` plus pore pressure ``+M``. + + PYTHONPATH=src python test/sem/specfem_compare/run_compare_poro_coupling.py \\ + --runs --out [--case pa_example ...] +""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +import numpy as np +import torch + +from sweep.equations.biot import biot_from_rock +from sweep.equations.biot_params import biot_elastic, biot_fluid +from sweep.sem import BiotVSSEM2DDeformed, PropSEM + +ROCK = biot_from_rock(k_s=12.2e9, rho_s=2650.0, k_f=1.985e9, rho_f=880.0, eta=1e-3, k_d=9.6e9, + mu=5.1e9, phi=0.1, kappa=1e-9, tortuosity=2.0) +PA_ROCK = biot_from_rock(k_s=16.0554e9, rho_s=2500.0, k_f=2.295e9, rho_f=1020.0, eta=0.0, + k_d=10.0e9, mu=9.63342e9, phi=0.4, kappa=1e-9, tortuosity=2.0) +WATER = biot_fluid(1020.0 * 1500.0**2, 1020.0) +SHALE = biot_elastic(4000.0, 2300.0, 2500.0) + +# name: box (LX, LZ, element), f0, dt, nt, media top-down (coefficients, thickness), +# source kind, source (x, z up), receivers (x, z up) +CASES = { + "pa_example": ((4800.0, 4800.0, 4800.0 / 220), 15.0, 4.2e-4, 2500, + [(WATER, 2400.0), (PA_ROCK, 2400.0)], "pressure", (1600.0, 2900.0), + [(2000.0, 2934.0), (2000.0, 1867.0)]), + # the example on elements half as large: its inviscid slow wave (1060 m/s) has ~1.2 + # elements per wavelength at 40 Hz on the example's own mesh + "pa_example_fine": ((4800.0, 4800.0, 4800.0 / 440), 15.0, 2.1e-4, 5000, + [(WATER, 2400.0), (PA_ROCK, 2400.0)], "pressure", (1600.0, 2900.0), + [(2000.0, 2934.0), (2000.0, 1867.0)]), + "pe_layers": ((2000.0, 1600.0, 5.0), 30.0, 1e-4, 2500, [(ROCK, 800.0), (SHALE, 800.0)], + "explosion", (1000.0, 1000.0), + [(1100.0, 1100.0), (1100.0, 850.0), (1050.0, 750.0), (1150.0, 650.0)]), + "ape_layers": ((2000.0, 1600.0, 5.0), 30.0, 1e-4, 5000, + [(WATER, 300.0), (ROCK, 500.0), (SHALE, 800.0)], "pressure", + (1000.0, 1450.0), + [(1100.0, 1400.0), (1100.0, 1250.0), (1050.0, 1000.0), (1150.0, 700.0)]), +} + + +def layered(mesh, media, dtype): + """GLL-native coefficients of horizontal layers, ``media`` top-down as + ``(coefficients, thickness)``; the last one fills the rest (the C-PML below).""" + zc = mesh.z_e.mean((1, 2)) + edges = torch.tensor(np.cumsum([t for _, t in media])[:-1], dtype=zc.dtype, device=zc.device) + k = torch.bucketize(zc, edges) + E, n = mesh.n_elements, mesh.ngll + return [torch.stack([torch.as_tensor(m[i], dtype=dtype) for m, _ in media]).to(zc.device)[k] + .view(E, 1, 1).expand(E, n, n).contiguous() for i in range(8)] + + +def in_water(case, z_up): + (_, lz, _), media = CASES[case][0], CASES[case][4] + return media[0][0] is WATER and lz - z_up < media[0][1] + + +def run_sem(case, stf, device, dtype, n_order=4): + (lx, lz, h), f0, dt, nt, media, kind, (xs, zs), recs = CASES[case] + dh = h / n_order + nz, nx = int(round(lz / dh)) + 1, int(round(lx / dh)) + 1 + kw = dict(shape=(nz, nx), dh=dh, dt=dt, n_order=n_order, free_surface=True, f_dominant=f0, + dev=device, impl="eager", dtype=dtype, receiver_type=["vx", "vz", "qx", "qz"]) + src = np.array([[xs / dh, (lz - zs) / dh]]) + rec = np.array([[[x / dh, (lz - z) / dh] for x, z in recs]]) + rate = np.gradient(stf, dt) + if kind == "explosion": + calls = ((["sxx", "szz"], -rate), (["p"], rate)) + else: + calls = ((["p"], -rate),) + out, models = 0.0, None + for names, w in calls: + prop = PropSEM(BiotVSSEM2DDeformed, source_type=names, **kw) + if models is None: + models = layered(prop.mesh, media, dtype) + with torch.no_grad(): + out = out + prop(torch.from_numpy(w).to(dtype), src, rec, + models=models)[0].double().cpu().numpy() + # an integrated source's velocity samples sit at (n + 1) dt + traces = np.stack([out[:, r, 2:4] if in_water(case, z) else out[:, r, 0:2] + for r, (_, z) in enumerate(recs)], 1) + traces[:, :, 1] *= -1.0 # z down -> SPECFEM's up + return (np.arange(nt) + 1.0) * dt, traces + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--runs", type=Path, required=True) + ap.add_argument("--out", type=Path, required=True) + ap.add_argument("--case", nargs="*", default=list(CASES)) + ap.add_argument("--device", default="cuda" if torch.cuda.is_available() else "cpu") + args = ap.parse_args() + args.out.mkdir(parents=True, exist_ok=True) + + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + for case in args.case: + run = args.runs / case / "OUTPUT_FILES" + t_stf, stf = np.loadtxt(run / "plot_source_time_function.txt", unpack=True) + t_sem, sem = run_sem(case, stf, args.device, torch.float32) + t_sem = t_sem + t_stf[0] # SPECFEM's time axis, t - t0 + np.savez(args.out / f"sem_{case}.npz", t=t_sem, sem=sem) + recs = CASES[case][7] + fig, axes = plt.subplots(len(recs), 2, figsize=(10, 2.3 * len(recs)), sharex=True, + squeeze=False) + for r, (x, z) in enumerate(recs): + med = "water" if in_water(case, z) else "solid" + for k, comp in enumerate(("BXX", "BXZ")): + ax = axes[r, k] + t_ref, ref = np.loadtxt(run / f"AA.S{r + 1:04d}.{comp}.semv", unpack=True) + ours = np.interp(t_ref, t_sem, sem[:, r, k]) + if not (np.linalg.norm(ref) > 0 and np.linalg.norm(ours) > 0): + print(f"{case} rec{r + 1} ({med}) {comp}: no arrival in the record", flush=True) + continue + err = np.linalg.norm(ours - ref) / np.linalg.norm(ref) + corr = float(ours @ ref) / float(np.linalg.norm(ours) * np.linalg.norm(ref)) + amp = np.linalg.norm(ours) / np.linalg.norm(ref) + print(f"{case} rec{r + 1} ({med}) {comp}: rel L2 {err:.4f} corr {corr:+.5f} " + f"amp {amp:.4f}", flush=True) + ax.plot(t_ref, ref, "k", lw=1.4, label="SPECFEM2D") + ax.plot(t_ref, ours, "r--", lw=0.9, label="SEM") + ax.plot(t_ref, 5 * (ours - ref), color="0.6", lw=0.7, label="5x diff") + ax.set_title(f"({'abcdefgh'[2 * r + k]}) rec {r + 1} {med} {comp}, " + f"rel. L2 {err:.4f}", fontsize=8) + axes[0, 0].legend(frameon=False, fontsize=7) + for ax in axes[-1]: + ax.set_xlabel("time (s)") + fig.tight_layout() + fig.savefig(args.out / f"sem_vs_specfem2d_{case}.png", dpi=220) + plt.close(fig) + print("wrote", args.out / f"sem_vs_specfem2d_{case}.png", flush=True) + + +if __name__ == "__main__": + main() diff --git a/test/sem/specfem_compare/run_compare_poro_features.py b/test/sem/specfem_compare/run_compare_poro_features.py new file mode 100644 index 0000000..3b01c57 --- /dev/null +++ b/test/sem/specfem_compare/run_compare_poro_features.py @@ -0,0 +1,137 @@ +"""BiotVSSEM2DDeformed / BiotAnisoPermSEM2DDeformed against SPECFEM2D runs of the poroelastic +features: a point force (``force_z``), an anisotropic permeability (``aniso_k``) and the +viscous attenuation (``visc``, ``visc_atten``), on the box ``make_specfem_poro_cases.py`` +builds -- 2000 m x 1600 m of 5 m Quad9 elements, free top, the source 600 m deep, four +receivers (two on the source's grid lines), a Gaussian of f0 = 30 Hz, dt = 1e-4, 2500 steps. +Solid velocity only: SPECFEM2D computes no pressure in poroelastic media. + +SPECFEM2D's z points up: its BXZ is -v_z here, and its permeability's kxz flips sign. Its +poroelastic force loads the solid and the fluid equations alike (our ``fz``) and, since its +May 2018 sign fix (``sourcearray = (+sin, -cos) hlagrange``), points down at angle 0 for a +positive source time function; its moment tensor ``+M : grad(phi)`` on both is our total +stress ``-M`` plus pore pressure ``+M``. + + PYTHONPATH=src python test/sem/specfem_compare/run_compare_poro_features.py \\ + --runs --out [--case force_z ...] +""" + +from __future__ import annotations + +import argparse +import math +from pathlib import Path + +import numpy as np +import torch + +from sweep.equations.biot import biot_from_rock +from sweep.equations.biot_params import biot_drag_tensor +from sweep.sem import BiotAnisoPermSEM2DDeformed, BiotVSSEM2DDeformed, PropSEM + +LX, LZ, H_EL, DT, NT, F0, FACTOR = 2000.0, 1600.0, 5.0, 1e-4, 2500, 30.0, 1e10 +SRC = (1000.0, 1000.0) # SPECFEM (x, z), z up +RECS = [(1100.0, 1100.0), (1100.0, 1000.0), (1000.0, 1100.0), (1150.0, 950.0)] +ROCK = dict(k_s=12.2e9, rho_s=2650.0, k_f=1.985e9, rho_f=880.0, k_d=9.6e9, mu=5.1e9, phi=0.1, + tortuosity=2.0) +CASES = { # source, eta, (kxx, kxz, kzz) z up, viscous_q + "force_z": ("force", 0.0, (1e-9, 0.0, 1e-9), None), + "aniso_k": ("explosion", 1e-3, (2e-9, 5e-10, 5e-10), None), + "visc": ("explosion", 1e-3, (1e-9, 0.0, 1e-9), None), + "visc_atten": ("explosion", 1e-3, (1e-9, 0.0, 1e-9), (5.0, 30.0)), +} + + +def stf(t): + """SPECFEM2D's time_function_type 3 with t0 = 1.2 / f0 taken off ``t``.""" + a = (math.pi * F0) ** 2 + return FACTOR * np.exp(-a * t**2) / (2.0 * a) + + +def run_sem(case, device, dtype, n_order=4): + source, eta, (kxx, kxz, kzz), viscous_q = CASES[case] + dh = H_EL / n_order + nz, nx = int(round(LZ / dh)) + 1, int(round(LX / dh)) + 1 + t0 = 1.2 / F0 + t = np.arange(NT) * DT + coef = biot_from_rock(eta=eta, kappa=kxx, **ROCK) + eq = BiotVSSEM2DDeformed + if kxz != 0.0 or kxx != kzz: + eq = BiotAnisoPermSEM2DDeformed + coef = list(coef[:7]) + list(biot_drag_tensor(eta, kxx, -kxz, kzz)) + models = [torch.full((nz, nx), float(c), dtype=dtype) for c in coef] + src = np.array([[SRC[0] / dh, (LZ - SRC[1]) / dh]]) + rec = np.array([[[x / dh, (LZ - z) / dh] for x, z in RECS]]) + kw = dict(shape=(nz, nx), dh=dh, dt=DT, n_order=n_order, free_surface=True, f_dominant=F0, + dev=device, impl="eager", dtype=dtype, receiver_type=["vx", "vz", "p"]) + if viscous_q is not None: + kw["equation_options"] = {"viscous_q": viscous_q} + if source == "force": + # a downward force stf(t): fz takes the force itself, down positive; a load at + # step n is the force at n dt, a velocity sample n is at (n + 1/2) dt + calls = ((["fz"], stf(t - t0)),) + shifts = (0.5, 0.5, 1.0) + else: + # the stress and pressure sources take the rate; samples at (n + 1) dt (velocity) + # and (n + 3/2) dt (pressure) + rate = np.gradient(stf(t - t0), t) + calls = ((["sxx", "szz"], -rate), (["p"], rate)) + shifts = (1.0, 1.0, 1.5) + out = 0.0 + for names, w in calls: + prop = PropSEM(eq, source_type=names, **kw) + with torch.no_grad(): + out = out + prop(torch.from_numpy(w).to(dtype), src, rec, + models=models)[0].double().cpu().numpy() + out[:, :, 1] *= -1.0 # v_z down -> SPECFEM's up + return [(np.arange(NT) + s) * DT - t0 for s in shifts], out + + +def read_trace(path): + d = np.loadtxt(path) + return d[:, 0], d[:, 1] + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--runs", type=Path, required=True) + ap.add_argument("--out", type=Path, required=True) + ap.add_argument("--case", nargs="*", default=list(CASES)) + ap.add_argument("--device", default="cuda" if torch.cuda.is_available() else "cpu") + args = ap.parse_args() + args.out.mkdir(parents=True, exist_ok=True) + + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + for case in args.case: + t_sem, sem = run_sem(case, args.device, torch.float32) + np.savez(args.out / f"sem_{case}.npz", t=np.stack(t_sem), sem=sem) + fig, axes = plt.subplots(len(RECS), 2, figsize=(10, 2.3 * len(RECS)), sharex=True) + for r in range(len(RECS)): + for k, comp in enumerate(("BXX", "BXZ")): + ax = axes[r, k] + t_ref, ref = read_trace(args.runs / case / "OUTPUT_FILES" + / f"AA.S{r + 1:04d}.{comp}.semv") + ours = np.interp(t_ref, t_sem[k], sem[:, r, k]) + err = np.linalg.norm(ours - ref) / np.linalg.norm(ref) + corr = float(ours @ ref) / float(np.linalg.norm(ours) * np.linalg.norm(ref)) + amp = np.linalg.norm(ours) / np.linalg.norm(ref) + print(f"{case} rec{r + 1} {comp}: rel L2 {err:.4f} corr {corr:+.5f} amp {amp:.4f}", + flush=True) + ax.plot(t_ref, ref, "k", lw=1.4, label="SPECFEM2D") + ax.plot(t_ref, ours, "r--", lw=0.9, label="SEM") + ax.plot(t_ref, 5 * (ours - ref), color="0.6", lw=0.7, label="5x diff") + ax.set_title(f"({'abcdefgh'[2 * r + k]}) rec {r + 1} {comp}, rel. L2 {err:.4f}", + fontsize=8) + axes[0, 0].legend(frameon=False, fontsize=7) + for ax in axes[-1]: + ax.set_xlabel("time (s)") + fig.tight_layout() + fig.savefig(args.out / f"sem_vs_specfem2d_{case}.png", dpi=220) + plt.close(fig) + print("wrote", args.out / f"sem_vs_specfem2d_{case}.png", flush=True) + + +if __name__ == "__main__": + main() diff --git a/test/sem/test_biot_sem.py b/test/sem/test_biot_sem.py new file mode 100644 index 0000000..816cebd --- /dev/null +++ b/test/sem/test_biot_sem.py @@ -0,0 +1,914 @@ +"""BiotVSSEM2DDeformed (2-D velocity-stress Biot SEM): the undrained limit, the analytic line +source, the open-pore surface, gradient, sources and receivers, refusals. + +The references are the FD ``Biot`` tests' (test/test_biot.py). Time levels: a velocity sample +``n`` is ``v^{n+1/2}``, a stress / pressure sample ``sigma^{n+1}``, and a stress or pressure +source's integrated wavelet leads its rate by ``dt / 2`` -- so the records sit ``dt`` (velocity) +and ``3 dt / 2`` (pressure) after the step's time ``n dt``. +""" + +import functools +import math + +import numpy as np +import pytest +import torch + +from sweep.equations.biot import biot_fast_p_velocity, biot_from_biot_constants, biot_from_rock +from sweep.sem import (BiotAnisoPermSEM2DDeformed, BiotVSSEM2DDeformed, + ElasticVSSEM2DDeformed, PropSEM) +from test_biot import DAI1, _dwn_halfspace, _layered_rock + +NO = 4 +# Two elements of C-PML instead of PropSEM's ten, for tests whose C-PML is not under test: a +# gradient against finite differences (the discrete adjoint, whatever the layer), records that +# end before the first echo; on these small models the layer is most of the elements. +THIN_ABCN = 2 * NO + + +def _prop(shape, dh, dt, src, rec, *, eq=BiotVSSEM2DDeformed, dtype=torch.float32, f0=30.0, + impl="eager", dev="cpu", **kw): + return PropSEM(eq, shape=shape, dh=dh, dt=dt, n_order=NO, f_dominant=f0, dev=dev, + impl=impl, dtype=dtype, source_type=src, receiver_type=rec, **kw) + + +def _fields(coef, shape, dtype=torch.float32): + """Coefficient list -> arrays of ``shape`` (scalars broadcast).""" + return [torch.as_tensor(c, dtype=dtype).expand(shape).contiguous() for c in coef] + + +def _gauss(f0, t0): + return lambda tt: np.exp(-(math.pi * f0 * (tt - t0)) ** 2) + + +def _rate(s_time, t, dtype=torch.float32): + """The wavelet of a stress / pressure source whose history is ``s_time``: PropSEM + integrates it (the moment rate), as PropTorch's stress injection does.""" + return torch.from_numpy(np.gradient(s_time(t), t)).to(dtype) + + +def _ricker(nt, dt, f0, dtype=torch.float32): + a = (math.pi * f0 * (np.arange(nt) * dt - 1.2 / f0)) ** 2 + return torch.from_numpy((1 - 2 * a) * np.exp(-a)).to(dtype) + + +# ------------------------------------------------------------------------------------------ # +# limits and analytic solutions +# ------------------------------------------------------------------------------------------ # +@pytest.mark.parametrize("hill", [False, True]) +@pytest.mark.parametrize("source_type, elastic_source", [ + (["sxx", "szz"], ["sxx", "szz"]), (["vz"], ["vz"]), (["fz"], ["vz"])]) +def test_no_relative_flow_limit_is_the_undrained_elastic(source_type, elastic_source, hill): + # b -> infinity freezes q: v and the total stress obey ElasticVSSEM2DDeformed in the + # undrained medium, open-pore surface and all (no fluid crosses it), flat or under a + # topography; a force on both the solid and the fluid (fz) is then a force on the + # solid. Two layers, the models given on the GLL nodes so that both equations see + # the same moduli. Not to round-off: the alpha^2 M part of the undrained modulus passes + # through the continuous p, a weighted mean of div v over the elements at a node, where + # the elastic equation keeps every element's own -- 4e-3 for the explosion on 40 m + # elements, 1e-4 on 10 m ones. + shape, dh, dt, nt = (41, 49), 10.0, 5e-4, 500 + f64 = torch.float64 + src = np.array([[24, 18]]) + rec = np.stack([np.arange(4, 46, 3), np.full(14, 30)], 1)[None] + w = _ricker(nt, dt, 12.0, f64) + rows = 2.0 + 4.0 * np.exp(-((np.arange(shape[1]) - 24.0) / 8.0) ** 2) if hill else None + kw = dict(abcn=0, free_surface=True, topography=rows, dtype=f64) # no C-PML for both + biot = _prop(shape, dh, dt, source_type, ["vx", "vz"], **kw) + deep = (biot.mesh.z_e.mean((1, 2), keepdim=True) > 200.0).to(f64).expand(biot.mesh.z_e.shape) + rock = dict(k_s=36e9, rho_s=2650.0, k_f=2.25e9, rho_f=1000.0, eta=1e-3, k_d=8e9 + 6e9 * deep, + mu=6e9 + 5e9 * deep, phi=0.25 - 0.1 * deep, kappa=1e-13, tortuosity=2.0) + coef = [torch.as_tensor(c, dtype=f64).expand(deep.shape).contiguous() + for c in biot_from_rock(**rock)] + coef[7] = torch.full(deep.shape, 1e16, dtype=f64) + h, mu, rho = coef[0], coef[1], coef[4] + db = biot(w, src, rec, models=coef) + de = _prop(shape, dh, dt, elastic_source, ["vx", "vz"], eq=ElasticVSSEM2DDeformed, **kw)( + w, src, rec, models=[torch.sqrt(h / rho), torch.sqrt(mu / rho), rho]) + err = float((db - de).abs().max() / de.abs().max()) + assert err < 1e-2, err + + +def analytic_line_source(models, t, x, z, s_time, force, b_of_w=None): + """``(vz, p)`` of a line source in a homogeneous Biot medium, the source as in + :func:`_analytic_line_source_vz` with forces ``force = (total, relative) * grad Phi``, + ``Phi = S(t) delta``, exp(+i w t). With ``u = grad Psi_u``, ``w = grad Psi_w`` and + ``p = -C div u - M div w``, mode ``k`` adds ``E_u a k S (z / r) (-w / 4) H1(k r)`` to vz + and ``k^2 (C E_u + M E_w) a S (-i/4) H0(k r)`` to p. ``b_of_w(w)`` replaces the drag + (complex, frequency dependent); default the constant ``b``.""" + from scipy.special import hankel2e + + h, mu, c, biot_m, rho, rho_f, m, b = (float(v) for v in models) + dt, nt = t[1] - t[0], len(t) + n = 16 * nt + s_hat = np.fft.rfft(s_time(np.arange(n) * dt)) + w = 2 * math.pi * np.fft.rfftfreq(n, dt) + r = math.hypot(x, z) + K = np.array([[h, c], [c, biot_m]], dtype=complex) + F = np.asarray(force, dtype=complex) + v_hat, p_hat = np.zeros_like(s_hat), np.zeros_like(s_hat) + for i in np.nonzero(np.abs(s_hat) > 1e-14 * np.abs(s_hat).max())[0]: + if w[i] == 0: + continue + bw = b if b_of_w is None else b_of_w(w[i]) + s2, E = np.linalg.eig(np.linalg.solve(K, np.array([[rho, rho_f], [rho_f, m - 1j * bw / w[i]]]))) + s = np.sqrt(s2) + s = np.where(s.real < 0, -s, s) + a = np.linalg.solve(K @ E, F) + k = w[i] * s + v_hat[i] = -(w[i] / 4) * (z / r) * s_hat[i] * np.sum( + E[0] * a * k * hankel2e(1, k * r) * np.exp(-1j * k * r)) + p_hat[i] = -0.25j * s_hat[i] * np.sum( + k**2 * (c * E[0] + biot_m * E[1]) * a * hankel2e(0, k * r) * np.exp(-1j * k * r)) + return np.fft.irfft(v_hat, n)[:nt], np.fft.irfft(p_hat, n)[:nt] + + +def _split_source(shape, dh, dt, nt, models, src, rec, rec_types, *, qpart, phi, f0, t0, + dev="cpu", **kw): + """The record of Dai's split source: total stress plus ``qpart / phi`` in the relative + flow, both of history ``exp(-(pi f0 (t - t0))^2)``, as two PropSEM calls.""" + t = np.arange(nt) * dt + w = _rate(_gauss(f0, t0), t) + coef = _fields(models, shape) + num = 0.0 + for fields, scale in ((["sxx", "szz"], 1.0), (["p"], -qpart / phi)): + prop = _prop(shape, dh, dt, fields, rec_types, dev=dev, **kw) + with torch.no_grad(): + num = num + prop(w * scale, src, rec, models=coef)[0].double().cpu().numpy() + return t, num + + +@functools.lru_cache(maxsize=None) +def line_source_case(b_dai, dev="cpu", viscous_q=None): + """``(t, sem (nt, 2), analytic (nt, 2))``: solid vz and pore pressure of the split source + in Dai's medium, 60 m off in x and z. The model keeps the C-PML's echoes out of the record + (they reach a few percent of the direct waves; six elements of it instead of ten add at + most 6e-4 to the errors). ``viscous_q`` turns on the viscous attenuation, the analytic drag + then ``b (1 + i w tau_e) / (1 + i w tau_s)``. Cached: the plain b = 7e4 case serves three + tests.""" + from sweep.sem.equations.biot_vs_deformed import viscous_relaxation_times + + dh, dt, nt, f0, t0, qpart, phi = 4.0, 2e-4, 800, 30.0, 0.04, 0.1, 0.3 + models = biot_from_biot_constants(**DAI1, b=b_dai, phi=phi) + kw = {} if viscous_q is None else {"equation_options": {"viscous_q": viscous_q}} + t, num = _split_source((121, 121), dh, dt, nt, models, np.array([[45, 45]]), + np.array([[[60, 60]]]), ["vz", "p"], qpart=qpart, phi=phi, f0=f0, + t0=t0, dev=dev, abcn=6 * NO, **kw) + b_of_w = None + if viscous_q is not None: + tau_e, tau_s = viscous_relaxation_times(*viscous_q) + b = float(models[7]) + b_of_w = lambda w: b * (1 + 1j * w * tau_e) / (1 + 1j * w * tau_s) + s_time = _gauss(f0, t0) + vz, _ = analytic_line_source(models, t, 60.0, 60.0, lambda tt: s_time(tt + dt), + (1.0, qpart / phi), b_of_w) + _, p = analytic_line_source(models, t, 60.0, 60.0, lambda tt: s_time(tt + 1.5 * dt), + (1.0, qpart / phi), b_of_w) + return t, num[:, 0], np.stack([vz, p], 1) + + +def open_pore_case(b_dai, dev="cpu"): + """``(t, sem (nt, 2), wavenumber integration (nt, 2))``: solid vx, vz on the open-pore + surface 60 m from the split source 20 m below it, the model just large enough that the + C-PML's echoes reach no more than the end of the record.""" + dh, dt, nt, f0, t0, qpart, phi = 2.0, 1e-4, 1200, 30.0, 0.04, 0.1, 0.3 + models = biot_from_biot_constants(**DAI1, b=b_dai, phi=phi) + xs, zs, xr = 65, 10, 95 + t, num = _split_source((85, 161), dh, dt, nt, models, np.array([[xs, zs]]), + np.array([[[xr, 0]]]), ["vx", "vz"], qpart=qpart, phi=phi, f0=f0, + t0=t0, dev=dev, free_surface=True) + s_time = _gauss(f0, t0) + refs = _dwn_halfspace(models, (1.0, qpart / phi), zs * dh, [((xr - xs) * dh, 0.0)], t, + lambda tt: s_time(tt + dt)) + return t, num[:, 0], np.stack([r[:, 0] for r in refs], 1) + + +def _rel(a, b): + return float(np.linalg.norm(a - b) / np.linalg.norm(b)) + + +@pytest.mark.parametrize("b_dai, tol", [(0.0, 0.01), (7e4, 0.005), (9e8, 0.005)]) +def test_viscous_line_source_matches_analytic(b_dai, tol): + # slow wave propagating, critically damped, diffusive + _, num, ana = line_source_case(b_dai) + for k, comp in enumerate(("vz", "p")): + assert _rel(num[:, k], ana[:, k]) < tol, (comp, _rel(num[:, k], ana[:, k])) + + +# In the Darcy limit the surface's drainage layer sqrt(2 D / w) is ~7 cm, far below the GLL +# spacing, and the nodes there are fully drained: 0.9% on 8 m elements 60 m from the source, +# 1.2% at 100 m (FD Biot there: 1.3% on 2 m cells). With the layer resolved (b = 7e5) the error +# falls to 1e-4 under refinement. +@pytest.mark.parametrize("b_dai, tol", [(0.0, 0.005), (7e5, 0.005), (9e8, 0.02)]) +def test_open_pore_surface_matches_wavenumber_integration(b_dai, tol): + # sigma . n = 0, p = 0 on the surface: slow wave propagating, diffusive, Darcy limit + _, num, ref = open_pore_case(b_dai) + for k in range(2): + assert _rel(num[:, k], ref[:, k]) < tol, (("vx", "vz")[k], _rel(num[:, k], ref[:, k])) + + +# the pore fluid on its own: C = 0 and rho_f = 0 leave an acoustic medium (modulus M, +# density m), excited by a pore-pressure source +FLUID = [2e10, 6e9, 0.0, 4e9, 2300.0, 0.0, 2000.0, 0.0] + + +def pressure_source_case(shape, src, rec, free_surface, dh=2.0, dt=2e-4, nt=900, dev="cpu"): + """``(t, sem p (nt, nrec), analytic (nt, nrec))`` of a pore pressure ``S(t)`` injected + in :data:`FLUID`: the free-space field, minus its image in an open-pore surface.""" + from scipy.special import hankel2e + + f0, t0 = 30.0, 0.04 + t = np.arange(nt) * dt + s_time = _gauss(f0, t0) + prop = _prop(shape, dh, dt, ["p"], ["p"], free_surface=free_surface, dev=dev, + abcn=THIN_ABCN) + with torch.no_grad(): + num = prop(_rate(s_time, t), np.array([src]), np.array([rec]), + models=_fields(FLUID, shape))[0, :, :, 0].double().cpu().numpy() + n = 16 * nt + om = 2 * math.pi * np.fft.rfftfreq(n, dt)[1:] # k^2 H0(k r) -> 0 at w = 0 + s_hat = np.fft.rfft(s_time(np.arange(n) * dt))[1:] * np.exp(1.5j * om * dt) + k = om * math.sqrt(FLUID[6] / FLUID[3]) + ana = [] + for xr, zr in rec: + r1 = dh * math.hypot(xr - src[0], zr - src[1]) + green = hankel2e(0, k * r1) * np.exp(-1j * k * r1) + if free_surface: + r2 = dh * math.hypot(xr - src[0], zr + src[1]) + green = green - hankel2e(0, k * r2) * np.exp(-1j * k * r2) + ana.append(np.fft.irfft(np.concatenate([[0.0], 0.25j * k**2 * green * s_hat]), n)[:nt]) + return t, num, np.stack(ana, 1) + + +def test_decoupled_pore_fluid_sees_a_pressure_release_surface(): + # p = 0 on the open-pore surface: the free-space field minus the image source's, on + # receivers one element under the surface and deeper + _, num, ana = pressure_source_case((101, 257), [94, 15], + [[124, 4], [154, 4], [124, 30], [154, 30]], True) + for i in range(num.shape[1]): + assert _rel(num[:, i], ana[:, i]) < 0.005, (i, _rel(num[:, i], ana[:, i])) + + +def test_pressure_has_no_grid_line_modes(): + # a point source excites nothing that runs along the grid lines through it: a receiver + # on the source's row records the free-space pressure like one off it (with a continuous + # flux and an element-local p, the row's record was 14x the analytic at dh = 2 and grew + # like 1 / dh) + _, num, ana = pressure_source_case((161, 161), [65, 80], [[95, 80], [95, 95]], False) + for i in range(2): + assert _rel(num[:, i], ana[:, i]) < 0.005, (i, _rel(num[:, i], ana[:, i])) + + +def test_matches_the_fd_biot(): + # PropTorch's Biot and PropSEM on the same call: the FD stress injection is a moment of + # dh^2 per unit wavelet sum, the SEM one of dt, and FD's vz sits half a cell down + from sweep.equations import Biot + from sweep.propagator.options import EagerOptions + from sweep.propagator.torch import PropTorch + + dh, dt, nt, f0 = 2.0, 2e-4, 900, 25.0 + shape = (81, 101) + models = biot_from_biot_constants(**DAI1, b=7e4, phi=0.3) + w = _ricker(nt, dt, f0) + src, rec = np.array([[30, 30]]), np.array([[[70, 50]]]) + # uncompiled: no other test steps this grid, which a cold compile makes 8 - 20 s, not 1.3 s + fd = PropTorch(Biot(device="cpu"), shape=shape, dev="cpu", dh=dh, dt=dt, source_type=["sxx", "szz"], + receiver_type=["vz"], abcn=40, backend="torch", impl="eager", + eager_options=EagerOptions(use_compile=False)) + sem = _prop(shape, dh, dt, ["sxx", "szz"], ["vz"]) + with torch.no_grad(): + a = fd(w, src, rec, models=_fields(models, shape))[0, :, 0, 0].double().numpy() + b = sem(w, src, rec + np.array([0, 0.5]), models=_fields(models, shape))[0, :, 0, 0].double().numpy() + b = b * dh**2 / dt + corr = float(a @ b) / float(np.linalg.norm(a) * np.linalg.norm(b)) + assert corr > 0.995, corr + assert np.linalg.norm(b) / np.linalg.norm(a) == pytest.approx(1.0, abs=0.03) + + +# ------------------------------------------------------------------------------------------ # +# gradient +# ------------------------------------------------------------------------------------------ # +# The finite-difference checks: a central difference of a 1e-5 relative step against the +# adjoint, to 1e-6 and no absolute floor (|J' dm| is 1e-21 - 1e-16 here: pytest.approx's default +# abs = 1e-12 would pass any gradient). Measured 1e-10 - 1e-7, the largest the step's own +# truncation (M in the water); a 1e-3 step was 1e-3 off, edge cells 2e-4, a moving C-PML 2e-3. +FD_STEP = 1e-5 + + +def _fd_matches(J, coef, which, grad, dm): + pred = float((grad * dm).sum()) + with torch.no_grad(): + fd = (float(J([c + dm if k == which else c for k, c in enumerate(coef)])) + - float(J([c - dm if k == which else c for k, c in enumerate(coef)]))) / 2.0 + assert fd == pytest.approx(pred, rel=1e-6, abs=0), (which, fd, pred) + + +def _interior_bump(shape): + """A Gaussian bump on the model grid, zero on its outermost cells: the C-PML copies those + outward and the gradient is the physical box's, so theirs leaves the copies out.""" + zz, xx = torch.meshgrid(*(torch.arange(n, dtype=torch.float64) for n in shape), indexing="ij") + bump = torch.exp(-((zz - 18.0) ** 2 + (xx - 16.0) ** 2) / 30.0) + bump[[0, -1]] = 0.0 + bump[:, [0, -1]] = 0.0 + return bump + + +@functools.lru_cache(maxsize=None) +def _gradient_case(free_surface=False, variant=None): + """``(coef, J, grads)``: the misfit of a layered rock and its gradient to every + coefficient, one backward for every coefficient checked. The C-PML is laid out for the + model's fast P speed, a float the gradient does not follow: held, or a step in M would + move the layer too.""" + shape, dh, dt, nt = (33, 41), 10.0, 5e-4, 360 + f64 = torch.float64 + rock = _layered_rock(shape) + rock["kappa"] = 1e-10 # gamma dt ~ 1: the drag's exponential matters + coef = _fields(biot_from_rock(**rock), shape, f64) + eq, kw = BiotVSSEM2DDeformed, {} + if variant == "tensor": + eq = BiotAnisoPermSEM2DDeformed + coef = coef[:7] + [coef[7], 0.3 * coef[7], 2.0 * coef[7]] + elif variant == "viscous": + kw["equation_options"] = {"viscous_q": (5.0, 12.0)} + src = np.array([[20, 10]]) + rec = np.stack([np.arange(6, 36, 4), np.full(8, 24)], 1)[None] + w = _ricker(nt, dt, 12.0, f64) + prop = _prop(shape, dh, dt, ["sxx", "szz"], ["vx", "vz", "p"], eq=eq, dtype=f64, f0=12.0, + free_surface=free_surface, abcn=THIN_ABCN, **kw) + c_max = prop._pml_speed([prop._interp_to_gll(c) for c in coef]) + prop._pml_speed = lambda models: c_max + obs = prop(w, src, rec, models=coef).detach() + + def J(models): + return 0.5 * (prop(w, src, rec, models=models) - 1.05 * obs).square().sum() + + m0 = [c.clone().requires_grad_() for c in coef] + return coef, J, torch.autograd.grad(J(m0), m0) + + +@pytest.mark.parametrize("free_surface", [False, True]) +@pytest.mark.parametrize("which", [1, 3, 7]) # mu, M, b +def test_gradient_matches_finite_difference(free_surface, which): + coef, J, grads = _gradient_case(free_surface) + dm = FD_STEP * coef[which] * _interior_bump(coef[which].shape) + _fd_matches(J, coef, which, grads[which], dm) + + +@pytest.mark.parametrize("variant", ["tensor", "viscous"]) +def test_gradient_through_the_drag_tensor_and_the_viscous_memory(variant): + # the per-node matrix exponential and the memory variable carry the gradient too + coef, J, grads = _gradient_case(variant=variant) + which = 8 if variant == "tensor" else 7 # b_xz, b + _fd_matches(J, coef, which, grads[which], FD_STEP * coef[7] * _interior_bump(coef[7].shape)) + + +# ------------------------------------------------------------------------------------------ # +# sources, receivers, wavefields +# ------------------------------------------------------------------------------------------ # +def _small(src=None, rec=None, **kw): + shape = (21, 25) + prop = _prop(shape, 10.0, 5e-4, src, rec, f0=12.0, **kw) + coef = _fields(biot_from_biot_constants(**DAI1, b=3e3, phi=0.3), shape) + return prop, coef + + +def test_source_names_are_fd_biots_and_specfems(): + # payloads over (v_x, v_z, p, q_x, q_z); the stress names are a moment tensor + p, _ = _small() + assert p.source_type == ["sxx", "szz"] and p._src_moment == (1.0, 1.0, 0.0) + assert p._src_mode == "load" + for names, payload, mode in ((["p"], (0, 0, 1, 0, 0), "velocity"), + (["vz", "p", "qx"], (0, 1, 1, 1, 0), "velocity"), + (["fz"], (0, 1, 0, 0, 1), "load"), + (["fx", "p"], (1, 0, 1, 1, 0), "load")): + p, _ = _small(names) + assert p._src_moment is None and p._source_vector() == payload, names + assert p._src_mode == mode, names + + +@pytest.mark.parametrize("names, match", [ + (["vz", "sxx"], "mixes a stress source"), + (["sxx", "szz", "p"], "mixes a stress source"), + (["vz", "v_z"], "both drive v_z"), + (["fz", "qz"], "both drive q_z"), + (["fx", "vz"], "mixes an injected velocity with a force"), + (["u_z"], "Unknown source_type entry 'u_z'"), + (["sxx", "sxx"], "duplicate"), +]) +def test_bad_source_type_is_refused(names, match): + with pytest.raises(ValueError, match=match): + _small(names) + + +def _one_step(names, field): + """One step from rest with a unit load on the node mid-mesh in source field ``field``.""" + p, coef = _small(names) + eq = p._build_eq(coef) + state = eq.init_state(p.dt, "cpu", torch.float32) + load = torch.zeros(eq.mesh.n_global).index_fill(0, torch.tensor([eq.mesh.n_global // 2]), 1.0) + src = [torch.zeros_like(load) for _ in range(5)] + src[field] = load + return p, eq, load, eq.step(state, src, p.dt) + + +@pytest.mark.parametrize("field", [0, 1]) +def test_a_velocity_source_moves_only_the_solid(field): + # vx / vz inject dt src / M_rho into v, the flux does not move, whatever the drag + p, eq, load, out = _one_step(["vz"], field) + assert eq.source_mode == "velocity" + torch.testing.assert_close(out[field], p.dt * eq._inv_M * load, rtol=0, atol=0) + torch.testing.assert_close(out[1 - field], torch.zeros_like(load), rtol=0, atol=0) + for k in (3, 4): # q_x, q_z + assert not out[k].any() + + +@pytest.mark.parametrize("field", [3, 4]) +def test_a_flux_source_moves_only_the_fluid(field): + # qx / qz: every element copy of the flux at the node jumps by dt src / A(w m), v not + p, eq, load, out = _one_step(["qz"], field) + want = (p.dt * eq._inv_Mm * load)[eq.mesh.ibool] + torch.testing.assert_close(out[field], want, rtol=0, atol=0) + assert not out[7 - field].any() and not out[0].any() and not out[1].any() + + +def test_a_force_loads_both_equations(): + # fz is SPECFEM2D's poroelastic point force: the same load on the total momentum and on + # the relative flow, i.e. on the solid and the fluid in proportion 1 - phi : phi + from sweep.sem.prop_base import _Acq + p, coef = _small(["fz"]) + eq = p._build_eq(coef) + assert eq.source_mode == "load" + w = _ricker(40, 5e-4, 12.0) + p._set_nt(w) + acq = _Acq(p, w, np.array([[12, 10]]), np.array([[[6, 6]]]), None) + (prop,) = p._build_props(eq, None, acq, "warn") + pats = prop._build_source_patterns(eq.mesh.n_global, "cpu", torch.float32) + assert pats[1].any() + torch.testing.assert_close(pats[4], pats[1], rtol=0, atol=0) + for k in (0, 2, 3): + assert not pats[k].any() + + +def test_a_stress_source_is_a_load(): + p, coef = _small(["sxz"]) + assert p._build_eq(coef).source_mode == "load" + + +def test_default_receivers_are_the_solid_velocity(): + p, coef = _small() + assert p.receiver_type == ["v_x", "v_z"] + out = p(_ricker(80, 5e-4, 12.0), np.array([[12, 10]]), np.array([[[6, 6], [18, 14]]]), models=coef) + assert out.shape == (1, 80, 2, 2) + + +def test_pressure_flux_and_stress_receivers_on_shots_and_snapshots(): + # the flux and the stress are element-local: they are gathered from their own layout, + # in the vmapped multi-shot batch as in the one-shot loop, and snapshotted on the grid + # with the continuous p and v. The batch runs batched kernels, so it agrees to rounding + # (torch 2.14 on CPU: 2e-6 of each field's peak), not bitwise. + p, coef = _small(["sxx", "szz"], ["p", "sxx", "qz", "vz"]) + w = _ricker(120, 5e-4, 12.0) + src = np.array([[12, 10], [8, 6]]) + rec = np.array([[[6, 7], [18, 14], [10, 3]]]) + with torch.no_grad(): + both = p(w, src, rec, models=coef) + one = torch.cat([p(w, src[i:i + 1], rec, models=coef) for i in range(2)]) + peak = one.abs().amax((0, 1, 2)) + assert both.shape == (2, 120, 3, 4) and bool((peak > 0).all()) + torch.testing.assert_close(both / peak, one / peak, rtol=0, atol=1e-5) + with torch.no_grad(): + seis, snaps = p(w, src[:1], rec, models=coef, return_wavefield=True, + snapshot_times=[40, 100]) + assert snaps.shape == (2, 4, 1, 1) + p._grid_shape() + torch.testing.assert_close(seis, one[:1], rtol=0, atol=0) + + +def test_boundary_saving_and_cuda_are_refused(): + from sweep.propagator.options import BoundarySaving + with pytest.raises(NotImplementedError, match="no boundary saving"): + _small(memory=BoundarySaving()) + with pytest.raises(NotImplementedError, match="no CUDA kernel"): + _small(impl="c") + + +def test_source_vector_takes_one_entry_per_source_field(): + p, coef = _small(["vz"]) + rec = np.array([[[6, 6]]]) + with pytest.raises(ValueError, match="one entry per source field"): + p(_ricker(40, 5e-4, 12.0), np.array([[12, 10]]), rec, models=coef, source_vector=(0.0, 1.0)) + a = p(_ricker(40, 5e-4, 12.0), np.array([[12, 10]]), rec, models=coef, + source_vector=(0.0, 1.0, 0.0, 0.0, 0.0)) + b = p(_ricker(40, 5e-4, 12.0), np.array([[12, 10]]), rec, models=coef) + torch.testing.assert_close(a, b, rtol=0, atol=0) + + +def test_the_pml_is_designed_for_the_fast_p_wave(): + # not for the first model, which is a modulus here + p, coef = _small() + assert p._pml_speed(coef) == pytest.approx(biot_fast_p_velocity(coef)) + assert 3000.0 < p._pml_speed(coef) < 3500.0 + + +# ------------------------------------------------------------------------------------------ # +# anisotropic permeability, viscous attenuation, parameter families, options +# ------------------------------------------------------------------------------------------ # +def _explosion(eq, models, rec, shape=(33, 33), src=(16, 16), nt=300, **kw): + f64 = torch.float64 + prop = _prop(shape, 4.0, 2e-4, ["sxx", "szz"], ["vx", "vz", "p"], eq=eq, dtype=f64, + abcn=THIN_ABCN, **kw) + with torch.no_grad(): + return prop(_ricker(nt, 2e-4, 30.0, f64), np.array([src]), np.array([rec]), + models=_fields(models, shape, f64))[0].numpy() + + +def test_an_isotropic_drag_tensor_is_the_scalar_drag(): + models = biot_from_biot_constants(**DAI1, b=7e4, phi=0.3) + b = models[7] + rec = [[24, 8], [8, 26]] + iso = _explosion(BiotVSSEM2DDeformed, models, rec) + tensor = _explosion(BiotAnisoPermSEM2DDeformed, list(models[:7]) + [b, 0.0, b], rec) + assert _rel(tensor, iso) < 1e-12, _rel(tensor, iso) + + +def test_the_drag_tensor_turns_with_the_medium(): + # mirroring about x = z swaps b_xx and b_zz, keeps b_xz, and swaps x and z of every + # record: the off-diagonal coupling of the flux copies included + from sweep.equations.biot_params import biot_drag_tensor + + models = list(biot_from_biot_constants(**DAI1, b=0.0, phi=0.3)[:7]) + b_xx, b_xz, b_zz = biot_drag_tensor(1e-3, 2e-9, -5e-10, 5e-10) + a = _explosion(BiotAnisoPermSEM2DDeformed, models + [b_xx, b_xz, b_zz], [[24, 10]]) + b = _explosion(BiotAnisoPermSEM2DDeformed, models + [b_zz, b_xz, b_xx], [[10, 24]]) + assert _rel(b[:, 0, 1], a[:, 0, 0]) < 1e-10 and _rel(b[:, 0, 0], a[:, 0, 1]) < 1e-10 + assert _rel(b[:, 0, 2], a[:, 0, 2]) < 1e-10 + # ... and the tensor is felt: the isotropic mean drag gives another record + mean = 0.5 * (b_xx + b_zz) + c = _explosion(BiotAnisoPermSEM2DDeformed, models + [mean, 0.0, mean], [[24, 10]]) + assert _rel(c[:, 0, 2], a[:, 0, 2]) > 1e-2 + + +def test_a_drag_tensor_must_be_positive_semidefinite(): + models = list(biot_from_biot_constants(**DAI1, b=0.0, phi=0.3)[:7]) + with pytest.raises(ValueError, match="positive semi-definite"): + _explosion(BiotAnisoPermSEM2DDeformed, models + [1e5, 2e5, 1e5], [[24, 10]]) + + +@pytest.mark.parametrize("q0", [1.0, 5.0]) +def test_viscous_attenuation_matches_the_analytic_relaxation(q0): + # SPECFEM2D's ATTENUATION_PORO_FLUID_PART: the drag a standard linear solid, Q0 at 30 Hz, + # against the analytic line source with b (1 + i w tau_e) / (1 + i w tau_s). With the + # slow wave critically damped (b = 7e4) it moves the pore pressure by 1.3% (Q0 = 1) and + # 0.3% (Q0 = 5) and the solid velocity by less than 0.1% -- so its own part, num - plain + # against ana - ana_plain, is checked on the pore pressure: 2.5% and 1.1% off + _, num, ana = line_source_case(7e4, viscous_q=(q0, 30.0)) + _, plain, ana0 = line_source_case(7e4) + for k, comp in enumerate(("vz", "p")): + assert _rel(num[:, k], ana[:, k]) < 0.005, (comp, _rel(num[:, k], ana[:, k])) + assert _rel(ana[:, 1], ana0[:, 1]) > 0.002 + err = _rel(num[:, 1] - plain[:, 1], ana[:, 1] - ana0[:, 1]) + assert err < 0.05, err + + +def test_equation_options_are_checked(): + p, coef = _small(equation_options={"viscous_q": (5.0, 30.0)}) + assert p._build_eq(coef).viscous_q == (5.0, 30.0) + with pytest.raises(TypeError, match="not options of BiotVSSEM2DDeformed.build"): + _small(equation_options={"viscous": (5.0, 30.0)}) + p, coef = _small(equation_options={"viscous_q": (0.0, 30.0)}) + with pytest.raises(ValueError, match="Q0 > 0"): + p._build_eq(coef) + + +@pytest.mark.parametrize("family", ["primary", "density", "wavespeed"]) +@pytest.mark.parametrize("tensor", [False, True]) +def test_parameter_families_round_trip_and_chain_rule(family, tensor): + from sweep.equations.biot_params import (BIOT_FAMILIES, biot_family_kernels, + biot_from_family, biot_to_family) + rocks = [biot_from_biot_constants(**DAI1, b=b, phi=phi) for b, phi in ((0.0, 0.3), (7e4, 0.25))] + coef = [torch.tensor([float(r[i]) for r in rocks], dtype=torch.float64) for i in range(8)] + if tensor: + coef = coef[:7] + [coef[7], 0.3 * coef[7], 2.0 * coef[7]] + y = biot_to_family(family, coef) + for got, want in zip(biot_from_family(family, y), coef): + assert float(((got - want).abs() / want.abs().clamp_min(1e-300)).max()) < 1e-12 + # a nonlinear misfit of the coefficients; its gradient taken to the family against + # central differences in the family's own coordinates, one medium at a time + torch.manual_seed(0) + weights = [torch.rand(2, dtype=torch.float64) / (c.abs() + 1.0) for c in coef] + + def misfit(cs): + return sum(((w * c) ** 2).sum() for w, c in zip(weights, cs)) + + cf = [c.clone().requires_grad_() for c in coef] + grads = torch.autograd.grad(misfit(cf), cf) + kern = biot_family_kernels(family, coef, grads, relative=False) + names = list(BIOT_FAMILIES[family]) + if tensor: + i = names.index("b") + names[i:i + 1] = ["b_xx", "b_xz", "b_zz"] + assert list(kern) == names + for j, name in enumerate(names): + for q in range(2): + h = 1e-6 * float(y[j][q].abs().clamp_min(1e-30)) + step = torch.zeros(2, dtype=torch.float64) + step[q] = h + up = [yy + (step if k == j else 0.0) for k, yy in enumerate(y)] + dn = [yy - (step if k == j else 0.0) for k, yy in enumerate(y)] + fd = float(misfit(biot_from_family(family, up)) + - misfit(biot_from_family(family, dn))) / (2 * h) + assert float(kern[name][q]) == pytest.approx(fd, rel=1e-5, abs=1e-300), (name, q) + rel = biot_family_kernels(family, coef, grads) + for j, name in enumerate(names): + assert torch.allclose(rel[name], y[j] * kern[name]) + + + +# ------------------------------------------------------------------------------------------ # +# fluid and dry media, their interfaces +# ------------------------------------------------------------------------------------------ # +WATER = (2.25e9, 1000.0) # k_f, rho_f +SHALE = (3000.0, 1600.0, 2200.0) # vp, vs, rho + + +def _layered(mesh, tops, media, dtype=torch.float64): + """GLL-native coefficients: an element takes ``media[k]`` when its centre lies between + the depths ``tops[k]`` and ``tops[k + 1]`` -- each element wholly one medium.""" + zc = mesh.z_e.mean((1, 2)) + k = torch.bucketize(zc, torch.tensor(tops[1:], dtype=zc.dtype, device=zc.device)) + E, n = mesh.n_elements, mesh.ngll + return [torch.stack([torch.as_tensor(m[i], dtype=dtype) for m in media]).to(zc.device)[k] + .view(E, 1, 1).expand(E, n, n).contiguous() for i in range(len(media[0]))] + + +def _three_media(mesh, b=0.0, tops=(0.0, 60.0, 140.0)): + """Water over Dai's porous rock (drag ``b``) over a dry shale.""" + from sweep.equations.biot_params import biot_elastic, biot_fluid + poro = biot_from_biot_constants(**DAI1, b=b, phi=0.3) + return _layered(mesh, list(tops), [biot_fluid(*WATER), poro, biot_elastic(*SHALE)]) + + +def _node(mesh, ez, ex, a=2, b=2): + """The global node at GLL ``(a, b)`` of element row ``ez``, column ``ex``.""" + return int(mesh.ibool[ez * mesh.nex + ex, a, b]) + + +def test_a_dry_solid_is_the_elastic_sem(): + # no pore fluid anywhere (M = 0): v and the stress obey ElasticVSSEM2DDeformed and the + # pressure and the flux stay zero + from sweep.equations.biot_params import biot_elastic + shape, dh, dt, nt = (33, 41), 10.0, 5e-4, 300 + f64 = torch.float64 + kw = dict(abcn=0, free_surface=True, dtype=f64) + rocks = [(2500.0, 1400.0, 2200.0), SHALE] + biot = _prop(shape, dh, dt, ["sxx", "szz"], ["vx", "vz", "p", "qz"], **kw) + elas = _prop(shape, dh, dt, ["sxx", "szz"], ["vx", "vz"], eq=ElasticVSSEM2DDeformed, **kw) + src, rec = np.array([[20, 10]]), np.stack([np.arange(4, 37, 4), np.full(9, 22)], 1)[None] + w = _ricker(nt, dt, 12.0, f64) + with torch.no_grad(): + db = biot(w, src, rec, models=_layered(biot.mesh, [0.0, 150.0], + [biot_elastic(*r) for r in rocks])) + de = elas(w, src, rec, models=_layered(elas.mesh, [0.0, 150.0], rocks)) + assert not db[..., 2:].any() + err = float((db[..., :2] - de).abs().max() / de.abs().max()) + assert err < 1e-11, err + + +def test_water_over_a_dry_solid_is_the_fluid_solid_sem(): + # AcousticElasticSEM2D (potential and displacement, the same FluidSolidInterface) and the + # Biot SEM with a fluid over a dry solid are one discrete system: p = -chi'' and v = u', + # a force in the solid, the free surface on the water, no C-PML + from sweep.equations.biot_params import biot_elastic, biot_fluid + from sweep.sem.equations import AcousticElasticSEM2D + f64 = torch.float64 + shape, dh, dt, nt = (41, 61), 5.0, 5e-4, 400 + prop = _prop(shape, dh, dt, ["fz"], ["vz"], abcn=0, free_surface=True, dtype=f64) + mesh = prop.mesh + tops = [0.0, 60.0] + eq_b = prop._build_eq(_layered(mesh, tops, [biot_fluid(*WATER), biot_elastic(*SHALE)])) + vp, vs, rho = _layered(mesh, tops, [(1500.0, 0.0, 1000.0), SHALE]) + eq_a = AcousticElasticSEM2D.build(mesh, vp, vs, rho, dirichlet=("top",)) + load = torch.zeros(mesh.n_global, dtype=f64) + load[_node(mesh, 5, 7, 1, 3)] = 1.0 + w = _ricker(nt, dt, 20.0, f64) + rec_v = [_node(mesh, 3, 4, 0, 0), _node(mesh, 6, 11)] # the seafloor, in the rock + rec_p = [_node(mesh, 1, 9), _node(mesh, 2, 3, 4, 1)] # the water, the seafloor + sb, sa = eq_b.init_state(dt, "cpu", f64), eq_a.init_state(dt, "cpu", f64) + zero = torch.zeros_like(load) + v, p, u, chi = [], [], [], [] + for n in range(nt): + sb = eq_b.step(sb, [zero, w[n] * load, zero, zero, zero], dt) + sa = eq_a.step(sa, (zero, zero, w[n] * load), dt) + v.append(sb[1][rec_v]); p.append(sb[2][rec_p]) + u.append(sa[3][rec_v]); chi.append(sa[0][rec_p]) + v, p, u, chi = (torch.stack(x) for x in (v, p, u, chi)) + u = torch.cat([torch.zeros_like(u[:1]), u]) # u^0 ... u^nt + chi = torch.cat([torch.zeros_like(chi[:1]), chi]) + v_a = (u[1:] - u[:-1]) / dt # v^{n+1/2} + p_a = -(chi[2:] - 2 * chi[1:-1] + chi[:-2]) / dt**2 # p^{n+1} = -chi''^{n+1} + assert float(v.abs().max()) > 0 and float(p.abs().max()) > 0 + torch.testing.assert_close(v, v_a, rtol=0, atol=1e-9 * float(v_a.abs().max())) + torch.testing.assert_close(p[:-1], p_a, rtol=0, atol=1e-9 * float(p_a.abs().max())) + + +def _green(eq, dt, nt, field, src_node, rec_field, rec_node, w): + """``rec_field`` at ``rec_node`` after each of ``nt`` steps from rest, a unit load + ``w[n]`` at ``src_node``: on the momentum (``field`` 0, 1), or a fluid volume on the + pressure equation (``field`` 2): ``p`` and, as the flow does, the total stress by + ``-alpha dp`` -- not the source ``p``, which leaves the total stress and so carries a + stress source in a porous rock.""" + state = eq.init_state(dt, "cpu", torch.float64) + zero = torch.zeros(eq.mesh.n_global, dtype=torch.float64) + ibool = eq.mesh.ibool + out = [] + for n in range(nt): + src = [zero] * 5 + if field == 2: + dp = zero.clone() + dp[src_node] = dt * w[n] * eq._inv_Mp[src_node] + state = list(state) + state[2] = state[2] + dp + for k in (5, 6): # sigma_xx, sigma_zz + state[k] = state[k] - eq._alpha * dp[ibool] + else: + src[field] = zero.index_fill(0, torch.tensor([src_node]), 1.0) * w[n] + state = eq.step(tuple(state), src, dt) + out.append(state[rec_field][rec_node]) + return torch.stack(out) + + +@pytest.mark.parametrize("b", [0.0, 7e5]) +def test_the_interfaces_are_reciprocal(b): + # the two interface integrals are exact transposes and the sealed interface is the + # natural one, so the discrete Green's functions are symmetric across every interface, + # with the drag too: pressure-pressure water <-> rock and seafloor <-> rock, + # velocity-velocity rock <-> shale and seafloor (x) <-> shale (z) + f64 = torch.float64 + shape, dh, dt, nt = (41, 49), 5.0, 2e-4, 500 + prop = _prop(shape, dh, dt, None, None, f0=40.0, abcn=0, free_surface=True, dtype=f64) + mesh = prop.mesh + eq = prop._build_eq(_three_media(mesh, b)) + w = _ricker(nt, dt, 40.0, f64) + a, c, g = _node(mesh, 1, 4), _node(mesh, 4, 8, 1, 3), _node(mesh, 2, 6, 4, 2) + d, e, f = _node(mesh, 5, 3), _node(mesh, 8, 7, 3, 1), _node(mesh, 3, 9, 0, 2) + for fw, bw in (((2, a, 2, c), (2, c, 2, a)), ((2, g, 2, c), (2, c, 2, g)), + ((1, d, 1, e), (1, e, 1, d)), ((0, f, 1, e), (1, e, 0, f))): + x, y = _green(eq, dt, nt, *fw, w), _green(eq, dt, nt, *bw, w) + assert float(x.abs().max()) > 0 + torch.testing.assert_close(x, y, rtol=0, atol=1e-10 * float(x.abs().max())) + + +def test_an_undrained_porous_layer_on_a_dry_solid_is_two_elastic_layers(): + # b -> inf: the porous layer is its undrained elastic solid and the sealed interface asks + # nothing more of it. Not to round-off, for the reason in + # test_no_relative_flow_limit_is_the_undrained_elastic + from sweep.equations.biot_params import biot_elastic + shape, dh, dt, nt = (41, 49), 10.0, 5e-4, 500 + f64 = torch.float64 + kw = dict(abcn=0, free_surface=True, dtype=f64) + biot = _prop(shape, dh, dt, ["sxx", "szz"], ["vx", "vz"], **kw) + rock = dict(k_s=36e9, rho_s=2650.0, k_f=2.25e9, rho_f=1000.0, eta=1e-3, k_d=8e9, mu=6e9, + phi=0.25, kappa=1e-13, tortuosity=2.0) + poro = [float(c) for c in biot_from_rock(**rock)] + poro[7] = 1e16 + coef = _layered(biot.mesh, [0.0, 200.0], [poro, biot_elastic(*SHALE)]) + h, mu, rho = poro[0], poro[1], poro[4] + el = _layered(biot.mesh, [0.0, 200.0], [((h / rho) ** 0.5, (mu / rho) ** 0.5, rho), SHALE]) + src, rec = np.array([[24, 18]]), np.stack([np.arange(4, 46, 3), np.full(14, 30)], 1)[None] + w = _ricker(nt, dt, 12.0, f64) + with torch.no_grad(): + db = biot(w, src, rec, models=coef) + de = _prop(shape, dh, dt, ["sxx", "szz"], ["vx", "vz"], eq=ElasticVSSEM2DDeformed, **kw)( + w, src, rec, models=el) + err = float((db - de).abs().max() / de.abs().max()) + assert err < 1e-2, err + + +@pytest.mark.parametrize("edit, match", [ + (lambda c, e: c[4].__setitem__((e, 0, 0), 1000.0), "on some of its nodes only"), + (lambda c, e: [c[k].__setitem__(e, 0.0) for k in (3, 4)], "neither a solid frame"), + (lambda c, e: c[3].__setitem__(e, -1.0), "M must be >= 0"), + (lambda c, e: c[1].__setitem__(e, 1e9), r"fluid element \(rho = 0\) has nonzero \['mu'\]"), + (lambda c, e: c[6].__setitem__(e, 0.0), "m must be > 0"), +]) +def test_media_are_read_off_the_coefficients(edit, match): + p = _prop((41, 49), 5.0, 2e-4, None, None, abcn=0, free_surface=True, dtype=torch.float64) + coef = _three_media(p.mesh) + edit(coef, 0) # element 0 is water + with pytest.raises(ValueError, match=match): + p._build_eq(coef) + + +def test_a_dry_element_takes_no_fluid_coefficients(): + p = _prop((41, 49), 5.0, 2e-4, None, None, abcn=0, free_surface=True, dtype=torch.float64) + coef = _three_media(p.mesh) + coef[2][-1] = 1e9 # C in the last (dry) element + with pytest.raises(ValueError, match=r"dry elastic element \(M = 0\) has nonzero \['C'\]"): + p._build_eq(coef) + + +def _patterns(names, src): + from sweep.sem.prop_base import _Acq + p = _prop((41, 49), 5.0, 2e-4, names, ["vz"], abcn=0, free_surface=True, dtype=torch.float64) + eq = p._build_eq(_three_media(p.mesh)) + w = _ricker(40, 2e-4, 40.0, torch.float64) + p._set_nt(w) + (prop,) = p._build_props(eq, None, _Acq(p, w, np.array([src]), np.array([[[6, 6]]]), None), + "warn") + return eq, prop._build_source_patterns(eq.mesh.n_global, "cpu", torch.float64) + + +@pytest.mark.parametrize("names, src", [(["sxx", "szz"], [21, 5]), (["vz"], [21, 5]), + (["p"], [21, 33]), (["qz"], [21, 33])]) +def test_a_source_needs_its_fields_in_its_medium(names, src): + # grid row 5 is in the water, row 33 in the shale + with pytest.raises(ValueError, match="none of which lives there"): + _patterns(names, src) + + +def test_a_force_in_the_water_moves_the_water_only(): + # fz just above the seafloor: no load reaches the rock's nodes on the seafloor + eq, pats = _patterns(["fz"], [21, 11]) + assert not pats[1].any() and pats[4].any() + eq, pats = _patterns(["fz"], [21, 13]) # just below: the rock, both + assert pats[1].any() and pats[4].any() + + +def test_a_stress_source_on_the_seafloor_loads_the_rock_alone(): + # on the water-rock edge the moment tensor takes the rock's one-sided derivative in full, + # as on the free surface, not half of it + from sweep.sem.propagator import _stamp_components, build_moment_stamp + eq, pats = _patterns(["sxx", "szz"], [21, 12]) + stamp = build_moment_stamp(eq.mesh, 21 * 5.0, 12 * 5.0) + rock = [s for s in stamp.parts if bool(eq._solid_e[s.elem_id])] + assert len(stamp.parts) == 2 and len(rock) == 1 + comps = _stamp_components(rock[0], "moment_tensor", (1.0, 1.0, 0.0)) + g = eq.mesh.ibool[rock[0].elem_id].reshape(-1) + for k in (0, 1): + want = torch.zeros(eq.mesh.n_global, dtype=torch.float64).index_add( + 0, g, comps[k].reshape(-1).double()) + torch.testing.assert_close(pats[k], want) + + +@pytest.mark.parametrize("b", [0.0, 7e5]) +def test_interfaces_through_the_pml_stay_stable(b): + # the seafloor and the rock-shale interface cross the side bands with the plain + # coupling: the records keep decaying (4 s runs: to 3% of the peak with the inviscid + # slow wave, 0.1% with b = 7e5, never growing); here 1.2 s through five elements of + # C-PML, the window maxima of the last 0.2 s against 0.4 - 0.6 s: 0.50 / 0.53 (b = 0), + # 0.19 / 0.05 (ten elements: 0.48 / 0.53, 0.28 / 0.04) + f64 = torch.float64 + shape, dh, dt, nt = (41, 49), 5.0, 2e-4, 6000 + prop = _prop(shape, dh, dt, ["p"], ["p", "vz"], f0=40.0, free_surface=True, dtype=f64, + abcn=5 * NO) + w = _ricker(nt, dt, 40.0, f64) + with torch.no_grad(): + out = prop(w, np.array([[24, 6]]), np.array([[[10, 6], [30, 20]]]), + models=_three_media(prop.mesh, b))[0] + win = out.abs().reshape(6, 1000, 2, 2).amax(1) + p_water, vz_rock = win[:, 0, 0], win[:, 1, 1] + for x in (p_water, vz_rock): + assert float(x[0]) > 0 and bool(torch.isfinite(x).all()) + assert float(x[-1]) < 0.75 * float(x[2]) and float(x[-1]) < 0.15 * float(x.max()), x + + +def test_gradient_through_the_three_media(): + # one coefficient of each medium; the coefficients a medium does not have get none, so + # an update leaves the media as they are + f64 = torch.float64 + shape, dh, dt, nt = (41, 49), 5.0, 2e-4, 500 + prop = _prop(shape, dh, dt, ["p"], ["p", "vx", "vz"], f0=40.0, free_surface=True, dtype=f64, + abcn=THIN_ABCN) + coef = _three_media(prop.mesh, 7e5) + src, rec = np.array([[24, 6]]), np.stack([np.arange(6, 44, 4), np.full(10, 20)], 1)[None] + w = _ricker(nt, dt, 40.0, f64) + obs = prop(w, src, rec, models=coef).detach() + E, n = prop.mesh.n_elements, prop.mesh.ngll + zc = prop.mesh.z_e.mean((1, 2)) + xc = prop._coords([24, 0])[0] # mid-model, the C-PML aside + + def J(models): + return 0.5 * (prop(w, src, rec, models=models) - 1.05 * obs).square().sum() + + m0 = [c.clone().requires_grad_() for c in coef] + J(m0).backward() + water, dry = zc < 60, zc > 140 + for k in (0, 1, 2, 4, 5, 7): # H, mu, C, rho, rho_f, b + assert not m0[k].grad[water].any(), k + for k in (2, 3, 5, 6, 7): # C, M, rho_f, m, b + assert not m0[k].grad[dry].any(), k + for which, (z0, z1) in ((3, (0, 60)), (1, (60, 140)), (0, (140, 1e9))): # M, mu, H + sel = ((zc > z0) & (zc < z1)).to(f64).view(E, 1, 1).expand(E, n, n) + bump = torch.exp(-((prop.mesh.x_e - xc) ** 2) / 800.0) * sel + _fd_matches(J, coef, which, m0[which].grad, FD_STEP * coef[which] * bump) + + +def test_the_pml_speed_takes_every_medium(): + # Dai's fast P wave (3210 m/s) over the shale's 3000 and the water's 1500 + p = _prop((41, 49), 5.0, 2e-4, None, None, abcn=0, free_surface=True, dtype=torch.float64) + coef = _three_media(p.mesh) + fast = biot_fast_p_velocity(biot_from_biot_constants(**DAI1, b=0.0, phi=0.3)) + assert fast > SHALE[0] and p._pml_speed(coef) == pytest.approx(fast) + zc = p.mesh.z_e.mean((1, 2)) + no_rock = (zc < 60) | (zc > 140) + assert p._pml_speed([c[no_rock] for c in coef]) == pytest.approx(SHALE[0]) + assert p._pml_speed([c[zc < 60] for c in coef]) == pytest.approx((WATER[0] / WATER[1]) ** 0.5) diff --git a/test/test_cli_show.py b/test/test_cli_show.py index 3ffb877..0ee0315 100644 --- a/test/test_cli_show.py +++ b/test/test_cli_show.py @@ -99,3 +99,14 @@ def test_a_typo_is_offered_the_sem_name(capsys): out = capsys.readouterr().out assert "Did you mean" in out and "AcousticPotentialSEM2D" in out assert "spectral-element" in out + + +def test_biot_shows_propsem_source_names(capsys): + """Biot names its PropSEM sources (a force, a moment, a pore-pressure or + flux injection), which are not its source fields.""" + assert list_wavefields("BiotVSSEM2DDeformed") == 0 + out = capsys.readouterr().out + line = next(l for l in out.splitlines() if "Source types" in l) + assert all(f"'{n}'" in line for n in ("fz", "p", "sxx", "qz")) + assert ("Needed models: ['H', 'mu', 'C', 'M', 'rho', 'rho_f', 'm', 'b']" + in out)