Merged PR 86: Version 0.1.0
This commit is contained in:
Gianluca Campanella
Родитель
a2dd264caf
Коммит
b99f5de00d
|
|
@ -0,0 +1,6 @@
|
|||
{
|
||||
"python.formatting.provider": "black",
|
||||
"python.linting.enabled": true,
|
||||
"python.linting.flake8Enabled": true,
|
||||
"python.linting.pylintEnabled": false,
|
||||
}
|
||||
|
|
@ -0,0 +1,9 @@
|
|||
# Microsoft Open Source Code of Conduct
|
||||
|
||||
This project has adopted the [Microsoft Open Source Code of Conduct](https://opensource.microsoft.com/codeofconduct/).
|
||||
|
||||
Resources:
|
||||
|
||||
- [Microsoft Open Source Code of Conduct](https://opensource.microsoft.com/codeofconduct/)
|
||||
- [Microsoft Code of Conduct FAQ](https://opensource.microsoft.com/codeofconduct/faq/)
|
||||
- Contact [opencode@microsoft.com](mailto:opencode@microsoft.com) with questions or concerns
|
||||
|
|
@ -1,6 +1,6 @@
|
|||
MIT License
|
||||
|
||||
Copyright (c) Microsoft Corporation. All rights reserved.
|
||||
Copyright (c) Microsoft Corporation.
|
||||
|
||||
Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
of this software and associated documentation files (the "Software"), to deal
|
||||
|
|
@ -1 +1,9 @@
|
|||
# DeepSeismic
|
||||
|
||||
## Contributing
|
||||
|
||||
This project welcomes contributions and suggestions. Most contributions require you to agree to a Contributor License Agreement (CLA) declaring that you have the right to, and actually do, grant us the rights to use your contribution. For details, visit https://cla.opensource.microsoft.com.
|
||||
|
||||
When you submit a pull request, a CLA bot will automatically determine whether you need to provide a CLA and decorate the PR appropriately (e.g., status check, comment). Simply follow the instructions provided by the bot. You will only need to do this once across all repos using our CLA.
|
||||
|
||||
This project has adopted the [Microsoft Open Source Code of Conduct](https://opensource.microsoft.com/codeofconduct/). For more information see the [Code of Conduct FAQ](https://opensource.microsoft.com/codeofconduct/faq/) or contact [opencode@microsoft.com](mailto:opencode@microsoft.com) with any additional questions or comments.
|
||||
|
|
|
|||
|
|
@ -0,0 +1,6 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
from deepseismic import cli
|
||||
|
||||
if __name__ == "__main__":
|
||||
cli.main()
|
||||
|
|
@ -0,0 +1,3 @@
|
|||
from . import cli, forward, velocity
|
||||
|
||||
__all__ = ["cli", "forward", "velocity"]
|
||||
|
|
@ -0,0 +1,21 @@
|
|||
from functools import partial
|
||||
|
||||
import click
|
||||
|
||||
from . import forward, velocity
|
||||
|
||||
click.option = partial(click.option, show_default=True)
|
||||
|
||||
|
||||
@click.group()
|
||||
@click.pass_context
|
||||
def cli(ctx):
|
||||
ctx.ensure_object(dict)
|
||||
|
||||
|
||||
cli.add_command(forward.fwd)
|
||||
cli.add_command(velocity.vp)
|
||||
|
||||
|
||||
def main():
|
||||
cli(obj={})
|
||||
|
|
@ -0,0 +1,123 @@
|
|||
from functools import partial
|
||||
|
||||
import click
|
||||
import h5py
|
||||
import numpy as np
|
||||
|
||||
from ..forward import Receiver, RickerSource, TimeAxis, VelocityModel
|
||||
|
||||
click.option = partial(click.option, show_default=True)
|
||||
|
||||
|
||||
@click.group()
|
||||
@click.argument("input", type=click.Path())
|
||||
@click.argument("output", type=click.Path())
|
||||
@click.option(
|
||||
"-d",
|
||||
"--duration",
|
||||
default=1000.0,
|
||||
type=float,
|
||||
help="Simulation duration (in ms)",
|
||||
)
|
||||
@click.option("-dt", default=2.0, type=float, help="Time increment (in ms)")
|
||||
@click.option(
|
||||
"--n-pml", default=10, type=int, help="PML size (in grid points)"
|
||||
)
|
||||
@click.option(
|
||||
"--n-receivers",
|
||||
default=11,
|
||||
type=int,
|
||||
help="Number of receivers per horizontal dimension",
|
||||
)
|
||||
@click.option("--space-order", default=2, type=int, help="Space order")
|
||||
@click.option(
|
||||
"--spacing", default=10.0, type=float, help="Spacing between grid points"
|
||||
)
|
||||
@click.pass_context
|
||||
def fwd(
|
||||
ctx,
|
||||
dt: float,
|
||||
duration: float,
|
||||
input: str,
|
||||
n_pml: int,
|
||||
n_receivers: int,
|
||||
output: str,
|
||||
space_order: int,
|
||||
spacing: float,
|
||||
):
|
||||
"""Forward modelling"""
|
||||
if dt:
|
||||
ctx.obj["dt"] = dt
|
||||
ctx.obj["duration"] = duration
|
||||
ctx.obj["input_file"] = h5py.File(input, mode="r")
|
||||
ctx.obj["n_pml"] = n_pml
|
||||
ctx.obj["n_receivers"] = n_receivers
|
||||
ctx.obj["output_file"] = h5py.File(output, mode="w")
|
||||
ctx.obj["space_order"] = space_order
|
||||
ctx.obj["spacing"] = spacing
|
||||
|
||||
|
||||
@fwd.command()
|
||||
@click.option(
|
||||
"-f0", default=0.01, type=float, help="Source peak frequency (in kHz)"
|
||||
)
|
||||
@click.pass_context
|
||||
def ricker(ctx, f0: float):
|
||||
"""Ricker source"""
|
||||
input_file = ctx.obj["input_file"]
|
||||
output_file = ctx.obj["output_file"]
|
||||
n = sum(len(x.values()) for x in input_file.values())
|
||||
with click.progressbar(length=n) as bar:
|
||||
for input_group_name, input_group in input_file.items():
|
||||
for dataset in input_group.values():
|
||||
first_dataset = dataset
|
||||
break
|
||||
model = VelocityModel(
|
||||
shape=first_dataset.shape,
|
||||
origin=tuple(0.0 for _ in first_dataset.shape),
|
||||
spacing=tuple(ctx.obj["spacing"] for _ in first_dataset.shape),
|
||||
vp=first_dataset[()],
|
||||
space_order=ctx.obj["space_order"],
|
||||
n_pml=ctx.obj["n_pml"],
|
||||
)
|
||||
time_range = TimeAxis(
|
||||
start=0.0, stop=ctx.obj["duration"], step=ctx.obj["dt"]
|
||||
)
|
||||
source = RickerSource(
|
||||
name="source",
|
||||
grid=model.grid,
|
||||
f0=f0,
|
||||
npoint=1,
|
||||
time_range=time_range,
|
||||
)
|
||||
source.coordinates.data[0, :] = np.array(model.domain_size) * 0.5
|
||||
source.coordinates.data[0, -1] = 0.0
|
||||
n_receivers = ctx.obj["n_receivers"]
|
||||
total_receivers = n_receivers ** (len(model.shape) - 1)
|
||||
receivers = Receiver(
|
||||
name="receivers",
|
||||
grid=model.grid,
|
||||
npoint=total_receivers,
|
||||
time_range=time_range,
|
||||
)
|
||||
receivers_coords = np.meshgrid(
|
||||
*(
|
||||
np.linspace(start=0, stop=s, num=n_receivers + 2)[1:-1]
|
||||
for s in model.domain_size[:-1]
|
||||
)
|
||||
)
|
||||
for d in range(len(receivers_coords)):
|
||||
receivers.coordinates.data[:, d] = receivers_coords[
|
||||
d
|
||||
].flatten()
|
||||
receivers.coordinates.data[:, -1] = 0.0
|
||||
output_group = output_file.create_group(input_group_name)
|
||||
for input_dataset_name, vp in input_group.items():
|
||||
model.vp = vp[()]
|
||||
seismograms = model.solve(
|
||||
source=source, receivers=receivers, time_range=time_range
|
||||
)
|
||||
output_group.create_dataset(
|
||||
input_dataset_name, data=seismograms
|
||||
)
|
||||
bar.update(1)
|
||||
|
|
@ -0,0 +1,96 @@
|
|||
from functools import partial
|
||||
from itertools import islice
|
||||
from typing import Tuple
|
||||
|
||||
import click
|
||||
import h5py
|
||||
|
||||
from ..velocity import RoethTarantolaGenerator
|
||||
|
||||
click.option = partial(click.option, show_default=True)
|
||||
|
||||
|
||||
@click.group()
|
||||
@click.argument("output", type=click.Path())
|
||||
@click.option(
|
||||
"--append/--no-append",
|
||||
default=False,
|
||||
help="Whether to append to output file",
|
||||
)
|
||||
@click.option("-n", default=1, type=int, help="Number of simulations")
|
||||
@click.option(
|
||||
"-nx",
|
||||
default=100,
|
||||
type=int,
|
||||
help="Number of grid points along the first dimension",
|
||||
)
|
||||
@click.option(
|
||||
"-ny",
|
||||
default=100,
|
||||
type=int,
|
||||
help="Number of grid points along the second dimension",
|
||||
)
|
||||
@click.option(
|
||||
"-nz", type=int, help="Number of grid points along the third dimension"
|
||||
)
|
||||
@click.option("-s", "--seed", default=42, type=int, help="Random seed")
|
||||
@click.pass_context
|
||||
def vp(
|
||||
ctx,
|
||||
append: bool,
|
||||
n: int,
|
||||
nx: int,
|
||||
ny: int,
|
||||
nz: int,
|
||||
output: str,
|
||||
seed: int,
|
||||
):
|
||||
"""Vp simulation"""
|
||||
shape = (nx, ny)
|
||||
if nz is not None:
|
||||
shape += (nz,)
|
||||
output_file = h5py.File(output, mode=("a" if append else "w"))
|
||||
output_group = output_file.create_group(
|
||||
str(max((int(x) for x in output_file.keys()), default=-1) + 1)
|
||||
)
|
||||
ctx.obj["n"] = n
|
||||
ctx.obj["output_file"] = output_file
|
||||
ctx.obj["output_group"] = output_group
|
||||
ctx.obj["seed"] = seed
|
||||
ctx.obj["shape"] = shape
|
||||
|
||||
|
||||
@vp.command()
|
||||
@click.option("--n-layers", default=8, type=int, help="Number of layers")
|
||||
@click.option(
|
||||
"--initial-vp",
|
||||
default=(1350.0, 1650.0),
|
||||
type=(float, float),
|
||||
help="Initial Vp (in km/s)",
|
||||
)
|
||||
@click.option(
|
||||
"--vp-perturbation",
|
||||
default=(-190.0, 570.0),
|
||||
type=(float, float),
|
||||
help="Per-layer Vp perturbation (in km/s)",
|
||||
)
|
||||
@click.pass_context
|
||||
def rt(
|
||||
ctx,
|
||||
initial_vp: Tuple[float, float],
|
||||
n_layers: int,
|
||||
vp_perturbation: Tuple[float, float],
|
||||
):
|
||||
"""Röth-Tarantola model"""
|
||||
model = RoethTarantolaGenerator(
|
||||
shape=ctx.obj["shape"],
|
||||
seed=ctx.obj["seed"],
|
||||
n_layers=n_layers,
|
||||
initial_vp=initial_vp,
|
||||
vp_perturbation=vp_perturbation,
|
||||
)
|
||||
group = ctx.obj["output_group"]
|
||||
with click.progressbar(length=ctx.obj["n"]) as bar:
|
||||
for i, data in enumerate(islice(model.generate_many(), ctx.obj["n"])):
|
||||
group.create_dataset(str(i), data=data, compression="gzip")
|
||||
bar.update(1)
|
||||
|
|
@ -0,0 +1,14 @@
|
|||
from .models import Model, VelocityModel
|
||||
from .sources import Receiver, RickerSource, WaveletSource
|
||||
from .time import TimeAxis
|
||||
from .types import Kernel
|
||||
|
||||
__all__ = [
|
||||
"Kernel",
|
||||
"Model",
|
||||
"Receiver",
|
||||
"RickerSource",
|
||||
"TimeAxis",
|
||||
"VelocityModel",
|
||||
"WaveletSource",
|
||||
]
|
||||
|
|
@ -0,0 +1,162 @@
|
|||
from typing import Optional, Tuple, Union
|
||||
|
||||
import numpy as np
|
||||
from devito import (
|
||||
Constant,
|
||||
Eq,
|
||||
Function,
|
||||
Grid,
|
||||
Operator,
|
||||
SubDomain,
|
||||
TimeFunction,
|
||||
logger,
|
||||
solve,
|
||||
)
|
||||
|
||||
from .sources import PointSource
|
||||
from .subdomains import PhysicalDomain
|
||||
from .time import TimeAxis
|
||||
from .types import Kernel
|
||||
|
||||
logger.set_log_level("WARNING")
|
||||
|
||||
|
||||
class Model(object):
|
||||
def __init__(
|
||||
self,
|
||||
shape: Tuple[int, ...],
|
||||
origin: Tuple[float, ...],
|
||||
spacing: Tuple[float, ...],
|
||||
n_pml: Optional[int] = 0,
|
||||
dtype: Optional[type] = np.float32,
|
||||
subdomains: Optional[Tuple[SubDomain]] = (),
|
||||
):
|
||||
shape = tuple(int(x) for x in shape)
|
||||
origin = tuple(dtype(x) for x in origin)
|
||||
n_pml = int(n_pml)
|
||||
subdomains = tuple(subdomains) + (PhysicalDomain(n_pml),)
|
||||
shape_pml = tuple(x + 2 * n_pml for x in shape)
|
||||
extent_pml = tuple(s * (d - 1) for s, d in zip(spacing, shape_pml))
|
||||
origin_pml = tuple(
|
||||
dtype(o - s * n_pml) for o, s in zip(origin, spacing)
|
||||
)
|
||||
self.grid = Grid(
|
||||
shape=shape_pml,
|
||||
extent=extent_pml,
|
||||
origin=origin_pml,
|
||||
dtype=dtype,
|
||||
subdomains=subdomains,
|
||||
)
|
||||
self.n_pml = n_pml
|
||||
self.pml = Function(name="pml", grid=self.grid)
|
||||
pml_data = np.pad(
|
||||
np.zeros(shape, dtype=dtype),
|
||||
[(n_pml,) * 2 for _ in range(self.pml.ndim)],
|
||||
mode="edge",
|
||||
)
|
||||
pml_coef = 1.5 * np.log(1000.0) / 40.0
|
||||
for d in range(self.pml.ndim):
|
||||
for i in range(n_pml):
|
||||
pos = np.abs((n_pml - i + 1) / n_pml)
|
||||
val = pml_coef * (pos - np.sin(2 * np.pi * pos) / (2 * np.pi))
|
||||
idx = [slice(0, x) for x in pml_data.shape]
|
||||
idx[d] = slice(i, i + 1)
|
||||
pml_data[tuple(idx)] += val / self.grid.spacing[d]
|
||||
idx[d] = slice(
|
||||
pml_data.shape[d] - i, pml_data.shape[d] - i + 1
|
||||
)
|
||||
pml_data[tuple(idx)] += val / self.grid.spacing[d]
|
||||
pml_data = np.pad(
|
||||
pml_data,
|
||||
[(i.left, i.right) for i in self.pml._size_halo],
|
||||
mode="edge",
|
||||
)
|
||||
self.pml.data_with_halo[:] = pml_data
|
||||
self.shape = shape
|
||||
|
||||
@property
|
||||
def domain_size(self) -> Tuple[float, ...]:
|
||||
return tuple((d - 1) * s for d, s in zip(self.shape, self.spacing))
|
||||
|
||||
@property
|
||||
def dtype(self) -> type:
|
||||
return self.grid.dtype
|
||||
|
||||
@property
|
||||
def spacing(self):
|
||||
return self.grid.spacing
|
||||
|
||||
@property
|
||||
def spacing_map(self):
|
||||
return self.grid.spacing_map
|
||||
|
||||
@property
|
||||
def time_spacing(self):
|
||||
return self.grid.stepping_dim.spacing
|
||||
|
||||
|
||||
class VelocityModel(Model):
|
||||
def __init__(
|
||||
self,
|
||||
shape: Tuple[int, ...],
|
||||
origin: Tuple[float, ...],
|
||||
spacing: Tuple[float, ...],
|
||||
vp: Union[float, np.ndarray],
|
||||
space_order: Optional[int] = None,
|
||||
n_pml: Optional[int] = 0,
|
||||
dtype: Optional[type] = np.float32,
|
||||
subdomains: Optional[Tuple[SubDomain]] = (),
|
||||
):
|
||||
super().__init__(shape, origin, spacing, n_pml, dtype, subdomains)
|
||||
if isinstance(vp, np.ndarray):
|
||||
assert space_order is not None
|
||||
self.m = Function(
|
||||
name="m", grid=self.grid, space_order=int(space_order)
|
||||
)
|
||||
else:
|
||||
self.m = Constant(name="m", value=1.0 / float(vp) ** 2.0)
|
||||
self.vp = vp
|
||||
|
||||
@property
|
||||
def vp(self) -> Union[float, np.ndarray]:
|
||||
return self._vp
|
||||
|
||||
@vp.setter
|
||||
def vp(self, vp: Union[float, np.ndarray]) -> None:
|
||||
self._vp = vp
|
||||
if isinstance(vp, np.ndarray):
|
||||
pad_widths = [
|
||||
(self.n_pml + i.left, self.n_pml + i.right)
|
||||
for i in self.m._size_halo
|
||||
]
|
||||
self.m.data_with_halo[:] = np.pad(
|
||||
1.0 / self.vp ** 2.0, pad_widths, mode="edge"
|
||||
)
|
||||
else:
|
||||
self.m.data = 1.0 / float(vp) ** 2.0
|
||||
|
||||
def solve(
|
||||
self,
|
||||
source: PointSource,
|
||||
receivers: PointSource,
|
||||
time_range: TimeAxis,
|
||||
space_order: Optional[int] = 4,
|
||||
kernel: Optional[Kernel] = Kernel.OT2,
|
||||
) -> np.ndarray:
|
||||
assert isinstance(kernel, Kernel)
|
||||
u = TimeFunction(
|
||||
name="u", grid=self.grid, time_order=2, space_order=space_order
|
||||
)
|
||||
H = u.laplace
|
||||
if kernel is Kernel.OT4:
|
||||
H += self.time_spacing ** 2 / 12 * u.laplace2(1 / self.m)
|
||||
eq = Eq(
|
||||
u.forward, solve(self.m * u.dt2 - H + self.pml * u.dt, u.forward)
|
||||
)
|
||||
src_term = source.inject(
|
||||
field=u.forward, expr=source * self.time_spacing ** 2 / self.m
|
||||
)
|
||||
rec_term = receivers.interpolate(expr=u)
|
||||
op = Operator([eq] + src_term + rec_term, subs=self.spacing_map)
|
||||
op(time=time_range.num - 1, dt=time_range.step)
|
||||
return receivers.data
|
||||
|
|
@ -0,0 +1,132 @@
|
|||
from typing import Optional
|
||||
|
||||
import numpy as np
|
||||
import sympy
|
||||
from devito.types import Dimension, SparseTimeFunction
|
||||
from devito.types.basic import _SymbolCache
|
||||
from scipy import interpolate
|
||||
|
||||
from .time import TimeAxis
|
||||
|
||||
|
||||
class PointSource(SparseTimeFunction):
|
||||
def __new__(cls, *args, **kwargs):
|
||||
if cls in _SymbolCache:
|
||||
options = kwargs.get("options", {})
|
||||
obj = sympy.Function.__new__(cls, *args, **options)
|
||||
obj._cached_init()
|
||||
return obj
|
||||
name = kwargs.pop("name")
|
||||
grid = kwargs.pop("grid")
|
||||
time_range = kwargs.pop("time_range")
|
||||
time_order = kwargs.pop("time_order", 2)
|
||||
p_dim = kwargs.pop("dimension", Dimension(name="p_%s" % name))
|
||||
npoint = kwargs.pop("npoint", None)
|
||||
coordinates = kwargs.pop(
|
||||
"coordinates", kwargs.pop("coordinates_data", None)
|
||||
)
|
||||
if npoint is None:
|
||||
assert (
|
||||
coordinates is not None
|
||||
), "Either `npoint` or `coordinates` must be provided"
|
||||
npoint = coordinates.shape[0]
|
||||
obj = SparseTimeFunction.__new__(
|
||||
cls,
|
||||
name=name,
|
||||
grid=grid,
|
||||
dimensions=(grid.time_dim, p_dim),
|
||||
npoint=npoint,
|
||||
nt=time_range.num,
|
||||
time_order=time_order,
|
||||
coordinates=coordinates,
|
||||
**kwargs
|
||||
)
|
||||
obj._time_range = time_range
|
||||
data = kwargs.get("data")
|
||||
if data is not None:
|
||||
obj.data[:] = data
|
||||
return obj
|
||||
|
||||
@property
|
||||
def time_range(self) -> TimeAxis:
|
||||
return self._time_range
|
||||
|
||||
@property
|
||||
def time_values(self) -> np.ndarray:
|
||||
return self._time_range.time_values
|
||||
|
||||
def resample(
|
||||
self,
|
||||
dt: Optional[float] = None,
|
||||
num: Optional[int] = None,
|
||||
rtol: Optional[float] = 1.0e-5,
|
||||
order: Optional[int] = 3,
|
||||
):
|
||||
assert (dt is not None) ^ (
|
||||
num is not None
|
||||
), "Exactly one of `dt` or `num` must be provided"
|
||||
start = self._time_range.start
|
||||
stop = self._time_range.stop
|
||||
dt0 = self._time_range.step
|
||||
if dt is not None:
|
||||
new_time_range = TimeAxis(start=start, stop=stop, step=dt)
|
||||
else:
|
||||
new_time_range = TimeAxis(start=start, stop=stop, num=num)
|
||||
dt = new_time_range.step
|
||||
if np.isclose(dt0, dt, rtol=rtol):
|
||||
return self
|
||||
n_traces = self.data.shape[1]
|
||||
new_traces = np.zeros(
|
||||
(new_time_range.num, n_traces), dtype=self.data.dtype
|
||||
)
|
||||
for j in range(n_traces):
|
||||
tck = interpolate.splrep(
|
||||
self._time_range.time_values, self.data[:, j], k=order
|
||||
)
|
||||
new_traces[:, j] = interpolate.splev(
|
||||
new_time_range.time_values, tck
|
||||
)
|
||||
return PointSource(
|
||||
name=self.name,
|
||||
grid=self.grid,
|
||||
time_range=new_time_range,
|
||||
coordinates=self.coordinates.data,
|
||||
data=new_traces,
|
||||
)
|
||||
|
||||
_pickle_kwargs = SparseTimeFunction._pickle_kwargs + ["time_range"]
|
||||
_pickle_kwargs.remove("nt") # Inferred from time_range
|
||||
|
||||
|
||||
class Receiver(PointSource):
|
||||
pass
|
||||
|
||||
|
||||
class WaveletSource(PointSource):
|
||||
def __new__(cls, *args, **kwargs):
|
||||
if cls in _SymbolCache:
|
||||
options = kwargs.get("options", {})
|
||||
obj = sympy.Function.__new__(cls, *args, **options)
|
||||
obj._cached_init()
|
||||
return obj
|
||||
npoint = kwargs.pop("npoint", 1)
|
||||
obj = PointSource.__new__(cls, npoint=npoint, **kwargs)
|
||||
obj.f0 = kwargs.get("f0")
|
||||
for p in range(npoint):
|
||||
obj.data[:, p] = obj.wavelet(obj.f0, obj.time_values)
|
||||
return obj
|
||||
|
||||
def __init__(self, *args, **kwargs):
|
||||
if not self._cached():
|
||||
super(WaveletSource, self).__init__(*args, **kwargs)
|
||||
|
||||
def wavelet(self, f0: float, t: np.ndarray) -> np.ndarray:
|
||||
raise NotImplementedError
|
||||
|
||||
_pickle_kwargs = PointSource._pickle_kwargs + ["f0"]
|
||||
|
||||
|
||||
class RickerSource(WaveletSource):
|
||||
def wavelet(self, f0: float, t: np.ndarray) -> np.ndarray:
|
||||
r = np.pi * f0 * (t - 1.0 / f0)
|
||||
return (1.0 - 2.0 * r ** 2.0) * np.exp(-r ** 2.0)
|
||||
|
|
@ -0,0 +1,16 @@
|
|||
from typing import Dict, Iterable, Tuple
|
||||
|
||||
from devito import Dimension, SubDomain
|
||||
|
||||
|
||||
class PhysicalDomain(SubDomain):
|
||||
name = "physical_domain"
|
||||
|
||||
def __init__(self, n_pml: int):
|
||||
super().__init__()
|
||||
self.n_pml = n_pml
|
||||
|
||||
def define(
|
||||
self, dimensions: Iterable[Dimension]
|
||||
) -> Dict[Dimension, Tuple[str, int, int]]:
|
||||
return {d: ("middle", self.n_pml, self.n_pml) for d in dimensions}
|
||||
|
|
@ -0,0 +1,34 @@
|
|||
from typing import Optional
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
class TimeAxis(object):
|
||||
def __init__(
|
||||
self,
|
||||
start: Optional[float] = None,
|
||||
stop: Optional[float] = None,
|
||||
num: Optional[int] = None,
|
||||
step: Optional[float] = None,
|
||||
dtype: Optional[type] = np.float32,
|
||||
):
|
||||
if start is None:
|
||||
start = step * (1 - num) + stop
|
||||
elif stop is None:
|
||||
stop = step * (num - 1) + start
|
||||
elif num is None:
|
||||
num = int(np.ceil((stop - start + step) / step))
|
||||
stop = step * (num - 1) + start
|
||||
elif step is None:
|
||||
step = (stop - start) / (num - 1)
|
||||
else:
|
||||
raise ValueError
|
||||
self.start = start
|
||||
self.stop = stop
|
||||
self.num = num
|
||||
self.step = step
|
||||
self.dtype = dtype
|
||||
|
||||
@property
|
||||
def time_values(self) -> np.ndarray:
|
||||
return np.linspace(self.start, self.stop, self.num, dtype=self.dtype)
|
||||
|
|
@ -0,0 +1,6 @@
|
|||
from enum import Enum, auto
|
||||
|
||||
|
||||
class Kernel(Enum):
|
||||
OT2 = auto()
|
||||
OT4 = auto()
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
from .generator import Generator
|
||||
from .roeth_tarantola import RoethTarantolaGenerator
|
||||
|
||||
__all__ = ["Generator", "RoethTarantolaGenerator"]
|
||||
|
|
@ -0,0 +1,22 @@
|
|||
from typing import Optional, Tuple
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
class Generator(object):
|
||||
def __init__(
|
||||
self,
|
||||
shape: Tuple[int, ...],
|
||||
dtype: Optional[type] = np.float32,
|
||||
seed: Optional[int] = None,
|
||||
):
|
||||
self.shape = shape
|
||||
self.dtype = dtype
|
||||
self._prng = np.random.RandomState(seed)
|
||||
|
||||
def generate(self) -> np.ndarray:
|
||||
raise NotImplementedError
|
||||
|
||||
def generate_many(self) -> np.ndarray:
|
||||
while True:
|
||||
yield self.generate()
|
||||
|
|
@ -0,0 +1,41 @@
|
|||
from typing import Optional, Tuple
|
||||
|
||||
import numpy as np
|
||||
|
||||
from .generator import Generator
|
||||
|
||||
|
||||
class RoethTarantolaGenerator(Generator):
|
||||
def __init__(
|
||||
self,
|
||||
shape: Tuple[int, ...],
|
||||
dtype: Optional[type] = np.float32,
|
||||
seed: Optional[int] = None,
|
||||
depth_dim: Optional[int] = -1,
|
||||
n_layers: Optional[int] = 8,
|
||||
initial_vp: Optional[Tuple[float, float]] = (1.35, 1.65),
|
||||
vp_perturbation: Optional[Tuple[float, float]] = (-0.19, 0.57),
|
||||
):
|
||||
super().__init__(shape, dtype, seed)
|
||||
self.depth_dim = depth_dim
|
||||
self.n_layers = n_layers
|
||||
self.initial_vp = initial_vp
|
||||
self.vp_perturbation = vp_perturbation
|
||||
|
||||
def generate(self) -> np.ndarray:
|
||||
vp = np.zeros(self.shape, dtype=self.dtype)
|
||||
dim = self.depth_dim
|
||||
layer_idx = np.round(
|
||||
np.linspace(0, self.shape[dim], self.n_layers + 1)
|
||||
).astype(np.int)
|
||||
vp_idx = [slice(0, x) for x in vp.shape]
|
||||
layer_vp = None
|
||||
for i in range(self.n_layers):
|
||||
vp_idx[dim] = slice(layer_idx[i], layer_idx[i + 1])
|
||||
layer_vp = (
|
||||
self._prng.uniform(*self.initial_vp)
|
||||
if layer_vp is None
|
||||
else layer_vp + self._prng.uniform(*self.vp_perturbation)
|
||||
)
|
||||
vp[tuple(vp_idx)] = layer_vp
|
||||
return vp
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
[tool.black]
|
||||
line-length = 79
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
[aliases]
|
||||
test=pytest
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
import setuptools
|
||||
|
||||
with open("README.md", "r") as f:
|
||||
long_description = f.read()
|
||||
|
||||
setuptools.setup(
|
||||
author="DeepSeismic Maintainers",
|
||||
author_email="deepseismic@microsoft.com",
|
||||
classifiers=[
|
||||
"Intended Audience :: Developers",
|
||||
"Intended Audience :: Science/Research",
|
||||
"License :: OSI Approved :: MIT License",
|
||||
"Operating System :: OS Independent",
|
||||
"Programming Language :: Python :: 3",
|
||||
"Programming Language :: Python :: 3.5",
|
||||
"Programming Language :: Python :: 3.6",
|
||||
"Programming Language :: Python :: 3.7",
|
||||
"Topic :: Scientific/Engineering",
|
||||
"Topic :: Software Development",
|
||||
],
|
||||
dependency_links=[
|
||||
"https://github.com/opesci/devito/archive/v3.5.tar.gz#egg=devito-3.5"
|
||||
],
|
||||
description="DeepSeismic",
|
||||
install_requires=[
|
||||
"click==7.0",
|
||||
"devito==3.5",
|
||||
"h5py==2.9.0",
|
||||
"numpy==1.17.0",
|
||||
"scipy==1.3.0",
|
||||
"sympy==1.4",
|
||||
],
|
||||
license="MIT",
|
||||
long_description=long_description,
|
||||
long_description_content_type="text/markdown",
|
||||
name="deepseismic",
|
||||
packages=setuptools.find_packages(
|
||||
include=["deepseismic", "deepseismic.*"]
|
||||
),
|
||||
platforms="any",
|
||||
python_requires=">= 3.5",
|
||||
scripts=["bin/ds"],
|
||||
setup_requires=["pytest-runner"],
|
||||
tests_require=["pytest"],
|
||||
url="https://github.com/microsoft/deepseismic",
|
||||
version="0.1.0",
|
||||
zip_safe=False,
|
||||
)
|
||||
Загрузка…
Ссылка в новой задаче