🐛 [BUG] - <AXISYM simulation missing factor $2\pi$>

[tianshi-liu]

Description

A $2\pi$ factor may be missing in the axisymmetric solver of SPECFEM2D. Attached is a comparison between SPECFEM2D AXISYM simulation and analytical solution(eq. 4.27, Aki & Richards 2002) . The fact that the result of SPECFEM2D divided by $2\pi$ exactly matches the analytical solution seems to confirm the issue. I've also roughly checked the code, and it seems that the $2\pi$ factor in eq. 5 of Nissen-Meyer et al (2007) is never taken care of. Based on commit history, I think @bottero or @danielpeter can help me address this issue.

Affected SPECFEM2D version

latest development version (bfcc9b6)

Your software and hardware environment

GCC 11.5.0, Open MPI 4.1.1

Reproduction steps

1. Setup simulation based on EXAMPLES/applications/axisymmetric_examples/axisymmetric_case_AXISYM_option/
2. Change DATA/Par_file at the following lines such that the model is homogeneous

82c82
< ATTENUATION_VISCOELASTIC        = .true.         # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
---
> ATTENUATION_VISCOELASTIC        = .false.         # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
292c292
< nbmodels                        = 2
---
> nbmodels                        = 1
309,310c309,310
< 1 1 2500.d0 3400.d0 1963.d0 0 0 100 100 0 0 0 0 0 0
< 2 1 1020.d0 1500.d0 0.d0 0 0 9999 9999 0 0 0 0 0 0
---
> 1 1 2500.d0 3400.d0 1963.d0 0 0 9999 9999 0 0 0 0 0 0
> #2 1 1020.d0 1500.d0 0.d0 0 0 9999 9999 0 0 0 0 0 0
327c327
< nbregions                       = 2              # then set below the different regions and model number for each region
---
> nbregions                       = 1              # then set below the different regions and model number for each region
329,330c329,330
< 1 180 1   60 1
< 1 180 61 120 2
---
> 1 180 1  120 1
> #1 180 61 120 2

3. Change DATA/interfaces_axisym.dat to remove the ocean layer

4c4
<  3
---
>  2
17,19c17,19
<  2
<     0 2400
<  6400 2400
---
> # 2
> #    0 2400
> # 6400 2400
32c32
<  60
---
>  120
36c36
<  60
---
> # 60

4. Change DATA/SOURCE such that the source time function is Gaussian instead of Ricker:

< time_function_type              = 1
---
> time_function_type              = 3

5. Run mesher and solver

bash run_this_example.sh

6. Run the following code for plotting

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf

xs = 0.0
zs = 2500.0
f0 = 10.0

rho = 2500.0
vp = 3400.0
vs = 1963.0

xr = 2500.0
zr = 1500.0

dt = 0.60e-3
nt = 4500

A = np.loadtxt("OUTPUT_FILES/AA.S0003.BXZ.semd")

t = A[:,0]
seis = A[:,1] * (1.0e-10) * (2.0 * np.pi * np.pi * f0 * f0) / (2.0*np.pi) # multiply (1.0e-10) * (2.0 * np.pi * np.pi * f0 * f0) to make the source time function peaked at 1N

plt.plot(t, seis, 'g', label='SPECFEM divided by (2*pi)')

def delta(i, j):
  if i==j:
    return 1.0
  else:
    return 0.0

def gauss_stf(t_arr, tau):
  return np.exp(-0.5*t_arr**2/tau**2)/tau/np.sqrt(2.0*np.pi)

def get_radiation_pattern(gamma, comp, field):
  #gamma = dict({'x': rx, 'y': ry, 'z': rz})
  if len(comp) == 2:
    n = comp[0]
    p = comp[1]
    if field == 'FP':
      return gamma[n] * gamma[p]
    elif field == 'FS':
      return delta(n, p) - gamma[n] * gamma[p]
    elif field  == 'NF':
      return 3.0 * gamma[n] * gamma[p] - delta(n, p)
    else:
      return None
  elif len(comp) == 3:
    n = comp[0]
    p = comp[1]
    q = comp[2]
    if field == 'FP':
      return gamma[n] * gamma[p] * gamma[q]
    elif field == 'FS':
      return (delta(n, p) - gamma[n] * gamma[p]) * gamma[q]
    elif field == 'MP':
      return 6.0 * gamma[n] * gamma[p] * gamma[q] \
             - gamma[n] * delta(p, q) \
             - gamma[p] * delta(n, q) \
             - gamma[q] * delta(n, p)
    elif field == 'MS':
      return - 6.0 * gamma[n] * gamma[p] * gamma[q] \
             + gamma[n] * delta(p, q) \
             + gamma[p] * delta(n, q) \
             + 2.0 * gamma[q] * delta(n, p)
    elif field == 'NF':
      return  15.0 * gamma[n] * gamma[p] * gamma[q] \
             - 3.0 * gamma[n] * delta(p, q) \
             - 3.0 * gamma[p] * delta(n, q) \
             - 3.0 * gamma[q] * delta(n, p)
    else:
      return None
  else:
    return None
def get_fullspace_solution_time_domain(x, y, z, fx, fy, fz, la, mu, rho, dt, npts, tau=0.0, nf=False, amp=1.0, tshift=0.0):
  t_arr = np.arange(npts) * dt - tshift
  r = np.sqrt(x**2+y**2+z**2)
  rx = x/r
  ry = y/r
  rz = z/r
  gamma = dict({'x': rx, 'y': ry, 'z': rz})
  vp = np.sqrt((la + mu * 2.0) / rho)
  vs = np.sqrt(mu / rho)
  amp_p = 0.25/np.pi/(la + mu * 2.0)*amp
  amp_s = 0.25/np.pi/mu*amp
  amp_nf = 0.25/np.pi/rho*amp
  tp = r / vp
  ts = r / vs
  spec_p = gauss_stf(t_arr-tp, tau)
  spec_s = gauss_stf(t_arr-ts, tau)
  spec_nf = t_arr/2.0*(erf((ts-t_arr)/np.sqrt(2.0)/tau)-erf((tp-t_arr)/np.sqrt(2.0)/tau)) - tau*tau*(gauss_stf(t_arr-ts, tau)-gauss_stf(t_arr-tp, tau))
  fx_x = amp_p/r * spec_p * get_radiation_pattern(gamma, 'xx', 'FP') + \
         amp_s/r * spec_s * get_radiation_pattern(gamma, 'xx', 'FS')
  fy_x = amp_p/r * spec_p * get_radiation_pattern(gamma, 'xy', 'FP') + \
         amp_s/r * spec_s * get_radiation_pattern(gamma, 'xy', 'FS')
  fz_x = amp_p/r * spec_p * get_radiation_pattern(gamma, 'xz', 'FP') + \
         amp_s/r * spec_s * get_radiation_pattern(gamma, 'xz', 'FS')

  fx_y = amp_p/r * spec_p * get_radiation_pattern(gamma, 'yx', 'FP') + \
         amp_s/r * spec_s * get_radiation_pattern(gamma, 'yx', 'FS')
  fy_y = amp_p/r * spec_p * get_radiation_pattern(gamma, 'yy', 'FP') + \
         amp_s/r * spec_s * get_radiation_pattern(gamma, 'yy', 'FS')
  fz_y = amp_p/r * spec_p * get_radiation_pattern(gamma, 'yz', 'FP') + \
         amp_s/r * spec_s * get_radiation_pattern(gamma, 'yz', 'FS')

  fx_z = amp_p/r * spec_p * get_radiation_pattern(gamma, 'zx', 'FP') + \
         amp_s/r * spec_s * get_radiation_pattern(gamma, 'zx', 'FS')
  fy_z = amp_p/r * spec_p * get_radiation_pattern(gamma, 'zy', 'FP') + \
         amp_s/r * spec_s * get_radiation_pattern(gamma, 'zy', 'FS')
  fz_z = amp_p/r * spec_p * get_radiation_pattern(gamma, 'zz', 'FP') + \
         amp_s/r * spec_s * get_radiation_pattern(gamma, 'zz', 'FS')

  if nf: # include near-field
    fx_x += amp_nf/r/r/r * spec_nf * get_radiation_pattern(gamma, 'xx', 'NF')
    fy_x += amp_nf/r/r/r * spec_nf * get_radiation_pattern(gamma, 'xy', 'NF')
    fz_x += amp_nf/r/r/r * spec_nf * get_radiation_pattern(gamma, 'xz', 'NF')

    fx_y += amp_nf/r/r/r * spec_nf * get_radiation_pattern(gamma, 'yx', 'NF')
    fy_y += amp_nf/r/r/r * spec_nf * get_radiation_pattern(gamma, 'yy', 'NF')
    fz_y += amp_nf/r/r/r * spec_nf * get_radiation_pattern(gamma, 'yz', 'NF')

    fx_z += amp_nf/r/r/r * spec_nf * get_radiation_pattern(gamma, 'zx', 'NF')
    fy_z += amp_nf/r/r/r * spec_nf * get_radiation_pattern(gamma, 'zy', 'NF')
    fz_z += amp_nf/r/r/r * spec_nf * get_radiation_pattern(gamma, 'zz', 'NF')

  u = np.zeros(shape=(3,npts),dtype=float)
  u[0,:] = fx_x * fx + fy_x * fy + fz_x * fz
  u[1,:] = fx_y * fx + fy_y * fy + fz_y * fz
  u[2,:] = fx_z * fx + fy_z * fy + fz_z * fz

  return u

la2mu = rho * vp * vp
mu = rho * vs * vs
la = la2mu - 2.0 * mu


tau = 1.0/(np.pi*f0*np.sqrt(2.0)) # Gaussian defined by exp(-0.5*(t/tau)**2)
seis_ref = get_fullspace_solution_time_domain(xr-xs, 0.0, zr-zs, 0.0, 0.0, -1.0, la, mu, rho, dt, nt, tau, nf=True, amp=np.sqrt(2.0*np.pi)*tau)

t_arr = np.arange(nt) * dt

plt.plot(t_arr, seis_ref[2,:], 'r--', label='analytical solution')
plt.legend()
plt.show()

Screenshots

![DESCRIPTION](LINK.png)

Logs

OS

No response

[danielpeter]

thanks for noticing! indeed, there seems to be some factor missing when computing an axisymmetric solution. however, let me double-check and see how to fix this.

also, in your python script to compute the analytical solution, you seem to use a different source time function definition than the one used in the 2D code - but your code is going to be helpful as well for comparisons, I might just amend it a bit for this test.

(while looking into this, I noticed that there are also some confusing source time function, vector forces and source term sign conventions used - i'll try to make these formulations more similar to the 3D Cartesian version).