Skip to content

Commit

Permalink
Computing fundamental frequency
Browse files Browse the repository at this point in the history 
Computing fundamental frequency
  • Loading branch information
@ruansava
ruansava committed Mar 4, 2025
1 parent 4ea0726 commit 2116b1a
Show file tree
Hide file tree
Showing 4 changed files with 93 additions and 91 deletions.
32 changes: 19 additions & 13 deletions spyro/habc/eik.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ def define_bcs(self, Wave):
None
'''

print('Defining Eikonal BCs')
print('\nDefining Eikonal BCs')

# Identify source locations
possou = Wave.sources.point_locations
Expand Down Expand Up @@ -354,7 +354,7 @@ def solve_eik(self, Wave, tol=1e-16, f_est=0.06):
vy = fire.TestFunction(Wave.funct_space_eik)

# Linear Eikonal
print('Solving Pre-Eikonal')
print('\nSolving Pre-Eikonal')
FeikL = self.linear_eik(Wave, u, vy)
J = fire.derivative(FeikL, yp)

Expand All @@ -378,8 +378,12 @@ def solve_eik(self, Wave, tol=1e-16, f_est=0.06):
# Solving LIN Eikonal
fire.solve(fire.lhs(FeikL) == fire.rhs(FeikL), yp,
bcs=self.bcs_eik, solver_parameters=pL, J=J)
print(
f"\nSolver Executed Successfully. AbsTol: {user_atol:.1e}")

solv_ok = "Solver Executed Successfully. "
print((solv_ok + 'AbsTol: {:.1e}').format(user_atol))

# print(
# f"\nSolver Executed Successfully. AbsTol: {user_atol:.1e}")
break

except Exception as e:
Expand All @@ -389,14 +393,14 @@ def solve_eik(self, Wave, tol=1e-16, f_est=0.06):
user_atol = user_atol * 10 if user_atol < 1e-5 \
else round(user_atol + 1e-5, 5)
if user_atol > 1e-4:
print("\nTolerance too high. Exiting.")
print("Tolerance too high. Exiting.")
break

# Clean numerical instabilities
data_eikL = self.clean_inst_num(yp.dat.data_with_halos[:])

# Nonlinear Eikonal
print('Solving Post-Eikonal')
print('\nSolving Post-Eikonal')
user_atol = tol
# vinewtonrsls, vinewtonssls, newtonls, qn, ncg, newtontr, ngs, ngmres
nl_solver = 'vinewtonssls'
Expand All @@ -421,9 +425,10 @@ def solve_eik(self, Wave, tol=1e-16, f_est=0.06):
solver_parameters=pNL, J=J)

# Final parameters
print(
f"\nSolver Executed Successfully. AbsTol: {user_atol:.1e}")
print(f"Solver Executed Successfully. Festab: {user_est:.2f}")
solv_ok = "Solver Executed Successfully. "
print((solv_ok + 'AbsTol: {:.1e}').format(user_atol))
print((solv_ok + 'Festab: {:.2f}').format(user_est))

self.yp = yp
break

Expand All @@ -438,7 +443,7 @@ def solve_eik(self, Wave, tol=1e-16, f_est=0.06):
user_atol = user_atol * 10 if user_atol < 1e-5 \
else round(user_atol + 1e-5, 5)
if user_atol > 1e-4:
print('\nHigh Tolerance. Exiting!')
print('High Tolerance. Exiting!')
break

# Save Eikonal results
Expand Down Expand Up @@ -498,15 +503,16 @@ def ident_crit_eik(self, Wave):

# Loop over boundaries
eik_bnd = []
print('\nIdentifying Critical Points on Boundaries')
for bnd, bnd_str in zip(self.bnds, bnds_str):

# Identify Eikonal minimum
eikmin, idxmin = self.ident_eik_on_bnd(bnd)
pt_cr = (self.z_data[idxmin], self.x_data[idxmin])
c_bnd = np.float64(Wave.c.at(pt_cr).item())
print('Min Eikonal on', bnd_str, '(ms):', round(1e3 * eikmin, 3),
'at (in km): (' + str(round(pt_cr[0], 3)) + ','
+ str(round(pt_cr[1], 3)) + ')')
eik_str = "Min Eikonal on {0:>16} (ms): {1:>7.3f} "
print((eik_str + 'at (in km): ({2:3.3f}, {3:3.3f})').format(
bnd_str, 1e3 * eikmin, pt_cr[0], pt_cr[1]))

# Identify closest source
lref_allsou = [np.linalg.norm(
Expand Down
81 changes: 15 additions & 66 deletions spyro/habc/habc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,10 +152,10 @@ def roundFL(self):
self.pad_len = self.F_L * lref

print('\nModifying Layer Size Based on the Element Size')
print('Modified Parameter Size F_L:', round(self.F_L, 4))
print('Modified Layer Size (km):', round(self.pad_len, 4))
print('Elements (' + str(self.lmin), 'km) in Layer:',
int(self.pad_len / self.lmin))
print("Modified Parameter Size F_L: {:.4f}".format(self.F_L))
print("Modified Layer Size (km): {:.4f}".format(self.pad_len))
print('Elements ({:.3f} km) in Layer: {}'.format(
self.lmin, int(self.pad_len / self.lmin)))

def size_habc_criterion(self, Eikonal, layer_based_on_mesh=False):
'''
Expand Down Expand Up @@ -313,21 +313,17 @@ def fundamental_frequency(self):
a_ptr, a_ind, a_val = A.petscmat.getValuesCSR()
Asp = ss.csr_matrix((a_val, a_ind, a_ptr), A.petscmat.size)

# Lsp0 = ss.linalg.eigs(Asp, k=2, M=Msp, which='SM', ncv=12,
# return_eigenvectors=False, Minv=None)

# Lsp1 = ss.linalg.eigs(Asp, k=2, M=Msp, which='SM', ncv=12,
# return_eigenvectors=False, Minv=Msp_inv)

ipdb.set_trace()
# Operator
Lsp = Asp.diagonal() * Msp_inv.diagonal()
print('\nSolving Eigenvalue Problem')
Lsp = ss.linalg.eigs(Asp, k=2, M=Msp, sigma=0.0,
return_eigenvectors=False)

for eigval in np.unique(Lsp):
print(eigval, np.sqrt(abs(eigval)) / (2 * np.pi))
# for eigval in np.unique(Lsp):
# print(np.sqrt(abs(eigval)) / (2 * np.pi))

min_eigval = np.amin(np.unique(Lsp[Lsp > 0.]))
self.fundamental_freq = np.sqrt(min_eigval) / (2 * np.pi)
min_eigval = np.unique(Lsp[(Lsp > 0.) & (np.imag(Lsp) == 0.)])[1]
self.fundam_freq = np.real(np.sqrt(min_eigval) / (2 * np.pi))
print("Fundamental Frequency (Hz): {0:.5f}".format(self.fundam_freq))

def damping_layer(self):
'''
Expand All @@ -343,9 +339,9 @@ def damping_layer(self):
'''

self.fundamental_frequency()
print(self.fundamental_freq)
# 55.95494944065594 dx = 0.05
# 281.0691908560122 dx = 0.01
# print(self.fundamental_freq)
# 0.45592718619481450 dx = 0.05
# 0.47524802941560956 dx = 0.01

# Homogeneous domain
# Dirichlet m n f Neumann Sommerfeld
Expand All @@ -365,50 +361,3 @@ def damping_layer(self):
# 1.5940 0.93186 0.93195
# 1.6356 0.95159 0.95177
# 1.7080 1.04420 1.04460


# from firedrake import *
# from slepc4py import SLEPc

# # Create mesh and function space
# mesh = UnitSquareMesh(32, 32)
# V = FunctionSpace(mesh, "CG", 1)

# # Define the bilinear and linear forms for the scalar wave equation
# u = TrialFunction(V)
# v = TestFunction(V)
# c = Constant(1.0) # Wave speed
# a = c**2 * inner(grad(u), grad(v)) * dx
# m = u * v * dx

# # Assemble the matrices
# A = assemble(a)
# M = assemble(m)

# # Set up the SLEPc eigensolver
# eigensolver = SLEPc.EPS().create()
# A_petsc = A.M.handle
# M_petsc = M.M.handle
# eigensolver.setOperators(A_petsc, M_petsc)
# eigensolver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
# eigensolver.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)
# eigensolver.setFromOptions()

# # Solve the eigenvalue problem
# eigensolver.solve()

# # Extract the eigenvalues
# nconv = eigensolver.getConverged()
# eigenvalues = []
# for i in range(nconv):
# eigenvalue = eigensolver.getEigenvalue(i)
# eigenvalues.append(eigenvalue)

# print("Eigenvalues:", eigenvalues)



# https://github.com/firedrakeproject/slepc
# https://github.com/firedrakeproject/firedrake/blob/master/docs/source/install.rst
# https://slepc.upv.es/documentation/instal.htm
# https://pypi.org/project/slepc4py/
27 changes: 15 additions & 12 deletions spyro/habc/lay_len.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,24 +148,24 @@ def calc_size_lay(Wave, nz=5, crtCR=1, tol_rel=1e-3, monitor=False):
print('\nComputing Size for Absorbing Layer')
# z_par: Inverse of minimum Eikonal (Equivalent to c_bound / lref)
z_par = Wave.eik_bnd[0][3]
aux0 = f"Parameter z (1/s): {round(z_par, 4)},"
aux0 = "Parameter z (1/s): {:.4f},".format(z_par)
# fref: Reference frequency
fref = Wave.frequency
aux1 = f"Reference Frequency (Hz): {round(fref, 4)}"
aux1 = "Reference Frequency (Hz): {:.4f}".format(fref)
print(aux0, aux1)

# lmin: Minimal dimension of finite element in mesh
lmin = Wave.lmin
aux2 = f"Minimum Mesh Length (km): {lmin},"
aux2 = "Minimum Mesh Length (km): {:.4f},".format(lmin)
# lref: Reference length for the size of the absorbing layer
lref = Wave.eik_bnd[0][4]
aux3 = f"Reference Length (km): {round(lref, 4)}"
aux3 = "Reference Length (km): {:.4f}".format(lref)
print(aux2, aux3)

a = z_par / fref # Adimensional parameter
aux4 = f"Parameter a: {round(a, 4)},"
aux4 = "Parameter a: {:.4f},".format(a)
FLmin = 0.5 * lmin / lref # Initial guess
aux5 = f"Initial Guess: {round(FLmin, 4)}"
aux5 = "Initial Guess: {:.4f}".format(FLmin)
print(aux4, aux5)

x = FLmin
Expand Down Expand Up @@ -199,11 +199,14 @@ def calc_size_lay(Wave, nz=5, crtCR=1, tol_rel=1e-3, monitor=False):
pad_len = F_L * lref

# Visualizing options for layer size
print('Selected Parameter Size F_L:', round(F_L, 4))
print('Options for F_L:', [round(float(x), 4) for x in FLpos])
print('Aproximated Number of Elements (' + str(Wave.lmin),
'km) in Layer:', [int(x * lref / lmin) for x in FLpos])
print('Options for CR:', CRpos)
print('Selected Layer Size (km):', round(pad_len, 4))
print("Selected Parameter Size F_L: {:.4f}".format(F_L))
format_FL = ', '.join(['{:.4f}'.format(float(x)) for x in FLpos])
print("Options for F_L: [{}]".format(format_FL))
format_CR = ', '.join(['{:.4f}'.format(x) for x in CRpos])
print("Options for CR: [{}]".format(format_CR))
format_ele = [int(x * lref / lmin) for x in FLpos]
print("Aprox. Number of Elements ({:.3f} km) in Layer: {}".format(
lmin, format_ele))
print("Selected Layer Size (km): {:.4f}".format(pad_len))

return F_L, pad_len
44 changes: 44 additions & 0 deletions test_modal.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
from firedrake import *
from slepc4py import SLEPc

# Create mesh and function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 1)

# Define the bilinear and linear forms for the scalar wave equation
u = TrialFunction(V)
v = TestFunction(V)
c = Constant(1.0) # Wave speed
a = c**2 * inner(grad(u), grad(v)) * dx
m = u * v * dx

# Assemble the matrices
A = assemble(a)
M = assemble(m)

# Set up the SLEPc eigensolver
eigensolver = SLEPc.EPS().create()
A_petsc = A.M.handle
M_petsc = M.M.handle
eigensolver.setOperators(A_petsc, M_petsc)
eigensolver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eigensolver.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)
eigensolver.setFromOptions()

# Solve the eigenvalue problem
eigensolver.solve()

# Extract the eigenvalues
nconv = eigensolver.getConverged()
eigenvalues = []
for i in range(nconv):
eigenvalue = eigensolver.getEigenvalue(i)
eigenvalues.append(eigenvalue)

print("Eigenvalues:", eigenvalues)


https://github.com/firedrakeproject/slepc
https://github.com/firedrakeproject/firedrake/blob/master/docs/source/install.rst
https://slepc.upv.es/documentation/instal.htm
https://pypi.org/project/slepc4py/

0 comments on commit 2116b1a

Please sign in to comment.