#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on 8 May 2025

@author: Patrick Kurzeja

Version 1.0
A simple prototype model for calculating the dispersion relations of a residually saturated porous medium with internal oscillators for the residual phase
Simplified example Python code from the original Matlab code, hope it may help your research (refer to [1], all rights reserved)

*** Literature for further reading and referencing (a limited selection): ***
The key reference for this code to use:
[1] P Kurzeja, H Steeb (2022): Acoustic waves in saturated porous media with gas bubbles, Philosophical Transactions of the Royal Society A 380 (2237), 20210370
The original idea for residual saturation with internal oscillators:
[2] M Frehner M, H Steeb, SM Schmalholz (2010): Wave velocity dispersion and attenuation in media exhibiting internal oscillations. In Wave propagation in materials for modern applications (ed A Petrin), 455–476. Vukovar, Croatia: In-Tech Education and Publishing
Dispersion relations of a gas-filled porous medium with residual liquid clusters (liquid and gas change roles, same mathematical structure, though, with more details on the derivation):
[3] H Steeb, PS Kurzeja, M Frehner, SM Schmalholz (2012): Phase velocity dispersion and attenuation of seismic waves due to trapped fluids in residual saturated porous media, Vadose Zone Journal 11 (3), vzj2011.0121
[4] P Kurzeja (2014): Waves in partially saturated porous media: An investigation on multiple scales, PhD thesis, Ruhr-University Bochum, Germany
Oscillation properties of gas bubbles:
[5] M Minnaert (1933): XVI. On musical air-bubbles and the sounds of running water, Philosophical Magazine 16, 235–248
[6] T Leighton (1994): Acoustic bubble detection - I: the detection of stable gas bodies. Environmental Engineering 7, 9–16
Oscillation properties of liquid clusters:
[7] PS Kurzeja, H Steeb (2014): Variational formulation of oscillating fluid clusters and oscillator-like classification. I. Theory, Physics of Fluids 26 (4), 042106
[8] PS Kurzeja, H Steeb (2014): Variational formulation of oscillating fluid clusters and oscillator-like classification. II. Numerical study of pinned liquid clusters, Physics of Fluids 26 (4), 042107
***
This work closely follows [1] with the mathematical structure adopted in [3,4]
***
"""



print("*******************************************************")
print("*** Start: Calculating dispersion relations ***")
print("*******************************************************")
print(" ")



# -------------------------------------------------------
# Modules
# -------------------------------------------------------
print(" ... loading modules")

import numpy as np
import matplotlib.pyplot as plt
# -------------------------------------------------------



# -------------------------------------------------------
# Material parameters
# -------------------------------------------------------
print(" ... defining material parameters")

# solid skeleton
phi_0   =0.19 # initial porosity
alpha_infty = 1 + 0.5*(1/phi_0 - 1) # tortuosity for inertial drag
G       =6.0e9 # shear modulus skeleton    
mu      =G
K       =8.0e9 # bulk moduls of skeleton
lamb  = K-(2/3)*mu
k_s     = 2e-13 # intrinsic permeability [m^2] 
rho_sR_0=2650 # true density solid [kg/m3]

# continuous fluid (water at ambient conditions)
rho_cR_0 = 1000                 # true density water [kg/m3]
eta_cR   = 0.001                  # effective dynamic viscosity water [Pa s]
K_c      = 2.00e9                 # bulk modulus of water [Pa]

# discontinuous fluid (air at ambient conditions)
rho_dR_0 = 1.02 # true density gas [kg/m3]
eta_dR = 17.1e-6 # effective dynamic viscosity gas [Pa s]
s_d_0    = 0.1 # gas saturation = vol_gas/vol_freeSpaceInSolid
gamma_ai = 1.4 # adiabatic index or heat capacity ratio or polytropic exponent
p0 = 10**5 # pressure inside the discontinuous clusters
K_d      = p0*gamma_ai # bulk modulus gas [Pa]
sigma_surf = 0.07 # surface tension [N/m]
r_osc = np.array([1e-3, 2e-3, 3e-3]) # radii of the distinct air bubbles
alpha_osc = np.array([0.3, 0.5, 0.2]) # volumetric content of the distinct air bubbles with respect to overall air volume; must sum up to 1

# derived properties
s_c_0   = 1 - s_d_0 # saturation continuous phase
n_s_0 = 1.0 - phi_0 # volume fraction solid
n_d_0 = s_d_0 * phi_0 # volume fraction discontinuous phase
n_c_0   = s_c_0 * phi_0 # volume fraction continuous phase

rho_s_0 = n_s_0 * rho_sR_0 # partial density solid phase
rho_c_0 = n_c_0 * rho_cR_0 # partial density continuous phase
rho_d_0 = n_d_0 * rho_dR_0 # partial density discontinuous phase
rho     = rho_s_0 + rho_c_0 + rho_d_0 # mixture's density

K_Reuss    =   1/(s_c_0/K_c + s_d_0/K_d)

# density matrix terms (frequency-independent parts)
rho_11  = (rho_s_0-(1-alpha_infty)*(rho_c_0+rho_d_0))
rho_12  = (1-alpha_infty)*(rho_c_0+rho_d_0)
rho_22  = (rho_c_0-(1-alpha_infty)*(rho_c_0+rho_d_0))

# damping matrix terms (frequency-independent parts)
b_0 = eta_cR * (s_c_0*phi_0)**2/k_s # viscous damping factor (for continuous fluid-phase)
om_Biot = b_0/rho_22 # Biot's characteristic frequency for continuous phase

# stiffness matrix terms (frequency-independent parts)
N = G
Q = (1-phi_0)*K_Reuss
R = phi_0*K_Reuss
P = K + 4/3*mu + K_Reuss*(1-phi_0)**2/phi_0

# eigenfrequency and damping of a freely oscillating gas bubble
#   of course, alternative eigenfrequencies and damping can be used for other oscillators
#   damping coefficient D determines if peaks are visible in 1/Q (D < 1 yields peaks with capillary oscillation dominating, otherwise convergence to quasi-continuous behaviour with viscous damping dominating); can be sensitive wrt to radius, material parameters at ambient/reservoir conditions, etc
#   see [1] for using an integrated statistical distribution instead of individual oscillators (will later affect the variable oscillator_influence below by an integration instead of a summation)
def omega0_distr(r, gamma_ai, p0, rho_cR_0, sigma_surf):
    return np.sqrt(3 * gamma_ai * p0 / rho_cR_0 + 2 * sigma_surf * (3 * gamma_ai - 1) / (rho_cR_0 * r)) / r
def D_distr(r, eta_cR, rho_dR_0, eta_dR, gamma_ai, p0, rho_cR_0, sigma_surf):
    return (3 / 2) * (eta_cR / rho_dR_0) * ((1 + 1.5 * eta_dR / eta_cR) / (1 + eta_dR / eta_cR)) * (1 / omega0_distr(r, gamma_ai, p0, rho_cR_0, sigma_surf)) * (1 / r**2)
#inertial coupling due to form-locking between bubble and sourrounding liquid
alpha_bubble = 3.0
om_osc = [omega0_distr(r, gamma_ai, p0, rho_cR_0, sigma_surf) for r in r_osc]
D_osc = [D_distr(r, eta_cR, rho_dR_0, eta_dR, gamma_ai, p0, rho_cR_0, sigma_surf) for r in r_osc]
# -------------------------------------------------------



# -------------------------------------------------------
# Frequency range
# -------------------------------------------------------
print(" ... defining frequency range")

freq_exp_min = 2 # exponent for lower boundary of frequency range
freq_exp_max = 8 # exponent for upper boundary of frequency range
freq_exp_delta = 0.001 # log discretization of frequency range
omega_range = (10 ** np.arange(freq_exp_min, freq_exp_max + freq_exp_delta, freq_exp_delta)) * 2 * np.pi # frequency range
# -------------------------------------------------------



# -------------------------------------------------------
# Calculating dispersion relations - setup and solve eigenvalue problem analytically
# -------------------------------------------------------
print(" ... calculating dispersion relations in frequency range")

cp_1 = []
cp_2 = []
Qp_1 = []
Qp_2 = []
cs = []
Qs = []

for om in omega_range:
    
    oscillator_influence = 0
    for i in range(len(r_osc)):
        oscillator_influence = oscillator_influence + alpha_osc[i]*(om_osc[i]**2+2.0j*om*om_osc[i]*D_osc[i]) + alpha_osc[i]*((om_osc[i]**2+2.0j*om*om_osc[i]*D_osc[i])**2)/(om**2-(om_osc[i]**2+2.0j*om*om_osc[i]*D_osc[i]))

    F_Biot = np.sqrt(1+0.5j*om/om_Biot)# correction of frequency-dependent flow profile

    # matrix containing inertia and damping terms (matrix with stiffness terms described by Q, P, R and S)
    a_s = (rho_11)*(om)**2 - 1.0j*b_0*om*F_Biot # matrix coefficient A11 for inertia and damping
    a_0 = (rho_12)*(om)**2 + 1.0j*b_0*om*F_Biot # matrix coefficient A12 for inertia and damping
    a_c = (rho_22-n_d_0*(rho_dR_0+alpha_bubble*rho_cR_0))*(om)**2 - 1.0j*b_0*om*F_Biot - oscillator_influence # matrix coefficient A12 for inertia and damping
    
    # P-wave: solution of longitudinal mode
    c_1_quadr = 1/2 * (2 * a_0* Q - P * a_c - R * a_s) / (P*R - Q**2)
    c_2_quadr = (a_s* a_c- a_0**2) / (P*R - Q**2)
    Xsi_1 = (-c_1_quadr + np.sqrt(c_1_quadr**2 - c_2_quadr))
    Xsi_2 = (-c_1_quadr - np.sqrt(c_1_quadr**2 - c_2_quadr))
    
    k_analytical = np.sort(np.sqrt([Xsi_1, Xsi_2])) # wave vector = sqrt(Xsi), analytically solved - 2 solutions

    Re_k_1 = np.real(k_analytical[0])
    Re_k_2 = np.real(k_analytical[1])
    Im_k_1 = np.imag(k_analytical[0])
    Im_k_2 = np.imag(k_analytical[1])
    
    cp_1.append( om / Re_k_1 )
    cp_2.append( om / Re_k_2 )
    Qp_1.append( abs(Re_k_1 / Im_k_1/2.) )
    Qp_2.append( abs(Re_k_2 / Im_k_2/2.) )
    
    # S-wave: solution of transversal/shear mode
    Xsi_1_trans = (a_s* a_c- a_0**2) / (N * a_c)

    k_analytical_trans = np.sqrt(Xsi_1_trans) # wave vector = sqrt(Xsi), analytically solved - 1 solution

    Re_k_s  = np.real(k_analytical_trans)
    Im_k_s  = np.imag(k_analytical_trans)
    
    cs.append( om / Re_k_s )
    Qs.append( abs(Re_k_s / Im_k_s/2.) )

# end of frequency loop
# -------------------------------------------------------



# -------------------------------------------------------
# Plotting
# -------------------------------------------------------
print(" ... plotting")

fig, axs = plt.subplots(2, 3, figsize=(12, 6))

axs[0, 0].plot(omega_range,cp_1)
axs[0, 0].set(ylabel ='$c$ [m/s]', xscale = "log")
axs[0, 0].set_title('P1-wave')

axs[0, 1].plot(omega_range,cp_2, color = "gray")
axs[0, 1].set(xscale = "log")
axs[0, 1].set_title('P2-wave')

axs[0, 2].plot(omega_range,cs, color = "gray")
axs[0, 2].set(xscale = "log")
axs[0, 2].set_title('S-wave')

axs[1, 0].plot(omega_range,1./np.array(Qp_1))
axs[1, 0].set(xlabel = '$\omega$ [1/s]', ylabel ='$Q^{-1} [-]$', xscale = "log", yscale = "log")

axs[1, 1].plot(omega_range,1./np.array(Qp_2), color = "gray")
axs[1, 1].set(xlabel = '$\omega$ [1/s]', xscale = "log", yscale = "log")

axs[1, 2].plot(omega_range,1./np.array(Qs), color = "gray")
axs[1, 2].set(xlabel = '$\omega$ [1/s]', xscale = "log", yscale = "log")


plt.show()

# print(" ... and saving plot")
# fig.savefig("disp_rel.pdf", bbox_inches='tight')
# -------------------------------------------------------



print(" ")
print("*******************************************************")
print("*** Finish: Calculating dispersion relations ***")
print("*******************************************************")