fix high order stability

Project.toml

@@ -1,7 +1,7 @@
 name = "JUDI"
 uuid = "f3b833dc-6b2e-5b9c-b940-873ed6319979"
 authors = ["Philipp Witte, Mathias Louboutin"]
-version = "3.3.5"
+version = "3.3.6"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

examples/scripts/modeling_basic_2D.jl

@@ -91,7 +91,7 @@ q = judiVector(srcGeometry, wavelet)
 #' condition for the propagation.
 
 # Setup options
-opt = Options(subsampling_factor=2, dt_comp=1.0)
+opt = Options(subsampling_factor=2, space_order=32)
 
 #' Linear Operators
 #' The core idea behind JUDI is to abstract seismic inverse problems in term of linear algebra. In its simplest form, seismic inversion can be formulated as

src/TimeModeling/Types/ModelStructure.jl

@@ -170,6 +170,7 @@ end
 
 copy(x::PhysicalParameter{T, N}) where {T<:Real, N} = PhysicalParameter(x.n, x.d, x.o, copy(x.data))
 unsafe_convert(::Type{Ptr{T}}, p::PhysicalParameter{T, N}) where {T<:Real, N} = unsafe_convert(Ptr{T}, p.data)
+Base.Vector{T}(m::PhysicalParameter) where T = Vector{T}(m.data[:])
 
 # Equality
 ==(A::PhysicalParameter{T, N}, B::PhysicalParameter{T, N}) where {T<:Real, N} = (A.data == B.data && A.o == B.o && A.d == B.d)

src/pysource/models.py

@@ -24,6 +24,9 @@ def getmax(f):
         return np.max(f)
 
 
+_thomsen = [('epsilon', 1), ('delta', 1), ('theta', 0), ('phi', 0)]
+
+
 class PhysicalDomain(SubDomain):
 
     name = 'nofsdomain'
@@ -157,7 +160,7 @@ class Model(object):
     dt: Float
         User provided computational time-step
     """
-    def __init__(self, origin, spacing, shape, space_order=2, nbl=40, dtype=np.float32,
+    def __init__(self, origin, spacing, shape, space_order=8, nbl=40, dtype=np.float32,
                  m=None, epsilon=None, delta=None, theta=None, phi=None, rho=None,
                  b=None, qp=None, lam=None, mu=None, dm=None, fs=False, **kwargs):
         # Setup devito grid
@@ -267,7 +270,13 @@ def physical_params(self, **kwargs):
 
     @property
     def zero_thomsen(self):
-        return {self.epsilon: 1, self.delta: 1, self.theta: 0, self.phi: 0}
+        out = {}
+        for (t, v) in _thomsen:
+            try:
+                out.update({getattr(self, t): v})
+            except AttributeError:
+                pass
+        return out
 
     @switchconfig(log_level='ERROR')
     def _gen_phys_param(self, field, name, space_order, is_param=False,
@@ -397,12 +406,14 @@ def _cfl_coeff(self):
         """
         # Elasic coefficient (see e.g )
         if self.is_elastic:
-            coeffs = fd_w(1, range(-4, 5), .5)
+            so = max(self._space_order // 2, 2)
+            coeffs = fd_w(1, range(-so, so), .5)
             c_fd = sum(np.abs(coeffs[-1][-1])) / 2
             return np.sqrt(self.dim) / self.dim / c_fd
         a1 = 4  # 2nd order in time
-        coeffs = fd_w(2, range(-8, 9), 0)[-1][-1]
-        return np.sqrt(a1/float(self.grid.dim * sum(np.abs(coeffs))))
+        so = max(self._space_order // 2, 4)
+        coeffs = fd_w(2, range(-so, so), 0)[-1][-1]
+        return .9 * np.sqrt(a1/float(self.grid.dim * sum(np.abs(coeffs))))
 
     @property
     def _thomsen_scale(self):
@@ -526,14 +537,15 @@ def __init__(self, tti, visco, elastic, spacing, fs, space_order, p_params):
         self.is_elastic = elastic
         self.spacing = spacing
         self.fs = fs
+        N = 2 * space_order + 1
         if fs:
-            fsdomain = FSDomain(space_order + 1)
-            physdomain = PhysicalDomain(space_order + 1, fs=fs)
+            fsdomain = FSDomain(N)
+            physdomain = PhysicalDomain(N, fs=fs)
             subdomains = (physdomain, fsdomain)
         else:
             subdomains = ()
-        self.grid = Grid(tuple([space_order+1]*len(spacing)),
-                         extent=[s*space_order for s in spacing],
+        self.grid = Grid(tuple([N]*len(spacing)),
+                         extent=[s*(N-1) for s in spacing],
                          subdomains=subdomains)
         self.dimensions = self.grid.dimensions
 
@@ -564,3 +576,13 @@ def dim(self):
         Spatial dimension of the problem and model domain.
         """
         return self.grid.dim
+
+    @property
+    def zero_thomsen(self):
+        out = {}
+        for (t, v) in _thomsen:
+            try:
+                out.update({getattr(self, t): v})
+            except AttributeError:
+                pass
+        return out

src/pysource/operators.py

@@ -94,7 +94,7 @@ def forward_op(p_params, tti, visco, elas, space_order, fw, spacing, save, t_sub
     u = wavefield(model, space_order, save=save, nt=nt, t_sub=t_sub, fw=fw)
 
     # Expression for saving wavefield if time subsampling is used
-    eq_save = save_subsampled(model, u, nt, t_sub)
+    eq_save = save_subsampled(model, u, nt, t_sub, space_order=space_order)
 
     # Extended source
     q = extented_src(model, wsrc, wavelet, q=q)
@@ -153,7 +153,7 @@ def born_op(p_params, tti, visco, elas, space_order, fw, spacing, save, pt_src,
     ul = wavefield(model, space_order, name="l", fw=fw)
 
     # Expression for saving wavefield if time subsampling is used
-    eq_save = save_subsampled(model, u, nt, t_sub)
+    eq_save = save_subsampled(model, u, nt, t_sub, space_order=space_order)
 
     # Add extended source
     q = extented_src(model, wsrc, wavelet)

src/pysource/propagators.py

@@ -44,7 +44,7 @@ def forward(model, src_coords, rcv_coords, wavelet, space_order=8, save=False,
     f0q = Constant('f0', value=f0) if model.is_viscoacoustic else None
 
     # Expression for saving wavefield if time subsampling is used
-    u_save = wavefield_subsampled(model, u, nt, t_sub)
+    u_save = wavefield_subsampled(model, u, nt, t_sub, space_order=space_order)
 
     # Illumination
     I = illumination(u, illum)
@@ -179,7 +179,7 @@ def born(model, src_coords, rcv_coords, wavelet, space_order=8, save=False,
     kw = base_kwargs(model.critical_dt)
     f0q = Constant('f0', value=f0) if model.is_viscoacoustic else None
     # Expression for saving wavefield if time subsampling is used
-    u_save = wavefield_subsampled(model, u, nt, t_sub)
+    u_save = wavefield_subsampled(model, u, nt, t_sub, space_order=space_order)
 
     # On-the-fly Fourier
     dft_m, fr = fourier_modes(u, freq_list)