diff --git a/analyze_temporal_convergence.py b/analyze_temporal_convergence.py new file mode 100644 index 000000000..a0800045e --- /dev/null +++ b/analyze_temporal_convergence.py @@ -0,0 +1,205 @@ +""" +Analyze temporal convergence of MMS results. +For a central difference scheme in time, we expect second-order convergence (rate = 2). +This means error should decrease by ~4x when dt is halved. +""" + +import numpy as np +import re +import sys +from collections import defaultdict + + +class TeeOutput: + """Write to both console and file simultaneously.""" + def __init__(self, filename): + self.file = open(filename, 'w') + self.stdout = sys.stdout + + def write(self, message): + self.stdout.write(message) + self.file.write(message) + + def flush(self): + self.stdout.flush() + self.file.flush() + + def close(self): + self.file.close() + + +def parse_mms_results(filename): + """Parse the MMS results file and organize data by h value.""" + data = defaultdict(lambda: {"dt": [], "e1": [], "e2": []}) + + with open(filename, 'r') as f: + lines = f.readlines() + + current_dt = None + i = 0 + while i < len(lines): + line = lines[i].strip() + + # Look for dt value + if line.startswith("dt = "): + dt_str = line.split("=")[1].strip() + current_dt = float(dt_str) + i += 1 + continue + + # Skip separator lines and headers + if line.startswith("-") or line.startswith("h ") or not line or line.startswith("="): + i += 1 + continue + + # Parse data lines (h, e1, e2) + if current_dt is not None: + parts = line.split() + if len(parts) >= 3 and parts[1] != "Failed" and parts[2] != "Failed": + h = float(parts[0]) + e1 = float(parts[1]) + e2 = float(parts[2]) + + data[h]["dt"].append(current_dt) + data[h]["e1"].append(e1) + data[h]["e2"].append(e2) + + i += 1 + + return data + + +def calculate_convergence_rate(dt_values, errors): + """ + Calculate convergence rate between consecutive refinements. + For second-order convergence: error ~ dt^p, where p should be ~2 + Rate = log(error1/error2) / log(dt1/dt2) + """ + rates = [] + for i in range(len(dt_values) - 1): + if errors[i] > 0 and errors[i+1] > 0: + rate = np.log(errors[i] / errors[i+1]) / np.log(dt_values[i] / dt_values[i+1]) + rates.append(rate) + else: + rates.append(np.nan) + return rates + + +def analyze_convergence(data, expected_rate=2.0, tolerance=0.1): + """ + Analyze convergence for each h value. + + Parameters: + ----------- + data : dict + Dictionary with h values as keys + expected_rate : float + Expected convergence rate (2.0 for second-order) + tolerance : float + Acceptable deviation from expected rate + """ + print("=" * 80) + print("TEMPORAL CONVERGENCE ANALYSIS") + print("Expected convergence rate: {:.1f} (second-order in time)".format(expected_rate)) + print("Tolerance: ± {:.1f}".format(tolerance)) + print("=" * 80) + print() + + h_values = sorted(data.keys()) + + for h in h_values: + print(f"h = {h}") + print("-" * 80) + + # Sort by dt (largest to smallest) + dt_vals = np.array(data[h]["dt"]) + e1_vals = np.array(data[h]["e1"]) + e2_vals = np.array(data[h]["e2"]) + + sorted_idx = np.argsort(dt_vals)[::-1] + dt_sorted = dt_vals[sorted_idx] + e1_sorted = e1_vals[sorted_idx] + e2_sorted = e2_vals[sorted_idx] + + # Calculate convergence rates + rates_e1 = calculate_convergence_rate(dt_sorted, e1_sorted) + rates_e2 = calculate_convergence_rate(dt_sorted, e2_sorted) + + # Print table + print(f"{'dt':<12} {'e1':<15} {'rate(e1)':<12} {'e2':<15} {'rate(e2)':<12} {'Status':<15}") + print("-" * 80) + + for i in range(len(dt_sorted)): + dt_str = f"{dt_sorted[i]:.5g}" + e1_str = f"{e1_sorted[i]:.6e}" + e2_str = f"{e2_sorted[i]:.6e}" + + if i < len(rates_e1): + rate_e1_str = f"{rates_e1[i]:.2f}" + rate_e2_str = f"{rates_e2[i]:.2f}" + + # Check if converging at expected rate + e1_converging = abs(rates_e1[i] - expected_rate) <= tolerance + e2_converging = abs(rates_e2[i] - expected_rate) <= tolerance + + if e1_converging and e2_converging: + status = "✓ PASS" + elif e1_converging or e2_converging: + status = "~ PARTIAL" + else: + status = "✗ FAIL" + else: + rate_e1_str = "-" + rate_e2_str = "-" + status = "-" + + print(f"{dt_str:<12} {e1_str:<15} {rate_e1_str:<12} {e2_str:<15} {rate_e2_str:<12} {status:<15}") + + # Summary for this h + if len(rates_e1) > 0: + valid_rates_e1 = [r for r in rates_e1 if not np.isnan(r)] + valid_rates_e2 = [r for r in rates_e2 if not np.isnan(r)] + + if valid_rates_e1: + avg_rate_e1 = np.mean(valid_rates_e1) + avg_rate_e2 = np.mean(valid_rates_e2) + + print() + print(f"Average convergence rate (e1): {avg_rate_e1:.2f}") + print(f"Average convergence rate (e2): {avg_rate_e2:.2f}") + + e1_ok = abs(avg_rate_e1 - expected_rate) <= tolerance + e2_ok = abs(avg_rate_e2 - expected_rate) <= tolerance + + if e1_ok and e2_ok: + print(f"Overall for h={h}: ✓ CONVERGING at expected second-order rate") + else: + print(f"Overall for h={h}: ✗ NOT converging at expected rate") + + print() + print() + + +if __name__ == "__main__": + input_filename = "mms_results.txt" + output_filename = "temporal_convergence_analysis.txt" + + # Parse data + data = parse_mms_results(input_filename) + + if not data: + print(f"Error: No valid data found in {input_filename}") + else: + # Redirect output to both console and file + tee = TeeOutput(output_filename) + sys.stdout = tee + + try: + # Analyze convergence + analyze_convergence(data, expected_rate=2.0, tolerance=0.1) + finally: + # Restore stdout and close file + sys.stdout = tee.stdout + tee.close() + + print(f"\nAnalysis complete. Results saved to {output_filename}") diff --git a/check_MMS.py b/check_MMS.py new file mode 100644 index 000000000..1acf2d189 --- /dev/null +++ b/check_MMS.py @@ -0,0 +1,205 @@ +import firedrake as fire +import math +import spyro +import matplotlib.pyplot as plt +from collections import defaultdict +import csv + +def mms_h_convergence(dt, h): + u1 = lambda x, t: (x[0]**2 + x[0])*(x[1]**2 - x[1])*t + u2 = lambda x, t: (2*x[0]**2 + 2*x[0])*(-x[1]**2 + x[1])*t + u = lambda x, t: fire.as_vector([u1(x, t), u2(x, t)]) + + b1 = lambda x, t: -(2*x[0]**2 + 6*x[1]**2 - 16*x[0]*x[1] + 10*x[0] - 14*x[1] + 4)*t + b2 = lambda x, t: -(-12*x[0]**2 - 4*x[1]**2 + 8*x[0]*x[1] - 16*x[0] + 8*x[1] - 2)*t + b = lambda x, t: fire.as_vector([b1(x, t), b2(x, t)]) + + fo = int(0.1/dt) + + dictionary = {} + dictionary["options"] = { + "cell_type": "Q", # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q) + "variant": "lumped", # lumped, equispaced or DG, default is lumped "method":"MLT", # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral) You can either specify a cell_type+variant or a method + "degree": 4, # p order + "dimension": 2, # dimension + } + dictionary["parallelism"] = { + "type": "automatic", # options: automatic (same number of cores for evey processor) or spatial + } + dictionary["mesh"] = { + "Lz": 1.0, # depth in km - always positive + "Lx": 1.0, # width in km - always positive + "Ly": 0.0, # thickness in km - always positive + "mesh_type": "firedrake_mesh", # options: firedrake_mesh or user_mesh + "mesh_file": None, # specify the mesh file + } + dictionary["acquisition"] = { + "source_type": "MMS", + "source_locations": [(-1.0, 1.0)], + "frequency": 5.0, + "delay": 1.5, + "receiver_locations": [(-0.0, 0.5)], + "body_forces": b, + } + dictionary["time_axis"] = { + "initial_condition": u, + "initial_time": 0.0, # Initial time for event + "final_time": 1.0, # Final time for event + "dt": dt, # timestep size + "amplitude": 1, # the Ricker has an amplitude of 1. + "output_frequency": fo, # how frequently to output solution to pvds + "gradient_sampling_frequency": 1, # how frequently to save solution to RAM + } + + dictionary["visualization"] = { + "forward_output": False, + "output_filename": "results/forward_output.pvd", + "fwi_velocity_model_output": False, + "velocity_model_filename": None, + "gradient_output": False, + "gradient_filename": None, + } + + dictionary["synthetic_data"] = { + "type": "object", + "density": 1, + "lambda": 1, + "mu": 1, + "real_velocity_file": None, + } + dictionary["boundary_conditions"] = [ + ("uz", "on_boundary", 0), + ("ux", "on_boundary", 0), + ] + + wave = spyro.IsotropicWave(dictionary) + wave.set_mesh(input_mesh_parameters={"edge_length": h}) + # Lets build de matrix operator outside of the forward solve so we can + # set previous timesteps for the MMS problem + wave._initialize_model_parameters() + wave.matrix_building() + x = wave.get_spatial_coordinates() + # wave.u_nm1.interpolate(u(x, 0.0 - 2*dt)) + # wave.u_n.interpolate(u(x, 0.0 - dt)) + # wave.u_nm2.interpolate(u(x, 0.0 - 2*dt)) + + wave.forward_solve(build_matrix_operator=False) + + u_an = fire.Function(wave.function_space) + t = wave.current_time + u_an.interpolate(u(x, t)) + + e1 = fire.errornorm(wave.u_n.sub(0), u_an.sub(0)) / fire.norm(u_an.sub(0)) + e2 = fire.errornorm(wave.u_n.sub(1), u_an.sub(1)) / fire.norm(u_an.sub(1)) + + if e1 > 1e3 or e2 > 1e3: + raise ValueError("ERROR") + + print(f"For dt: {dt}, h: {h}, e1 = {e1}, e2 = {e2}") + return e1, e2 + +if __name__ == "__main__": + dts = [1e-3, 1e-4, 1e-5] + # dts = [1e-2, 1e-3] + hs = [0.125]#, 0.1, 0.05, 0.025] + + # Dictionary to store results: results[dt][h] = (e1, e2) or "Error" + results = defaultdict(dict) + + for dt in dts: + for h in hs: + # try: + e1, e2 = mms_h_convergence(dt, h) + results[dt][h] = (e1, e2) + # except Exception as e: + # print(f"Error at dt: {dt}, h: {h}: {str(e)}") + # results[dt][h] = "Error" + + # Print results in table format + print("\n" + "="*80) + print("RESULTS TABLE") + print("="*80) + + for dt in dts: + print(f"\ndt = {dt}") + print("-" * 60) + print(f"{'h':<10} {'e1':<25} {'e2':<25}") + print("-" * 60) + for h in hs: + if h in results[dt]: + if results[dt][h] == "Error": + print(f"{h:<10} {'Failed':<25} {'Failed':<25}") + else: + e1, e2 = results[dt][h] + print(f"{h:<10} {e1:<25.6e} {e2:<25.6e}") + print() + + # Save results to CSV file + print("\nSaving results to CSV file...") + with open('mms_results.csv', 'w', newline='') as csvfile: + csvwriter = csv.writer(csvfile) + csvwriter.writerow(['dt', 'h', 'e1', 'e2', 'status']) + + for dt in dts: + for h in sorted(results[dt].keys()): + if results[dt][h] == "Error": + csvwriter.writerow([dt, h, 'Error', 'Error', 'Failed']) + else: + e1, e2 = results[dt][h] + csvwriter.writerow([dt, h, f'{e1:.6e}', f'{e2:.6e}', 'Success']) + + print("Saved: mms_results.csv") + + # Save results to TXT file + print("Saving results to TXT file...") + with open('mms_results.txt', 'w') as txtfile: + txtfile.write("="*80 + "\n") + txtfile.write("MMS CONVERGENCE TEST RESULTS\n") + txtfile.write("="*80 + "\n\n") + + for dt in dts: + txtfile.write(f"dt = {dt}\n") + txtfile.write("-" * 60 + "\n") + txtfile.write(f"{'h':<10} {'e1':<25} {'e2':<25}\n") + txtfile.write("-" * 60 + "\n") + for h in sorted(results[dt].keys()): + if results[dt][h] == "Error": + txtfile.write(f"{h:<10} {'Failed':<25} {'Failed':<25}\n") + else: + e1, e2 = results[dt][h] + txtfile.write(f"{h:<10} {e1:<25.6e} {e2:<25.6e}\n") + txtfile.write("\n") + + print("Saved: mms_results.txt") + + # Create plots for each dt + print("\nGenerating convergence plots...") + + for dt in dts: + # Filter out errors and get successful runs + successful_h = [] + successful_e1 = [] + successful_e2 = [] + + for h in sorted(results[dt].keys()): + if results[dt][h] != "Error": + e1, e2 = results[dt][h] + successful_h.append(h) + successful_e1.append(e1) + successful_e2.append(e2) + + if successful_h: # Only plot if we have data + plt.figure(figsize=(10, 6)) + plt.loglog(successful_h, successful_e1, 'o-', label='e1 error', linewidth=2, markersize=8) + plt.loglog(successful_h, successful_e2, 's-', label='e2 error', linewidth=2, markersize=8) + plt.xlabel('Mesh size h', fontsize=12) + plt.ylabel('Error', fontsize=12) + plt.title(f'MMS Convergence for dt = {dt}', fontsize=14) + plt.legend(fontsize=11) + plt.grid(True, which="both", ls="-", alpha=0.3) + plt.tight_layout() + plt.savefig(f'mms_convergence_dt_{dt}.png', dpi=150) + print(f"Saved plot: mms_convergence_dt_{dt}.png") + plt.close() + + print("\nAll plots generated successfully!") diff --git a/check_MMS_elastic_as_acoustic.py b/check_MMS_elastic_as_acoustic.py new file mode 100644 index 000000000..4bab953eb --- /dev/null +++ b/check_MMS_elastic_as_acoustic.py @@ -0,0 +1,209 @@ +import firedrake as fire +import math +import spyro +import matplotlib.pyplot as plt +from collections import defaultdict +import csv + +def mms_h_convergence(dt, h): + # Acoustic-like solution: both components have the same scalar solution + scalar_solution = lambda x, t: x[0] * (x[0] + 1) * x[1] * (x[1] - 1) * t + u1 = lambda x, t: scalar_solution(x, t) + u2 = lambda x, t: scalar_solution(x, t) + u = lambda x, t: fire.as_vector([u1(x, t), u2(x, t)]) + + # Body forces for acoustic-like case (both components get the same source) + scalar_source = lambda x, t: -(x[0]**2 + x[0] + x[1]**2 - x[1]) * 2 * t + b1 = lambda x, t: scalar_source(x, t) + b2 = lambda x, t: scalar_source(x, t) + b = lambda x, t: fire.as_vector([b1(x, t), b2(x, t)]) + + fo = int(0.1/dt) + + dictionary = {} + dictionary["options"] = { + "cell_type": "Q", # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q) + "variant": "lumped", # lumped, equispaced or DG, default is lumped "method":"MLT", # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral) You can either specify a cell_type+variant or a method + "degree": 4, # p order + "dimension": 2, # dimension + } + dictionary["parallelism"] = { + "type": "automatic", # options: automatic (same number of cores for evey processor) or spatial + } + dictionary["mesh"] = { + "Lz": 1.0, # depth in km - always positive + "Lx": 1.0, # width in km - always positive + "Ly": 0.0, # thickness in km - always positive + "mesh_type": "firedrake_mesh", # options: firedrake_mesh or user_mesh + "mesh_file": None, # specify the mesh file + } + dictionary["acquisition"] = { + "source_type": "MMS", + "source_locations": [(-1.0, 1.0)], + "frequency": 5.0, + "delay": 1.5, + "receiver_locations": [(-0.0, 0.5)], + "body_forces": b, + } + dictionary["time_axis"] = { + "initial_condition": u, + "initial_time": 0.0, # Initial time for event + "final_time": 1.0, # Final time for event + "dt": dt, # timestep size + "amplitude": 1, # the Ricker has an amplitude of 1. + "output_frequency": fo, # how frequently to output solution to pvds + "gradient_sampling_frequency": 1, # how frequently to save solution to RAM + } + + dictionary["visualization"] = { + "forward_output": False, + "output_filename": "results/forward_output.pvd", + "fwi_velocity_model_output": False, + "velocity_model_filename": None, + "gradient_output": False, + "gradient_filename": None, + } + + dictionary["synthetic_data"] = { + "type": "object", + "density": 1, + "lambda": 1, + "mu": 0, # Set to 0 to make elastic equivalent to acoustic + "real_velocity_file": None, + } + dictionary["boundary_conditions"] = [ + ("uz", "on_boundary", 0), + ("ux", "on_boundary", 0), + ] + + wave = spyro.IsotropicWave(dictionary) + wave.set_mesh(input_mesh_parameters={"edge_length": h}) + # Lets build de matrix operator outside of the forward solve so we can + # set previous timesteps for the MMS problem + wave._initialize_model_parameters() + wave.matrix_building() + x = wave.get_spatial_coordinates() + # wave.u_nm1.interpolate(u(x, 0.0 - 2*dt)) + # wave.u_n.interpolate(u(x, 0.0 - dt)) + # wave.u_nm2.interpolate(u(x, 0.0 - 2*dt)) + + wave.forward_solve(build_matrix_operator=False) + + u_an = fire.Function(wave.function_space) + t = wave.current_time + u_an.interpolate(u(x, t)) + + e1 = fire.errornorm(wave.u_n.sub(0), u_an.sub(0)) / fire.norm(u_an.sub(0)) + e2 = fire.errornorm(wave.u_n.sub(1), u_an.sub(1)) / fire.norm(u_an.sub(1)) + + if e1 > 1e3 or e2 > 1e3: + raise ValueError("ERROR") + + print(f"For dt: {dt}, h: {h}, e1 = {e1}, e2 = {e2}") + return e1, e2 + +if __name__ == "__main__": + dts = [1e-3, 1e-4, 1e-5] + # dts = [1e-2, 1e-3] + hs = [0.125]#, 0.1, 0.05, 0.025] + + # Dictionary to store results: results[dt][h] = (e1, e2) or "Error" + results = defaultdict(dict) + + for dt in dts: + for h in hs: + # try: + e1, e2 = mms_h_convergence(dt, h) + results[dt][h] = (e1, e2) + # except Exception as e: + # print(f"Error at dt: {dt}, h: {h}: {str(e)}") + # results[dt][h] = "Error" + + # Print results in table format + print("\n" + "="*80) + print("RESULTS TABLE") + print("="*80) + + for dt in dts: + print(f"\ndt = {dt}") + print("-" * 60) + print(f"{'h':<10} {'e1':<25} {'e2':<25}") + print("-" * 60) + for h in hs: + if h in results[dt]: + if results[dt][h] == "Error": + print(f"{h:<10} {'Failed':<25} {'Failed':<25}") + else: + e1, e2 = results[dt][h] + print(f"{h:<10} {e1:<25.6e} {e2:<25.6e}") + print() + + # Save results to CSV file + print("\nSaving results to CSV file...") + with open('mms_results.csv', 'w', newline='') as csvfile: + csvwriter = csv.writer(csvfile) + csvwriter.writerow(['dt', 'h', 'e1', 'e2', 'status']) + + for dt in dts: + for h in sorted(results[dt].keys()): + if results[dt][h] == "Error": + csvwriter.writerow([dt, h, 'Error', 'Error', 'Failed']) + else: + e1, e2 = results[dt][h] + csvwriter.writerow([dt, h, f'{e1:.6e}', f'{e2:.6e}', 'Success']) + + print("Saved: mms_results.csv") + + # Save results to TXT file + print("Saving results to TXT file...") + with open('mms_results.txt', 'w') as txtfile: + txtfile.write("="*80 + "\n") + txtfile.write("MMS CONVERGENCE TEST RESULTS\n") + txtfile.write("="*80 + "\n\n") + + for dt in dts: + txtfile.write(f"dt = {dt}\n") + txtfile.write("-" * 60 + "\n") + txtfile.write(f"{'h':<10} {'e1':<25} {'e2':<25}\n") + txtfile.write("-" * 60 + "\n") + for h in sorted(results[dt].keys()): + if results[dt][h] == "Error": + txtfile.write(f"{h:<10} {'Failed':<25} {'Failed':<25}\n") + else: + e1, e2 = results[dt][h] + txtfile.write(f"{h:<10} {e1:<25.6e} {e2:<25.6e}\n") + txtfile.write("\n") + + print("Saved: mms_results.txt") + + # Create plots for each dt + print("\nGenerating convergence plots...") + + for dt in dts: + # Filter out errors and get successful runs + successful_h = [] + successful_e1 = [] + successful_e2 = [] + + for h in sorted(results[dt].keys()): + if results[dt][h] != "Error": + e1, e2 = results[dt][h] + successful_h.append(h) + successful_e1.append(e1) + successful_e2.append(e2) + + if successful_h: # Only plot if we have data + plt.figure(figsize=(10, 6)) + plt.loglog(successful_h, successful_e1, 'o-', label='e1 error', linewidth=2, markersize=8) + plt.loglog(successful_h, successful_e2, 's-', label='e2 error', linewidth=2, markersize=8) + plt.xlabel('Mesh size h', fontsize=12) + plt.ylabel('Error', fontsize=12) + plt.title(f'MMS Convergence for dt = {dt}', fontsize=14) + plt.legend(fontsize=11) + plt.grid(True, which="both", ls="-", alpha=0.3) + plt.tight_layout() + plt.savefig(f'mms_convergence_dt_{dt}.png', dpi=150) + print(f"Saved plot: mms_convergence_dt_{dt}.png") + plt.close() + + print("\nAll plots generated successfully!") diff --git a/check_acousticMMS_inElastic.py b/check_acousticMMS_inElastic.py new file mode 100644 index 000000000..62c3a0044 --- /dev/null +++ b/check_acousticMMS_inElastic.py @@ -0,0 +1,276 @@ +import firedrake as fire +import math +import numpy as np +import spyro +import matplotlib.pyplot as plt +from collections import defaultdict +import csv + + +def mms_h_convergence(dt, h): + c = lambda x: 1 + fire.sin(fire.pi*-x[0])*fire.sin(fire.pi*x[1]) + B = lambda x: x[0]*(x[0] + 1.0)*x[1]*(x[1] - 1.0) + Q = lambda x: 1.0 / (c(x)**2) + + dQ_dx0 = lambda x: (2.0 * fire.pi * fire.cos(fire.pi * x[0]) * fire.sin(fire.pi * x[1])) / (c(x)**3) + dQ_dx1 = lambda x: (2.0 * fire.pi * fire.sin(fire.pi * x[0]) * fire.cos(fire.pi * x[1])) / (c(x)**3) + + u1 = lambda x, t: -t * ( + (2.0*x[0] + 1.0) * x[1] * (x[1] - 1.0) * Q(x) + + B(x) * dQ_dx0(x) + ) + + u2 = lambda x, t: -t * ( + x[0] * (x[0] + 1.0) * (2.0*x[1] - 1.0) * Q(x) + + B(x) * dQ_dx1(x) + ) + + u = lambda x, t: fire.as_vector((u1(x, t), u2(x, t))) + + # Corresponding MMS pressure (for verification) + p = lambda x, t: x[0]*(x[0]+1)*x[1]*(x[1]-1)*t + + # Elastic MMS body force (μ = 0): + b1 = lambda x, t: (2*x[0] + 1) * x[1] * (x[1] - 1) * t + b2 = lambda x, t: (2*x[1] - 1) * x[0] * (x[0] + 1) * t + b = lambda x, t: fire.as_vector([b1(x, t), b2(x, t)]) + + fo = int(0.1/dt) + + lbda = lambda x: c(x)*c(x) + + dictionary = {} + dictionary["options"] = { + "cell_type": "Q", # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q) + "variant": "lumped", # lumped, equispaced or DG, default is lumped "method":"MLT", # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral) You can either specify a cell_type+variant or a method + "degree": 4, # p order + "dimension": 2, # dimension + } + dictionary["parallelism"] = { + "type": "automatic", # options: automatic (same number of cores for evey processor) or spatial + } + dictionary["mesh"] = { + "Lz": 1.0, # depth in km - always positive + "Lx": 1.0, # width in km - always positive + "Ly": 0.0, # thickness in km - always positive + "mesh_type": "firedrake_mesh", # options: firedrake_mesh or user_mesh + "mesh_file": None, # specify the mesh file + } + dictionary["acquisition"] = { + "source_type": "MMS", + "source_locations": [(-1.0, 1.0)], + "frequency": 5.0, + "delay": 1.5, + "receiver_locations": [(-0.0, 0.5)], + "body_forces": b, + } + dictionary["time_axis"] = { + "initial_condition": u, + "initial_time": 0.0, # Initial time for event + "final_time": 1.0, # Final time for event + "dt": dt, # timestep size + "amplitude": 1, # the Ricker has an amplitude of 1. + "output_frequency": fo, # how frequently to output solution to pvds + "gradient_sampling_frequency": 1, # how frequently to save solution to RAM + } + + dictionary["visualization"] = { + "forward_output": False, + "output_filename": "results/forward_output.pvd", + "fwi_velocity_model_output": False, + "velocity_model_filename": None, + "gradient_output": False, + "gradient_filename": None, + } + + dictionary["synthetic_data"] = { + "type": "object", + "density": 1, + "mu": 1, + "real_velocity_file": None, + } + dictionary["boundary_conditions"] = [ + ("uz", "on_boundary", 0), + ("ux", "on_boundary", 0), + ] + + wave = spyro.IsotropicWave(dictionary) + wave.set_mesh(input_mesh_parameters={"edge_length": h}) + # Lets build de matrix operator outside of the forward solve so we can + # set previous timesteps for the MMS problem + # wave._initialize_model_parameters() + # wave.matrix_building() + x = wave.get_spatial_coordinates() + dictionary["synthetic_data"]["lambda"] = lbda(x) + # wave.u_nm1.interpolate(u(x, 0.0 - 2*dt)) + # wave.u_n.interpolate(u(x, 0.0 - dt)) + # wave.u_nm2.interpolate(u(x, 0.0 - 2*dt)) + + wave.forward_solve() + + u_an = fire.Function(wave.function_space) + t = wave.current_time + u_an.interpolate(u(x, t)) + + V = wave.function_space.sub(0) + + # For elastic waves: ∂p/∂t + λ ∇·u = 0 + # The correct relationship is: p = -λ * ∇·u + p_num = fire.project(-lbda(x) * fire.div(wave.u_n), V) + p_numdat = p_num.dat.data[:] + spyro.plots.debug_pvd(p_num, filename="numerical.pvd") + + p_an = fire.Function(wave.function_space.sub(0)).interpolate(p(x, t)) + p_andat = p_an.dat.data[:] + + spyro.plots.debug_pvd(p_an, filename="analitical.pvd") + + # Debug the issue with nan errors + print(f"p_num stats: min={min(p_numdat)}, max={max(p_numdat)}, mean={sum(p_numdat)/len(p_numdat)}") + print(f"p_an stats: min={min(p_andat)}, max={max(p_andat)}, mean={sum(p_andat)/len(p_andat)}") + + # Check norms + norm_p_num = fire.norm(p_num) + norm_p_an = fire.norm(p_an) + print(f"Norms: p_num={norm_p_num}, p_an={norm_p_an}") + + # Check for special values + has_nan_pnum = np.any(np.isnan(p_numdat)) + has_nan_pan = np.any(np.isnan(p_andat)) + has_inf_pnum = np.any(np.isinf(p_numdat)) + has_inf_pan = np.any(np.isinf(p_andat)) + print(f"NaN/Inf check: p_num has NaN={has_nan_pnum}, has Inf={has_inf_pnum}") + print(f"NaN/Inf check: p_an has NaN={has_nan_pan}, has Inf={has_inf_pan}") + + # Manual L2 error calculation since fire.norm sometimes gives inf + manual_norm_p_an = np.sqrt(np.sum(p_andat**2)) + manual_norm_diff = np.sqrt(np.sum((p_numdat - p_andat)**2)) + + if manual_norm_p_an > 1e-15: + err_p = manual_norm_diff / manual_norm_p_an + print(f"Pressure recovery error: {err_p:.6f}") + else: + err_p = float('inf') + print("Cannot compute pressure error - p_an norm too small") + + # Compute displacement errors + e1 = fire.errornorm(wave.u_n, u_an) / fire.norm(u_an) + print("Displacement MMS error:", e1) + + # Use pressure error as second metric + e2 = err_p + + print(f"p_num max: {max(abs(p_numdat))}, p_an max: {max(abs(p_andat))}") + print(f"p_num shape: {p_numdat.shape}, p_an shape: {p_andat.shape}") + + if e1 > 1e3 or e2 > 1e3: + raise ValueError("ERROR") + + print(f"For dt: {dt}, h: {h}, e1 = {e1}, e2 = {e2}") + return e1, e2 + +if __name__ == "__main__": + dts = [1e-3] # , 1e-4, 1e-5] + # dts = [1e-2, 1e-3] + hs = [0.125]#, 0.1, 0.05, 0.025] + + # Dictionary to store results: results[dt][h] = (e1, e2) or "Error" + results = defaultdict(dict) + + for dt in dts: + for h in hs: + # try: + e1, e2 = mms_h_convergence(dt, h) + results[dt][h] = (e1, e2) + # except Exception as e: + # print(f"Error at dt: {dt}, h: {h}: {str(e)}") + # results[dt][h] = "Error" + + # Print results in table format + print("\n" + "="*80) + print("RESULTS TABLE") + print("="*80) + + for dt in dts: + print(f"\ndt = {dt}") + print("-" * 60) + print(f"{'h':<10} {'e1':<25} {'e2':<25}") + print("-" * 60) + for h in hs: + if h in results[dt]: + if results[dt][h] == "Error": + print(f"{h:<10} {'Failed':<25} {'Failed':<25}") + else: + e1, e2 = results[dt][h] + print(f"{h:<10} {e1:<25.6e} {e2:<25.6e}") + print() + + # Save results to CSV file + print("\nSaving results to CSV file...") + with open('mms_results.csv', 'w', newline='') as csvfile: + csvwriter = csv.writer(csvfile) + csvwriter.writerow(['dt', 'h', 'e1', 'e2', 'status']) + + for dt in dts: + for h in sorted(results[dt].keys()): + if results[dt][h] == "Error": + csvwriter.writerow([dt, h, 'Error', 'Error', 'Failed']) + else: + e1, e2 = results[dt][h] + csvwriter.writerow([dt, h, f'{e1:.6e}', f'{e2:.6e}', 'Success']) + + print("Saved: mms_results.csv") + + # Save results to TXT file + print("Saving results to TXT file...") + with open('mms_results.txt', 'w') as txtfile: + txtfile.write("="*80 + "\n") + txtfile.write("MMS CONVERGENCE TEST RESULTS\n") + txtfile.write("="*80 + "\n\n") + + for dt in dts: + txtfile.write(f"dt = {dt}\n") + txtfile.write("-" * 60 + "\n") + txtfile.write(f"{'h':<10} {'e1':<25} {'e2':<25}\n") + txtfile.write("-" * 60 + "\n") + for h in sorted(results[dt].keys()): + if results[dt][h] == "Error": + txtfile.write(f"{h:<10} {'Failed':<25} {'Failed':<25}\n") + else: + e1, e2 = results[dt][h] + txtfile.write(f"{h:<10} {e1:<25.6e} {e2:<25.6e}\n") + txtfile.write("\n") + + print("Saved: mms_results.txt") + + # Create plots for each dt + print("\nGenerating convergence plots...") + + for dt in dts: + # Filter out errors and get successful runs + successful_h = [] + successful_e1 = [] + successful_e2 = [] + + for h in sorted(results[dt].keys()): + if results[dt][h] != "Error": + e1, e2 = results[dt][h] + successful_h.append(h) + successful_e1.append(e1) + successful_e2.append(e2) + + if successful_h: # Only plot if we have data + plt.figure(figsize=(10, 6)) + plt.loglog(successful_h, successful_e1, 'o-', label='e1 error', linewidth=2, markersize=8) + plt.loglog(successful_h, successful_e2, 's-', label='e2 error', linewidth=2, markersize=8) + plt.xlabel('Mesh size h', fontsize=12) + plt.ylabel('Error', fontsize=12) + plt.title(f'MMS Convergence for dt = {dt}', fontsize=14) + plt.legend(fontsize=11) + plt.grid(True, which="both", ls="-", alpha=0.3) + plt.tight_layout() + plt.savefig(f'mms_convergence_dt_{dt}.png', dpi=150) + print(f"Saved plot: mms_convergence_dt_{dt}.png") + plt.close() + + print("\nAll plots generated successfully!") diff --git a/check_ricker.py b/check_ricker.py new file mode 100644 index 000000000..e69de29bb diff --git a/example_run.py b/example_run.py new file mode 100644 index 000000000..722c2954a --- /dev/null +++ b/example_run.py @@ -0,0 +1,199 @@ +import numpy as np +import os +import spyro +import time + +from firedrake.petsc import PETSc +from math import pi as PI +from mpi4py import MPI +from numpy.linalg import norm +from scipy.integrate import quad + +from spyro.utils.file_utils import mkdir_timestamp + +opts = PETSc.Options() + +moment_tensor = opts.getBool("moment_tensor", False) +length = 1500 +h = 25 +L = opts.getReal("L", default=length) # Length (m) +h = opts.getReal("h", default=h) # Number of elements in each direction + +alpha = opts.getReal("alpha", default=1500) # P-wave velocity +beta = opts.getReal("beta", default=150) # S-wave velocity +rho = opts.getReal("rho", default=150) # Density (kg/m3) + +smag = opts.getReal("amplitude", default=1) # Source amplitude +f0 = opts.getReal("frequency", default=5) # Frequency (Hz) +t0 = 1/f0 # Time delay + +final_time = 1.0 +dt = 0.0001 +tn = opts.getReal("total_time", default=final_time) # Simulation time (s) +nt = opts.getInt("time_steps", default=final_time/dt+1) # Number of time steps +time_step = (tn/(nt - 1)) +out_freq = int(100) + +x_s = np.r_[-L/2, L/2] # Source location (m) +x_r = np.r_[opts.getReal("receiver_x", default=0.0), # Receiver relative location (m) + opts.getReal("receiver_z", default=150.0)] + + +def analytical_solution(i, j): + t = np.linspace(0, tn, nt) + + u_near = np.zeros(nt) # near field contribution + P_far = np.zeros(nt) # P-wave far-field + S_far = np.zeros(nt) # S-wave far field + + r = np.linalg.norm(x_r) + gamma_i = x_r[i]/r + gamma_j = x_r[j]/r + delta_ij = 1 if i == j else 0 + + def X0(t): + a = PI*f0*(t - t0) + return (1 - 2*a**2)*np.exp(-a**2) + + for k in range(nt): + res = quad(lambda tau: tau*X0(t[k] - tau), r/alpha, r/beta) + u_near[k] = smag*(1./(4*PI*rho))*(3*gamma_i*gamma_j - delta_ij)*(1./r**3)*res[0] + P_far[k] = smag*(1./(4*PI*rho*alpha**2))*gamma_i*gamma_j*(1./r)*X0(t[k] - r/alpha) + S_far[k] = smag*(1./(4*PI*rho*beta**2))*(gamma_i*gamma_j - delta_ij)*(1./r)*X0(t[k] - r/beta) + + return u_near + P_far - S_far + + +def explosive_source_analytical(i): + t = np.linspace(0, tn, nt) + + P_mid = np.zeros(nt) # P wave intermediate field + P_far = np.zeros(nt) # P wave far field + + r = np.linalg.norm(x_r) + gamma_i = x_r[i]/r + + def w(t): + a = PI*f0*(t - t0) + return (t - t0)*np.exp(-a**2) + + def w_dot(t): + a = PI*f0*(t - t0) + return (1 - 2*a**2)*np.exp(-a**2) + + for k in range(nt): + P_mid[k] = smag*(gamma_i/(4*PI*rho*alpha**2))*(1./r**2)*w(t[k] - r/alpha) + P_far[k] = smag*(gamma_i/(4*PI*rho*alpha**3))*(1./r)*w_dot(t[k] - r/alpha) + + return P_mid + P_far + + +def numerical_solution(j): + if moment_tensor: + A0 = smag*np.eye(2) + else: + A0 = np.zeros(2) + A0[j] = smag + + d = {} + + d["options"] = { + "cell_type": "T", + "variant": "lumped", + "degree": 4, + "dimension": 2, + } + + d["parallelism"] = { + "type": "automatic", + } + + d["mesh"] = { + "Lz": L, + "Lx": L, + "Ly": 0.0, + "mesh_type": "firedrake_mesh", + } + + d["acquisition"] = { + "source_type": "ricker", + "source_locations": [x_s.tolist()], + "frequency": f0, + "delay": t0, + "delay_type": "time", + "amplitude": A0, + "receiver_locations": [(x_s + x_r).tolist()], + } + + d["synthetic_data"] = { + "type": "object", + "density": rho, + "p_wave_velocity": alpha, + "s_wave_velocity": beta, + "real_velocity_file": None, + } + + d["time_axis"] = { + "initial_time": 0, + "final_time": tn, + "dt": time_step, + "output_frequency": 100, + "gradient_sampling_frequency": 100, + } + + d["visualization"] = { + "forward_output": True, + "forward_output_filename": "results/elastic_analytical.pvd", + "fwi_velocity_model_output": False, + "gradient_output": False, + "adjoint_output": False, + "debug_output": False, + } + + d["absorving_boundary_conditions"] = { + "status": False, + } + + wave = spyro.IsotropicWave(d) + wave.set_mesh(input_mesh_parameters={'edge_length': h, 'periodic': True}) + wave.forward_solve() + spyro.plots.plot_receiver(wave, xi=0, filename="rec_at_dir0.png") + spyro.plots.plot_receiver(wave, xi=1, filename="rec_at_dir1.png") + + # return wave.receivers_output.reshape(nt + 1, 2)[0:-1, :] + return wave.receivers_output + + +def err(u_a, u_n): + return norm(u_a - u_n)/norm(u_a) + + +if __name__ == "__main__": + rank = MPI.COMM_WORLD.Get_rank() + if rank == 0: + start = time.time() + + j = 0 + U_x = analytical_solution(0, j) + U_z = analytical_solution(1, j) + + u_n = numerical_solution(j) + u_x = u_n[:, 0, 0] + u_z = u_n[:, 0, 1] + + if rank == 0: + end = time.time() + print(f"Elapsed time: {end - start:.2f} s") + print(f"err_x = {err(U_x, u_x)}") + print(f"err_z = {err(U_z, u_z)}") + + if moment_tensor: + basename = "ExplosiveSource" + else: + basename = "ForceSource" + output_dir = mkdir_timestamp(basename) + np.save(os.path.join(output_dir, "x_analytical.npy"), U_x) + np.save(os.path.join(output_dir, "z_analytical.npy"), U_z) + np.save(os.path.join(output_dir, "numerical.npy"), u_n) + + print(f"Data saved to: {output_dir}") diff --git a/example_run_gmsh.py b/example_run_gmsh.py new file mode 100644 index 000000000..1bd33b6fe --- /dev/null +++ b/example_run_gmsh.py @@ -0,0 +1,274 @@ +import numpy as np +import os +import spyro +import time +import gmsh + +from firedrake.petsc import PETSc +from math import pi as PI +from mpi4py import MPI +from numpy.linalg import norm +from scipy.integrate import quad + +from spyro.utils.file_utils import mkdir_timestamp + +opts = PETSc.Options() + +moment_tensor = opts.getBool("moment_tensor", False) +length = 1500 +h = 25 +L = opts.getReal("L", default=length) # Length (m) +h = opts.getReal("h", default=h) # Number of elements in each direction + +alpha = opts.getReal("alpha", default=1500) # P-wave velocity +beta = opts.getReal("beta", default=150) # S-wave velocity +rho = opts.getReal("rho", default=150) # Density (kg/m3) + +smag = opts.getReal("amplitude", default=1) # Source amplitude +f0 = opts.getReal("frequency", default=5) # Frequency (Hz) +t0 = 1/f0 # Time delay + +final_time = 3.0 +dt = 0.0001 +tn = opts.getReal("total_time", default=final_time) # Simulation time (s) +nt = opts.getInt("time_steps", default=final_time/dt+1) # Number of time steps +time_step = (tn/(nt - 1)) +out_freq = int(100) + +x_s = np.r_[-L/2, L/2] # Source location (m) +x_r = np.r_[opts.getReal("receiver_x", default=0.0), # Receiver relative location (m) + opts.getReal("receiver_z", default=150.0)] + + + +def generate_gmsh_mesh(Lx, Ly, edge_size, radius=200, + center=None, + output_file="mesh.msh"): + + gmsh.initialize() + gmsh.model.add("refined_mesh") + + # ---------------------------- + # 1. Center of refinement + # ---------------------------- + if center is None: + cx, cy = -Lx/2, Ly/2 + else: + cx, cy = center + + # ---------------------------- + # 2. Create rectangle (OCC) + # ---------------------------- + p1 = gmsh.model.occ.addPoint(-Lx, 0, 0) + p2 = gmsh.model.occ.addPoint(0, 0, 0) + p3 = gmsh.model.occ.addPoint(0, Ly, 0) + p4 = gmsh.model.occ.addPoint(-Lx, Ly, 0) + + l1 = gmsh.model.occ.addLine(p1, p2) + l2 = gmsh.model.occ.addLine(p2, p3) + l3 = gmsh.model.occ.addLine(p3, p4) + l4 = gmsh.model.occ.addLine(p4, p1) + + loop = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4]) + surf = gmsh.model.occ.addPlaneSurface([loop]) + + # Refinement target point (must be embedded) + p_center = gmsh.model.occ.addPoint(cx, cy, 0) + gmsh.model.occ.synchronize() + gmsh.model.mesh.embed(0, [p_center], 2, surf) + + # ---------------------------- + # 3. Attractor field + # ---------------------------- + f_attr = gmsh.model.mesh.field.add("Attractor") + gmsh.model.mesh.field.setNumbers(f_attr, "NodesList", [p_center]) + + # ---------------------------- + # 4. Threshold field (refinement) + # ---------------------------- + f_thresh = gmsh.model.mesh.field.add("Threshold") + gmsh.model.mesh.field.setNumber(f_thresh, "InField", f_attr) + gmsh.model.mesh.field.setNumber(f_thresh, "SizeMin", edge_size * 0.2) + gmsh.model.mesh.field.setNumber(f_thresh, "SizeMax", edge_size) + gmsh.model.mesh.field.setNumber(f_thresh, "DistMin", radius * 0.3) + gmsh.model.mesh.field.setNumber(f_thresh, "DistMax", radius) + + gmsh.model.mesh.field.setAsBackgroundMesh(f_thresh) + + # Safe options that exist in all modern Gmsh versions + gmsh.option.setNumber("Mesh.CharacteristicLengthMin", edge_size * 0.2) + gmsh.option.setNumber("Mesh.CharacteristicLengthMax", edge_size) + + # ---------------------------- + # 5. Generate mesh + # ---------------------------- + gmsh.model.mesh.generate(2) + gmsh.model.mesh.optimize("Netgen") + + gmsh.write(output_file) + gmsh.finalize() + + print(f"Mesh saved: {output_file}") + + +def analytical_solution(i, j): + t = np.linspace(0, tn, nt) + + u_near = np.zeros(nt) # near field contribution + P_far = np.zeros(nt) # P-wave far-field + S_far = np.zeros(nt) # S-wave far field + + r = np.linalg.norm(x_r) + gamma_i = x_r[i]/r + gamma_j = x_r[j]/r + delta_ij = 1 if i == j else 0 + + def X0(t): + a = PI*f0*(t - t0) + return (1 - 2*a**2)*np.exp(-a**2) + + for k in range(nt): + res = quad(lambda tau: tau*X0(t[k] - tau), r/alpha, r/beta) + u_near[k] = smag*(1./(4*PI*rho))*(3*gamma_i*gamma_j - delta_ij)*(1./r**3)*res[0] + P_far[k] = smag*(1./(4*PI*rho*alpha**2))*gamma_i*gamma_j*(1./r)*X0(t[k] - r/alpha) + S_far[k] = smag*(1./(4*PI*rho*beta**2))*(gamma_i*gamma_j - delta_ij)*(1./r)*X0(t[k] - r/beta) + + return u_near + P_far - S_far + + +def explosive_source_analytical(i): + t = np.linspace(0, tn, nt) + + P_mid = np.zeros(nt) # P wave intermediate field + P_far = np.zeros(nt) # P wave far field + + r = np.linalg.norm(x_r) + gamma_i = x_r[i]/r + + def w(t): + a = PI*f0*(t - t0) + return (t - t0)*np.exp(-a**2) + + def w_dot(t): + a = PI*f0*(t - t0) + return (1 - 2*a**2)*np.exp(-a**2) + + for k in range(nt): + P_mid[k] = smag*(gamma_i/(4*PI*rho*alpha**2))*(1./r**2)*w(t[k] - r/alpha) + P_far[k] = smag*(gamma_i/(4*PI*rho*alpha**3))*(1./r)*w_dot(t[k] - r/alpha) + + return P_mid + P_far + + +def numerical_solution(j): + if moment_tensor: + A0 = smag*np.eye(2) + else: + A0 = np.zeros(2) + A0[j] = smag + + # Generate mesh using gmsh + mesh_file = "elastic_mesh.msh" + generate_gmsh_mesh(L, L, h, output_file=mesh_file) + + d = {} + + d["options"] = { + "cell_type": "T", + "variant": "lumped", + "degree": 4, + "dimension": 2, + } + + d["parallelism"] = { + "type": "automatic", + } + + d["mesh"] = { + "Lz": L, + "Lx": L, + "Ly": 0.0, + "mesh_type": "file", + "mesh_file": mesh_file, + } + + d["acquisition"] = { + "source_type": "ricker", + "source_locations": [x_s.tolist()], + "frequency": f0, + "delay": t0, + "delay_type": "time", + "amplitude": A0, + "receiver_locations": [(x_s + x_r).tolist()], + } + + d["synthetic_data"] = { + "type": "object", + "density": rho, + "p_wave_velocity": alpha, + "s_wave_velocity": beta, + "real_velocity_file": None, + } + + d["time_axis"] = { + "initial_time": 0, + "final_time": tn, + "dt": time_step, + "output_frequency": 100, + "gradient_sampling_frequency": 100, + } + + d["visualization"] = { + "forward_output": True, + "forward_output_filename": "results/elastic_analytical.pvd", + "fwi_velocity_model_output": False, + "gradient_output": False, + "adjoint_output": False, + "debug_output": False, + } + + d["absorving_boundary_conditions"] = { + "status": False, + } + + wave = spyro.IsotropicWave(d) + wave.forward_solve() + spyro.plots.plot_receiver(wave, xi=0, filename="rec_at_dir0.png") + spyro.plots.plot_receiver(wave, xi=1, filename="rec_at_dir1.png") + + return wave.receivers_output + + +def err(u_a, u_n): + return norm(u_a - u_n)/norm(u_a) + + +if __name__ == "__main__": + rank = MPI.COMM_WORLD.Get_rank() + if rank == 0: + start = time.time() + + j = 0 + U_x = analytical_solution(0, j) + U_z = analytical_solution(1, j) + + u_n = numerical_solution(j) + u_x = u_n[:, 0, 0] + u_z = u_n[:, 0, 1] + + if rank == 0: + end = time.time() + print(f"Elapsed time: {end - start:.2f} s") + print(f"err_x = {err(U_x, u_x)}") + print(f"err_z = {err(U_z, u_z)}") + + if moment_tensor: + basename = "ExplosiveSource" + else: + basename = "ForceSource" + output_dir = mkdir_timestamp(basename) + np.save(os.path.join(output_dir, "x_analytical.npy"), U_x) + np.save(os.path.join(output_dir, "z_analytical.npy"), U_z) + np.save(os.path.join(output_dir, "numerical.npy"), u_n) + + print(f"Data saved to: {output_dir}") diff --git a/mms_results.csv b/mms_results.csv new file mode 100644 index 000000000..1256deadc --- /dev/null +++ b/mms_results.csv @@ -0,0 +1,5 @@ +dt,h,e1,e2,status +0.0001,0.025,1.935904e-08,1.456845e-08,Success +0.0001,0.05,1.935909e-08,1.456846e-08,Success +0.0001,0.1,1.935333e-08,1.456966e-08,Success +0.0001,0.5,1.786389e-08,1.425244e-08,Success diff --git a/mms_results.txt b/mms_results.txt new file mode 100644 index 000000000..6af5c6576 --- /dev/null +++ b/mms_results.txt @@ -0,0 +1,54 @@ +================================================================================ +MMS CONVERGENCE TEST RESULTS +================================================================================ + +dt = 0.01 +------------------------------------------------------------ +h e1 e2 +------------------------------------------------------------ +0.01 Failed Failed +0.025 Failed Failed +0.05 Failed Failed +0.1 Failed Failed +0.5 1.612843e-04 1.226219e-04 + +dt = 0.001 +------------------------------------------------------------ +h e1 e2 +------------------------------------------------------------ +0.01 Failed Failed +0.025 1.915329e-06 1.436019e-06 +0.05 1.915312e-06 1.436026e-06 +0.1 1.914611e-06 1.436140e-06 +0.5 1.767806e-06 1.404177e-06 + +dt = 0.0001 +------------------------------------------------------------ +h e1 e2 +------------------------------------------------------------ +0.01 1.935890e-08 1.456830e-08 +0.025 1.935901e-08 1.456855e-08 +0.05 1.935914e-08 1.456846e-08 +0.1 1.935333e-08 1.456966e-08 +0.5 1.786389e-08 1.425244e-08 + +dt = 1e-05 +------------------------------------------------------------ +h e1 e2 +------------------------------------------------------------ +0.025 3.304003e-10 4.908308e-10 +0.05 5.036554e-10 5.879650e-10 +0.1 5.685461e-10 5.463560e-10 +0.5 6.978163e-10 9.519788e-10 + +Exemplo mais simples (idea do Ruben) e se ater ao grau 2 inicialmente: +- Simulação de uma barra de um único elemento, fazer convergence em h em uma unica direção da barra (parecido com artigo do Ruben de HABC). Espessura reduzida. +- Depois converter 1d (fake2d) para teste. +- Modelo com solução analítica seno e cosseno. +- Se tiver problema com 2+ aumentar uma fileira de elementos. +- Usar solução acústica +- Olhar frequencia natural +- Voltar para as origens!!! +- Revisar materia antiga de elastico +- Placa retangular para testar também. +- Verificar convergencia de freq natural. diff --git a/mms_results_1.csv b/mms_results_1.csv new file mode 100644 index 000000000..11df7b6ab --- /dev/null +++ b/mms_results_1.csv @@ -0,0 +1,9 @@ +dt,h,e1,e2,status +0.01,0.025,Error,Error,Failed +0.01,0.05,Error,Error,Failed +0.01,0.1,Error,Error,Failed +0.01,0.5,1.612843e-04,1.226219e-04,Success +0.001,0.025,1.915329e-06,1.436019e-06,Success +0.001,0.05,1.915312e-06,1.436026e-06,Success +0.001,0.1,1.914611e-06,1.436140e-06,Success +0.001,0.5,1.767806e-06,1.404177e-06,Success diff --git a/mms_results_1.txt b/mms_results_1.txt new file mode 100644 index 000000000..cc1980137 --- /dev/null +++ b/mms_results_1.txt @@ -0,0 +1,32 @@ +================================================================================ +MMS CONVERGENCE TEST RESULTS +================================================================================ + +dt = 0.01 +------------------------------------------------------------ +h e1 e2 +------------------------------------------------------------ +0.025 Failed Failed +0.05 Failed Failed +0.1 Failed Failed +0.5 1.612843e-04 1.226219e-04 + +dt = 0.001 +------------------------------------------------------------ +h e1 e2 +------------------------------------------------------------ +0.025 1.915329e-06 1.436019e-06 +0.05 1.915312e-06 1.436026e-06 +0.1 1.914611e-06 1.436140e-06 +0.5 1.767806e-06 1.404177e-06 + +dt = 0.0001 +------------------------------------------------------------ +h e1 e2 +------------------------------------------------------------ +0.025 1.935904e-08 1.456845e-08 +0.05 1.935909e-08 1.456846e-08 +0.1 1.935333e-08 1.456966e-08 +0.5 1.786389e-08 1.425244e-08 + + diff --git a/spyro/examples/camembert_elastic.py b/spyro/examples/camembert_elastic.py index c95891682..458cbe0e1 100644 --- a/spyro/examples/camembert_elastic.py +++ b/spyro/examples/camembert_elastic.py @@ -24,7 +24,7 @@ time_step = 2e-4 # [s] final_time = 1.5 # [s] -out_freq = int(0.01/time_step) +out_freq = int(0.1/time_step) n = 20 mesh = fire.RectangleMesh(n, n, 0, L, originX=-L, diagonal='crossed') @@ -57,7 +57,7 @@ "source_type": "ricker", "source_locations": source_locations, "frequency": freq, - "delay": 0, + "delay": 1/freq, "delay_type": "time", "amplitude": np.array([0, smag]), # "amplitude": smag * np.eye(2), @@ -78,7 +78,7 @@ "final_time": final_time, "dt": time_step, "output_frequency": out_freq, - "gradient_sampling_frequency": 1, + "gradient_sampling_frequency": out_freq, } d["visualization"] = { diff --git a/spyro/examples/elastic_analytical.py b/spyro/examples/elastic_analytical.py new file mode 100644 index 000000000..8e4b1b58a --- /dev/null +++ b/spyro/examples/elastic_analytical.py @@ -0,0 +1,205 @@ +''' +The analytical solutions are provided by the book: +Quantitative Seismology (2nd edition) from Aki and Richards +''' +import numpy as np +import os +import spyro +import time + +from firedrake.petsc import PETSc +from math import pi as PI +from mpi4py import MPI +from numpy.linalg import norm +from scipy.integrate import quad + +from spyro.utils.file_utils import mkdir_timestamp + +opts = PETSc.Options() + +moment_tensor = opts.getBool("moment_tensor", False) + +L = opts.getReal("L", default=450) # Length (m) +N = opts.getInt("N", default=30) # Number of elements in each direction +h = L/N # Element size (m) + +alpha = opts.getReal("alpha", default=1500) # P-wave velocity +beta = opts.getReal("beta", default=1000) # S-wave velocity +rho = opts.getReal("rho", default=2000) # Density (kg/m3) + +smag = opts.getReal("amplitude", default=1e3) # Source amplitude +f0 = opts.getReal("frequency", default=20) # Frequency (Hz) +t0 = 1/f0 # Time delay + +tn = opts.getReal("total_time", default=0.3) # Simulation time (s) +nt = opts.getInt("time_steps", default=750) # Number of time steps +time_step = (tn/nt) +out_freq = int(0.01/time_step) + +x_s = np.r_[-L/2, L/2, L/2] # Source location (m) +x_r = np.r_[opts.getReal("receiver_x", default=100), # Receiver relative location (m) + opts.getReal("receiver_y", default=0), + opts.getReal("receiver_z", default=100)] + + +def analytical_solution(i, j): + t = np.linspace(0, tn, nt) + + u_near = np.zeros(nt) # near field contribution + P_far = np.zeros(nt) # P-wave far-field + S_far = np.zeros(nt) # S-wave far field + + r = np.linalg.norm(x_r) + gamma_i = x_r[i]/r + gamma_j = x_r[j]/r + delta_ij = 1 if i == j else 0 + + def X0(t): + a = PI*f0*(t - t0) + return (1 - 2*a**2)*np.exp(-a**2) + + for k in range(nt): + res = quad(lambda tau: tau*X0(t[k] - tau), r/alpha, r/beta) + u_near[k] = smag*(1./(4*PI*rho))*(3*gamma_i*gamma_j - delta_ij)*(1./r**3)*res[0] + P_far[k] = smag*(1./(4*PI*rho*alpha**2))*gamma_i*gamma_j*(1./r)*X0(t[k] - r/alpha) + S_far[k] = smag*(1./(4*PI*rho*beta**2))*(gamma_i*gamma_j - delta_ij)*(1./r)*X0(t[k] - r/beta) + + return u_near + P_far - S_far + + +def explosive_source_analytical(i): + t = np.linspace(0, tn, nt) + + P_mid = np.zeros(nt) # P wave intermediate field + P_far = np.zeros(nt) # P wave far field + + r = np.linalg.norm(x_r) + gamma_i = x_r[i]/r + + def w(t): + a = PI*f0*(t - t0) + return (t - t0)*np.exp(-a**2) + + def w_dot(t): + a = PI*f0*(t - t0) + return (1 - 2*a**2)*np.exp(-a**2) + + for k in range(nt): + P_mid[k] = smag*(gamma_i/(4*PI*rho*alpha**2))*(1./r**2)*w(t[k] - r/alpha) + P_far[k] = smag*(gamma_i/(4*PI*rho*alpha**3))*(1./r)*w_dot(t[k] - r/alpha) + + return P_mid + P_far + + +def numerical_solution(j): + if moment_tensor: + A0 = smag*np.eye(3) + else: + A0 = np.zeros(3) + A0[j] = smag + + d = {} + + d["options"] = { + "cell_type": "T", + "variant": "lumped", + "degree": 3, + "dimension": 3, + } + + d["parallelism"] = { + "type": "automatic", + } + + d["mesh"] = { + "Lz": L, + "Lx": L, + "Ly": L, + "mesh_file": None, + "mesh_type": "firedrake_mesh", + } + + d["acquisition"] = { + "source_type": "ricker", + "source_locations": [x_s.tolist()], + "frequency": f0, + "delay": t0, + "delay_type": "time", + "amplitude": A0, + "receiver_locations": [(x_s + x_r).tolist()], + } + + d["synthetic_data"] = { + "type": "object", + "density": rho, + "p_wave_velocity": alpha, + "s_wave_velocity": beta, + "real_velocity_file": None, + } + + d["time_axis"] = { + "initial_time": 0, + "final_time": tn, + "dt": time_step, + "output_frequency": out_freq, + "gradient_sampling_frequency": 1, + } + + d["visualization"] = { + "forward_output": True, + "forward_output_filename": "results/elastic_analytical.pvd", + "fwi_velocity_model_output": False, + "gradient_output": False, + "adjoint_output": False, + "debug_output": False, + } + + d["absorving_boundary_conditions"] = { + "status": True, + "damping_type": "local", + } + + wave = spyro.IsotropicWave(d) + wave.set_mesh(mesh_parameters={'edge_length': h}) + wave.forward_solve() + + return wave.receivers_output.reshape(nt + 1, 3)[0:-1, :] + + +def err(u_a, u_n): + return norm(u_a - u_n)/norm(u_a) + + +if __name__ == "__main__": + rank = MPI.COMM_WORLD.Get_rank() + if rank == 0: + start = time.time() + + j = 0 + U_x = analytical_solution(0, j) + U_y = analytical_solution(1, j) + U_z = analytical_solution(2, j) + + u_n = numerical_solution(j) + u_x = u_n[:, 0] + u_y = u_n[:, 1] + u_z = u_n[:, 2] + + if rank == 0: + end = time.time() + print(f"Elapsed time: {end - start:.2f} s") + print(f"err_x = {err(U_x, u_x)}") + print(f"err_y = {err(U_y, u_y)}") + print(f"err_z = {err(U_z, u_z)}") + + if moment_tensor: + basename = "ExplosiveSource" + else: + basename = "ForceSource" + output_dir = mkdir_timestamp(basename) + np.save(os.path.join(output_dir, "x_analytical.npy"), U_x) + np.save(os.path.join(output_dir, "y_analytical.npy"), U_y) + np.save(os.path.join(output_dir, "z_analytical.npy"), U_z) + np.save(os.path.join(output_dir, "numerical.npy"), u_n) + + print(f"Data saved to: {output_dir}") diff --git a/spyro/meshing/meshing_functions.py b/spyro/meshing/meshing_functions.py index c28dafca4..4964a8c05 100644 --- a/spyro/meshing/meshing_functions.py +++ b/spyro/meshing/meshing_functions.py @@ -222,8 +222,7 @@ def create_seismicmesh_mesh(self): if self.dimension == 2: return self.create_seimicmesh_2d_mesh() elif self.dimension == 3: - raise NotImplementedError("Not implemented yet") - # return self.create_seismicmesh_3D_mesh() + return self.create_seismicmesh_3d_mesh() else: raise ValueError("dimension is not supported") @@ -356,6 +355,68 @@ def create_seismicmesh_2D_mesh_homogeneous(self): return fire.Mesh(self.output_file_name) + def create_seismicmesh_3d_mesh(self): + """ + Creates a 3D mesh based on SeismicMesh meshing utilities. + """ + if self.velocity_model is None: + return self.create_seismicmesh_3D_mesh_homogeneous() + else: + raise NotImplementedError("Not implemented yet") + + def create_seismicmesh_3D_mesh_homogeneous(self): + """ + Creates a 3D mesh based on SeismicMesh meshing utilities, with homogeneous velocity model. + """ + Lz = self.length_z + Lx = self.length_x + Ly = self.length_y + pad = self.abc_pad + real_lz = Lz + pad + real_lx = Lx + 2 * pad + real_ly = Ly + 2 * pad + + bbox = (-real_lz, 0.0, + -pad, real_lx - pad, + -pad, real_ly - pad) + cube = SeismicMesh.Cube(bbox) + + points, cells = SeismicMesh.generate_mesh( + domain=cube, + edge_length=self.edge_length, + verbose=0, + max_iter=75, + ) + + points, cells = SeismicMesh.sliver_removal( + points=points, + domain=cube, + edge_length=self.edge_length, + preserve=True, + ) + + meshio.write_points_cells( + self.output_file_name, + points, + [("tetra", cells)], + file_format="gmsh22", + binary=False, + ) + + meshio.write_points_cells( + self.output_file_name + ".vtk", + points, + [("tetra", cells)], + file_format="vtk", + ) + + return fire.Mesh( + self.output_file_name, + distribution_parameters={ + "overlap_type": (fire.DistributedMeshOverlapType.NONE, 0) + } + ) + def calculate_edge_length(cpw, minimum_velocity, frequency): v_min = minimum_velocity @@ -365,25 +426,6 @@ def calculate_edge_length(cpw, minimum_velocity, frequency): edge_length = lbda_min/cpw return edge_length -# def create_firedrake_3D_mesh_based_on_parameters(dx, cell_type): -# nx = int(self.length_x / dx) -# nz = int(self.length_z / dx) -# ny = int(self.length_y / dx) -# if self.cell_type == "quadrilateral": -# quadrilateral = True -# else: -# quadrilateral = False - -# return spyro.BoxMesh( -# nz, -# nx, -# ny, -# self.length_z, -# self.length_x, -# self.length_y, -# quadrilateral=quadrilateral, -# ) - def RectangleMesh(nx, ny, Lx, Ly, pad=None, comm=None, quadrilateral=False): """Create a rectangle mesh based on the Firedrake mesh. @@ -460,7 +502,7 @@ def PeriodicRectangleMesh( else: pad = 0 mesh = fire.PeriodicRectangleMesh( - nx, ny, Lx, Ly, quadrilateral=quadrilateral, comm=comm + nx, ny, Lx, Ly, quadrilateral=quadrilateral, comm=comm, diagonal="crossed", ) mesh.coordinates.dat.data[:, 0] *= -1.0 mesh.coordinates.dat.data[:, 1] -= pad diff --git a/spyro/plots/__init__.py b/spyro/plots/__init__.py index e3dc7c3b3..740667f0b 100644 --- a/spyro/plots/__init__.py +++ b/spyro/plots/__init__.py @@ -1,5 +1,5 @@ from .plots import plot_shots, plot_mesh_sizes, plot_model, plot_function, plot_model_in_p1 -from .plots import debug_plot, debug_pvd +from .plots import debug_plot, debug_pvd, plot_receiver __all__ = [ "plot_shots", @@ -9,4 +9,5 @@ "debug_plot", "debug_pvd", "plot_model_in_p1", + "plot_receiver", ] diff --git a/spyro/plots/plots.py b/spyro/plots/plots.py index f9960f36c..b86897299 100644 --- a/spyro/plots/plots.py +++ b/spyro/plots/plots.py @@ -567,3 +567,21 @@ def plot_model_in_p1(Wave_object, dx=0.01, filename="model.png", abc_points=None new_wave_obj.set_initial_velocity_model(conditional=Wave_object.initial_velocity_model) return plot_model(new_wave_obj, filename=filename, abc_points=abc_points, show=show, flip_axis=flip_axis) + + +def plot_receiver(wave, receiver_id=0, show=False, filename="receiver.png", xi=None): + nt = int(wave.final_time/wave.dt) + 1 + time_vector = np.linspace(0.0, wave.final_time, nt) + plt.close() + if xi is not None: + title = f"Receiver {receiver_id} at direction {xi}" + plt.plot(time_vector, wave.receivers_output[:, receiver_id, xi]) + else: + title = f"Receiver {receiver_id}" + plt.plot(time_vector, wave.receivers_output[:, receiver_id]) + plt.title(title) + if filename is not None: + plt.savefig(filename) + if show: + plt.show() + plt.close() diff --git a/spyro/receivers/dirac_delta_projector.py b/spyro/receivers/dirac_delta_projector.py index d132870f8..91acc73c8 100644 --- a/spyro/receivers/dirac_delta_projector.py +++ b/spyro/receivers/dirac_delta_projector.py @@ -492,10 +492,8 @@ def get_hexa_real_cell_node_map(V, mesh): cell_node_map = np.zeros( (layers * cells_per_layer, nodes_per_cell), dtype=int ) - print(f"cnm size : {np.shape(cell_node_map)}", flush=True) for layer in range(layers): - print(f"layer : {layer} of {layers}", flush=True) for cell in range(cells_per_layer): cnm_base = weird_cnm[cell] cell_id = layer + layers * cell diff --git a/spyro/solvers/acoustic_wave.py b/spyro/solvers/acoustic_wave.py index b47f81b40..606823661 100644 --- a/spyro/solvers/acoustic_wave.py +++ b/spyro/solvers/acoustic_wave.py @@ -210,3 +210,10 @@ def rhs_no_pml(self): return self.B.sub(0) else: return self.B + + @override + def check_stability(self): + assert ( + fire.norm(self.get_function()) < 1 + ), "Numerical instability. Try reducing dt or building the " \ + "mesh differently" diff --git a/spyro/solvers/elastic_wave/isotropic_wave.py b/spyro/solvers/elastic_wave/isotropic_wave.py index fd9b11f65..574e6a1a0 100644 --- a/spyro/solvers/elastic_wave/isotropic_wave.py +++ b/spyro/solvers/elastic_wave/isotropic_wave.py @@ -1,7 +1,7 @@ import numpy as np -from firedrake import (assemble, Constant, curl, DirichletBC, div, Function, - FunctionSpace, project, VectorFunctionSpace) +from firedrake import (assemble, norm, Constant, curl, DirichletBC, div, Function, + FunctionSpace, project, VectorFunctionSpace, interpolate) from .elastic_wave import ElasticWave from .forms import (isotropic_elastic_without_pml, @@ -51,17 +51,18 @@ def __init__(self, dictionary, comm=None): @override def initialize_model_parameters_from_object(self, synthetic_data_dict: dict): - def constant_wrapper(value): + def variable_wrapper(value, function_space=None): if np.isscalar(value): return Constant(value) else: return value - def get_value(key, default=None): - return constant_wrapper(synthetic_data_dict.get(key, default)) + def get_value(key, default=None, function_space=None): + return variable_wrapper(synthetic_data_dict.get(key, default), function_space=function_space) self.rho = get_value("density") - self.lmbda = get_value("lambda", default=get_value("lame_first")) + Q = FunctionSpace(self.mesh, self.function_space.ufl_element().sub_elements[0]) + self.lmbda = interpolate(synthetic_data_dict.get("lambda"), Q) self.mu = get_value("mu", get_value("lame_second")) self.c = get_value("p_wave_velocity") self.c_s = get_value("s_wave_velocity") @@ -80,7 +81,7 @@ def get_value(key, default=None): not bool(self.mu) if option_1: - self.c = ((self.lmbda + 2*self.mu)/self.rho)**0.5 + self.c = self.lmbda # ((self.lmbda + 2*self.mu)/self.rho)**0.5 self.c_s = (self.mu/self.rho)**0.5 elif option_2: self.mu = self.rho*self.c_s**2 @@ -234,3 +235,10 @@ def update_s_wave(self): self.s_wave.assign(project(curl(self.get_function()), self.C_h)) return self.s_wave + + @override + def check_stability(self): + assert ( + np.isfinite(norm(self.get_function())) + ), "Numerical instability. Try reducing dt or building the " \ + "mesh differently" diff --git a/spyro/solvers/time_integration_central_difference.py b/spyro/solvers/time_integration_central_difference.py index 882956d4e..2c5447d3a 100644 --- a/spyro/solvers/time_integration_central_difference.py +++ b/spyro/solvers/time_integration_central_difference.py @@ -39,7 +39,6 @@ def central_difference(wave, source_ids=[0]): usol_recv = [] save_step = 0 for step in range(nt): - # Basic way of applying sources wave.update_source_expression(t) fire.assemble(wave.rhs, tensor=wave.B) @@ -61,10 +60,7 @@ def central_difference(wave, source_ids=[0]): save_step += 1 if (step - 1) % wave.output_frequency == 0: - assert ( - fire.norm(wave.get_function()) < 1 - ), "Numerical instability. Try reducing dt or building the " \ - "mesh differently" + wave.check_stability() wave.field_logger.log(t) helpers.display_progress(wave.comm, t) diff --git a/spyro/solvers/wave.py b/spyro/solvers/wave.py index 2601d3f7a..73e6c089c 100644 --- a/spyro/solvers/wave.py +++ b/spyro/solvers/wave.py @@ -102,7 +102,7 @@ def __init__(self, dictionary=None, comm=None): self.field_logger.add_field("forward", self.get_function_name(), lambda: self.get_function()) - def forward_solve(self): + def forward_solve(self, build_matrix_operator=True): """Solves the forward problem.""" print("\nSolving Forward Problem") @@ -112,7 +112,9 @@ def forward_solve(self): if self.abc_boundary_layer_type != "hybrid": self._initialize_model_parameters() - self.matrix_building() + + if build_matrix_operator: + self.matrix_building() self.wave_propagator() def force_rebuild_function_space(self): @@ -432,3 +434,7 @@ def rhs_no_pml(self): the DOFs associated with the subspace of the original problem). ''' pass + + @abstractmethod + def check_stability(self): + pass diff --git a/spyro/sources/Sources.py b/spyro/sources/Sources.py index e8e4dafb0..4280d2ae2 100644 --- a/spyro/sources/Sources.py +++ b/spyro/sources/Sources.py @@ -70,11 +70,13 @@ def __init__(self, wave_object): self.amplitude = wave_object.amplitude self.is_local = [0] * self.number_of_points self.current_sources = None - self.update_wavelet(wave_object) if np.isscalar(self.amplitude) or (self.amplitude.size <= 3): + self.integral = False self.build_maps(order=0) else: + self.integral = True self.build_maps(order=1) + self.update_wavelet(wave_object) def update_wavelet(self, wave_object): self.wavelet = full_ricker_wavelet( @@ -83,6 +85,7 @@ def update_wavelet(self, wave_object): frequency=wave_object.frequency, delay=wave_object.delay, delay_type=wave_object.delay_type, + integral=self.integral ) def apply_source(self, rhs_forcing, step): @@ -113,17 +116,16 @@ def apply_source(self, rhs_forcing, step): return rhs_forcing -def timedependentSource(model, t, freq=None, amp=1, delay=1.5): +def timedependentSource(model, t, freq=None, delay=1.5, amplitude=1.0): if model["acquisition"]["source_type"] == "Ricker": - return ricker_wavelet(t, freq, amp, delay=delay) - # elif model["acquisition"]["source_type"] == "MMS": - # return MMS_time(t) + return ricker_wavelet(t, freq, delay=delay, amplitude=amplitude) else: raise ValueError("source not implemented") def ricker_wavelet( - t, freq, amp=1.0, delay=1.5, delay_type="multiples_of_minimun" + t, freq, delay=1.5, delay_type="multiples_of_minimun", + integral=False, amplitude=1.0, ): """Creates a Ricker source function with a delay in term of multiples of the distance @@ -135,8 +137,6 @@ def ricker_wavelet( Time freq: float Frequency of the wavelet - amp: float - Amplitude of the wavelet delay: float Delay in term of multiples of the distance between the minimums. @@ -157,7 +157,10 @@ def ricker_wavelet( t = t - time_delay # t = t - delay / freq tt = (math.pi * freq * t) ** 2 - return amp * (1.0 - (2.0) * tt) * math.exp((-1.0) * tt) + if integral: + return t*math.exp((-1.0) * tt) * amplitude + else: + return (1.0 - (2.0) * tt) * math.exp((-1.0) * tt) * amplitude def full_ricker_wavelet( @@ -167,6 +170,7 @@ def full_ricker_wavelet( cutoff=None, delay=1.5, delay_type="multiples_of_minimun", + integral=False ): """Compute the Ricker wavelet optionally applying low-pass filtering using cutoff frequency in Hertz. @@ -199,7 +203,7 @@ def full_ricker_wavelet( full_wavelet = np.zeros((nt,)) for t in range(nt): full_wavelet[t] = ricker_wavelet( - time, frequency, 1, delay=delay, delay_type=delay_type + time, frequency, delay=delay, delay_type=delay_type, integral=integral ) time += dt if cutoff is not None: diff --git a/spyro/utils/file_utils.py b/spyro/utils/file_utils.py new file mode 100644 index 000000000..81b90ccf1 --- /dev/null +++ b/spyro/utils/file_utils.py @@ -0,0 +1,10 @@ +import os + +from time import strftime, gmtime + + +def mkdir_timestamp(basename): + name = basename + "_" + strftime("%Y%m%d_%H%M%S", gmtime()) + if not os.path.exists(name): + os.mkdir(name) + return name diff --git a/temporal_convergence_analysis.txt b/temporal_convergence_analysis.txt new file mode 100644 index 000000000..3ee091910 --- /dev/null +++ b/temporal_convergence_analysis.txt @@ -0,0 +1,66 @@ +================================================================================ +TEMPORAL CONVERGENCE ANALYSIS +Expected convergence rate: 2.0 (second-order in time) +Tolerance: ± 0.1 +================================================================================ + +h = 0.01 +-------------------------------------------------------------------------------- +dt e1 rate(e1) e2 rate(e2) Status +-------------------------------------------------------------------------------- +0.0001 1.935890e-08 - 1.456830e-08 - - + + +h = 0.025 +-------------------------------------------------------------------------------- +dt e1 rate(e1) e2 rate(e2) Status +-------------------------------------------------------------------------------- +0.001 1.915329e-06 2.00 1.436019e-06 1.99 ✓ PASS +0.0001 1.935901e-08 1.77 1.456855e-08 1.47 ✗ FAIL +1e-05 3.304003e-10 - 4.908308e-10 - - + +Average convergence rate (e1): 1.88 +Average convergence rate (e2): 1.73 +Overall for h=0.025: ✗ NOT converging at expected rate + + +h = 0.05 +-------------------------------------------------------------------------------- +dt e1 rate(e1) e2 rate(e2) Status +-------------------------------------------------------------------------------- +0.001 1.915312e-06 2.00 1.436026e-06 1.99 ✓ PASS +0.0001 1.935914e-08 1.58 1.456846e-08 1.39 ✗ FAIL +1e-05 5.036554e-10 - 5.879650e-10 - - + +Average convergence rate (e1): 1.79 +Average convergence rate (e2): 1.69 +Overall for h=0.05: ✗ NOT converging at expected rate + + +h = 0.1 +-------------------------------------------------------------------------------- +dt e1 rate(e1) e2 rate(e2) Status +-------------------------------------------------------------------------------- +0.001 1.914611e-06 2.00 1.436140e-06 1.99 ✓ PASS +0.0001 1.935333e-08 1.53 1.456966e-08 1.43 ✗ FAIL +1e-05 5.685461e-10 - 5.463560e-10 - - + +Average convergence rate (e1): 1.76 +Average convergence rate (e2): 1.71 +Overall for h=0.1: ✗ NOT converging at expected rate + + +h = 0.5 +-------------------------------------------------------------------------------- +dt e1 rate(e1) e2 rate(e2) Status +-------------------------------------------------------------------------------- +0.01 1.612843e-04 1.96 1.226219e-04 1.94 ✓ PASS +0.001 1.767806e-06 2.00 1.404177e-06 1.99 ✓ PASS +0.0001 1.786389e-08 1.41 1.425244e-08 1.18 ✗ FAIL +1e-05 6.978163e-10 - 9.519788e-10 - - + +Average convergence rate (e1): 1.79 +Average convergence rate (e2): 1.70 +Overall for h=0.5: ✗ NOT converging at expected rate + + diff --git a/test_3d/test_elastic_isotropic.py b/test_3d/test_elastic_isotropic.py new file mode 100644 index 000000000..7d5fa77bb --- /dev/null +++ b/test_3d/test_elastic_isotropic.py @@ -0,0 +1,39 @@ +import math +import numpy as np +import pytest + +from spyro.examples.elastic_analytical import (analytical_solution, + numerical_solution) +from numpy.linalg import norm + +def err(u_a, u_n): + return norm(u_a - u_n)/norm(u_a) + + +def check_err(u_a, u_n, err_max, norm_max): + assert (err_max is None) or (norm_max is None) # Sanity check + if err_max is None: + assert math.isclose(norm(u_a), 0.0) + assert norm(u_n) < norm_max + else: + assert err(u_a, u_n) < err_max + + +@pytest.mark.parametrize("j,err_x_max,err_y_max,err_z_max,norm_x_max,norm_y_max,norm_z_max",[ + (0, 0.03, None, 0.03, None, 2e-12, None), + (1, None, 0.03, None, 2e-12, None, 2e-12), + (2, 0.03, None, 0.03, None, 2e-12, None)]) +def test_force(j, err_x_max, err_y_max, err_z_max, + norm_x_max, norm_y_max, norm_z_max): + U_x = analytical_solution(0, j) + U_y = analytical_solution(1, j) + U_z = analytical_solution(2, j) + + u_n = numerical_solution(j) + u_x = u_n[:, 0] + u_y = u_n[:, 1] + u_z = u_n[:, 2] + + check_err(U_x, u_x, err_x_max, norm_x_max) + check_err(U_y, u_y, err_y_max, norm_y_max) + check_err(U_z, u_z, err_z_max, norm_z_max) diff --git a/tests/on_one_core/test_MMS.py b/tests/on_one_core/test_MMS.py index 455324ea9..f624ccfa3 100644 --- a/tests/on_one_core/test_MMS.py +++ b/tests/on_one_core/test_MMS.py @@ -1,9 +1,10 @@ import math +import numpy as np from copy import deepcopy from firedrake import * import spyro -from .model import dictionary as model +from tests.on_one_core.model import dictionary as model model["acquisition"]["source_type"] = "MMS" @@ -86,8 +87,115 @@ def test_isotropic_wave_2D(): assert math.isclose(e2, 0.0, abs_tol=1e-7) +def test_initial_conditions_elastic(): + + h = 0.25 + dt = 0.001 + # Acoustic-like solution: both components have the same scalar solution + scalar_solution = lambda x, t: x[0] * (x[0] + 1) * x[1] * (x[1] - 1) * t + u1 = lambda x, t: scalar_solution(x, t) + u2 = lambda x, t: scalar_solution(x, t) + u = lambda x, t: as_vector([u1(x, t), u2(x, t)]) + + # Body forces for acoustic-like case (both components get the same source) + scalar_source = lambda x, t: -(x[0]**2 + x[0] + x[1]**2 - x[1]) * 2 * t + b1 = lambda x, t: scalar_source(x, t) + b2 = lambda x, t: scalar_source(x, t) + b = lambda x, t: as_vector([b1(x, t), b2(x, t)]) + + fo = int(0.1/dt) + + dictionary = {} + dictionary["options"] = { + "cell_type": "Q", # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q) + "variant": "lumped", # lumped, equispaced or DG, default is lumped "method":"MLT", # (MLT/spectral_quadrilateral/DG_triangle/DG_quadrilateral) You can either specify a cell_type+variant or a method + "degree": 4, # p order + "dimension": 2, # dimension + } + dictionary["parallelism"] = { + "type": "automatic", # options: automatic (same number of cores for evey processor) or spatial + } + dictionary["mesh"] = { + "Lz": 1.0, # depth in km - always positive + "Lx": 1.0, # width in km - always positive + "Ly": 0.0, # thickness in km - always positive + "mesh_type": "firedrake_mesh", # options: firedrake_mesh or user_mesh + "mesh_file": None, # specify the mesh file + } + dictionary["acquisition"] = { + "source_type": "MMS", + "source_locations": [(-1.0, 1.0)], + "frequency": 5.0, + "delay": 1.5, + "receiver_locations": [(-0.0, 0.5)], + "body_forces": b, + } + dictionary["time_axis"] = { + "initial_condition": u, + "initial_time": 0.0, # Initial time for event + "final_time": 1.0, # Final time for event + "dt": dt, # timestep size + "amplitude": 1, # the Ricker has an amplitude of 1. + "output_frequency": fo, # how frequently to output solution to pvds + "gradient_sampling_frequency": 1, # how frequently to save solution to RAM + } + + dictionary["visualization"] = { + "forward_output": False, + "output_filename": "results/forward_output.pvd", + "fwi_velocity_model_output": False, + "velocity_model_filename": None, + "gradient_output": False, + "gradient_filename": None, + } + + dictionary["synthetic_data"] = { + "type": "object", + "density": 1, + "lambda": 1, + "mu": 0, # Set to 0 to make elastic equivalent to acoustic + "real_velocity_file": None, + } + dictionary["boundary_conditions"] = [ + ("uz", "on_boundary", 0), + ("ux", "on_boundary", 0), + ] + + wave = spyro.IsotropicWave(dictionary) + wave.set_mesh(input_mesh_parameters={"edge_length": h}) + # Lets build de matrix operator outside of the forward solve so we can + # set previous timesteps for the MMS problem + wave._initialize_model_parameters() + wave.matrix_building() + x = wave.get_spatial_coordinates() + V = wave.function_space + test_function = Function(V) + + # testing u_nm1 + test_function.interpolate(u(x, 0.0 - 2*dt)) + max_test = np.amax(test_function.dat.data[:]) + min_test = np.amin(test_function.dat.data[:]) + max_func = np.amax(wave.u_nm1.dat.data[:]) + min_func = np.amin(wave.u_nm1.dat.data[:]) + passed = math.isclose(max_test, max_func) and math.isclose(min_test, min_func) + print(f"Test of elastic MMS u_nm1 got correct starting value: {passed}", flush=True) + assert passed + + # testing u_n + test_function.interpolate(u(x, 0.0 - dt)) + max_test = np.amax(test_function.dat.data[:]) + min_test = np.amin(test_function.dat.data[:]) + max_func = np.amax(wave.u_n.dat.data[:]) + min_func = np.amin(wave.u_n.dat.data[:]) + passed = math.isclose(max_test, max_func) and math.isclose(min_test, min_func) + print(f"Test of elastic MMS u_n got correct starting value: {passed}", flush=True) + assert passed + + if __name__ == "__main__": test_method_triangles_lumped() test_method_quads_lumped() + test_initial_conditions_elastic() + test_isotropic_wave_2D() print("END") diff --git a/tests/on_one_core/test_sources.py b/tests/on_one_core/test_sources.py index aaeaf3b34..77e514973 100644 --- a/tests/on_one_core/test_sources.py +++ b/tests/on_one_core/test_sources.py @@ -21,7 +21,7 @@ def test_ricker_varies_in_time(): # tests if ricker starts at zero delay = 1.5 * math.sqrt(6.0) / (math.pi * frequency) t = 0.0 - r0 = spyro.sources.timedependentSource(modelRicker, t, frequency, amplitude) + r0 = spyro.sources.timedependentSource(modelRicker, t, frequency, amplitude=amplitude) test1 = math.isclose( r0, 0, @@ -32,7 +32,7 @@ def test_ricker_varies_in_time(): minimum = -amplitude * 2 / math.exp(3.0 / 2.0) t = 0.0 + delay + math.sqrt(6.0) / (2.0 * math.pi * frequency) rmin1 = spyro.sources.timedependentSource( - modelRicker, t, frequency, amplitude + modelRicker, t, frequency, amplitude=amplitude ) test2 = math.isclose( rmin1, @@ -41,14 +41,14 @@ def test_ricker_varies_in_time(): t = 0.0 + delay - math.sqrt(6.0) / (2.0 * math.pi * frequency) rmin2 = spyro.sources.timedependentSource( - modelRicker, t, frequency, amplitude + modelRicker, t, frequency, amplitude=amplitude ) test3 = math.isclose(rmin2, minimum) # tests if maximum value in correct and occurs at correct location t = 0.0 + delay rmax = spyro.sources.timedependentSource( - modelRicker, t, frequency, amplitude + modelRicker, t, frequency, amplitude=amplitude ) test4 = math.isclose( rmax,