From 2116b1a8fbcc30bdd89e95a46e92da8ee81eaec8 Mon Sep 17 00:00:00 2001 From: ruansava Date: Tue, 4 Mar 2025 20:26:34 -0300 Subject: [PATCH] Computing fundamental frequency Computing fundamental frequency --- spyro/habc/eik.py | 32 ++++++++++------- spyro/habc/habc.py | 81 ++++++++----------------------------------- spyro/habc/lay_len.py | 27 ++++++++------- test_modal.py | 44 +++++++++++++++++++++++ 4 files changed, 93 insertions(+), 91 deletions(-) create mode 100644 test_modal.py diff --git a/spyro/habc/eik.py b/spyro/habc/eik.py index 71d0a747..6c604530 100644 --- a/spyro/habc/eik.py +++ b/spyro/habc/eik.py @@ -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 @@ -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) @@ -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: @@ -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' @@ -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 @@ -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 @@ -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( diff --git a/spyro/habc/habc.py b/spyro/habc/habc.py index 0e91c489..e76898ad 100644 --- a/spyro/habc/habc.py +++ b/spyro/habc/habc.py @@ -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): ''' @@ -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): ''' @@ -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 @@ -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/ diff --git a/spyro/habc/lay_len.py b/spyro/habc/lay_len.py index 5be16e11..49576c36 100644 --- a/spyro/habc/lay_len.py +++ b/spyro/habc/lay_len.py @@ -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 @@ -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 diff --git a/test_modal.py b/test_modal.py new file mode 100644 index 00000000..3df2a5f9 --- /dev/null +++ b/test_modal.py @@ -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/ \ No newline at end of file