Can I fd.tripcolor a part of the data in pressure field function?

[acse-yw11823]

Hello, I am struggling ploting the pressure field. I have my model code as below:

dictionary = {}

dictionary["options"] = {
    "cell_type": "T",  # simplexes such as triangles or tetrahedra (T) or quadrilaterals (Q)
    "variant": "lumped",  # lumped, equispaced or DG, default is lumped
    "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   # Como ver isso sem ler a malha?
    "Lx": 1.0,  # width in km - always positive
    "Ly": 0.0,  # thickness in km - always positive
    "mesh_file": None,
    "mesh_type": "firedrake_mesh", # default use triangler mesh 
}
# dictionary["acquisition"] = {
#     "source_type": "ricker",
#     "source_locations": [(-0.1, 0.5)],
#     "frequency": 5.0,
#     "delay": 1.5,
#     # "delay_type": "time",
#     "receiver_locations":  spyro.create_transect((-0.9, 0), (-0.9, 1), 101),
#     # "receiver_locations": [(-0.5, 0.5)]
# }
dictionary["acquisition"] = {
    "source_type": "ricker",
    "source_locations": [(-0.1, 0.5)], 
    "frequency": 10.0,
    "receiver_locations":  spyro.create_transect((-0.9, 0), (-0.9, 1), 101),
    "delay": 1.5
}
dictionary["time_axis"] = {
    "initial_time": 0.0,  # Initial time for event
    "final_time": 1.0,  # Final time for event
    "dt": 0.0005,  # timestep size
    "output_frequency": 100,  # how frequently to output solution to pvds - Perguntar Daiane ''post_processing_frequnecy'
    "gradient_sampling_frequency": 1,  # how frequently to save solution to RAM    - Perguntar Daiane 'gradient_sampling_frequency'
}
dictionary["visualization"] = {
    "forward_output": True,
    "forward_output_filename": "results/forward_output_00.pvd",
    "fwi_velocity_model_output": False,
    "velocity_model_filename": None,
    "gradient_output": False,
    "gradient_filename": "results/Gradient_00.pvd",
    "adjoint_output": False,
    "adjoint_filename": None,
    "debug_output": False,
}
# Add absorbing boundary conditions with a pad length of 0.1 km.
dictionary["absorving_boundary_conditions"] = {
    "status": True, 
    "damping_type": "PML",
    "exponent": 2,
    "cmax": 4.5,
    "R": 1e-6,
    "pad_length": 0.8,   # in km
}

Wave_obj = spyro.AcousticWave(dictionary=dictionary)
Wave_obj.set_mesh(mesh_parameters={"dx": 0.02, "periodic": True})

# Set initial velocity model
Wave_obj.set_initial_velocity_model(constant=1.5)

# Calculate time
start_time = time.time()  # 记录开始时间
Wave_obj.forward_solve()  # 调用 forward_solve() 方法
end_time = time.time()  # 记录结束时间
execution_time = end_time - start_time
print(f"Execution time for forward_solve: {execution_time:.6f} seconds")

# define some items
src_00 = Wave_obj.sources
src_id_00 = src_00.current_source
num_src_00 = src_00.number_of_points
src_position_00 = src_00.point_locations
receivers_00 = Wave_obj.receivers
num_recvs_00 = Wave_obj.number_of_receivers
# comm = Wave_obj.comm

dt_00 = Wave_obj.dt
tf_00 = Wave_obj.final_time
t0_00 = Wave_obj.initial_time

p_field_00 = Wave_obj.forward_solution
p_rec_00 = Wave_obj.forward_solution_receivers

src_wavelet_00 = Wave_obj.wavelet
# Save wavelet as a .dat file
np.savetxt("src_spyro.dat", Wave_obj.wavelet[:])

and I want to plot the p-field using:

fd.tripcolor(p_field_00[1000]) 

The picture looks like the following:

. However, I only want to display a part of this image, specifically the region from 0 to -1 instead of the entire range from 0 to -1.75. How can I use fd.tripcolor to crop this section of the graph?

Many thanks for the help!!!!

[Ig-dolci]

I think ParaView will simplify the section crop for you. Just write the following code:

output_file = VTKFile("file_name.pvd")
output_file.write(your_function)

[ksagiyam]

Firedrake plotting uses matplotlib, so you can do everything that matplotlib allows you to do. Something like the following should work:

import matplotlib.pyplot as plt
fig, axes = plt.subplots()
tripcolor(p_field_00[1000], axes=axes)
axes.set_xlim(left=-1.0, right=0.0)

You can find many use cases in the notebooks: notebooks

[acse-yw11823]

many thanks!!!🙏 Sent from Outlook for iOS<https://aka.ms/o0ukef>

…

________________________________ From: ksagiyam ***@***.***> Sent: Friday, August 2, 2024 3:30:17 PM To: firedrakeproject/firedrake ***@***.***> Cc: Wu, Yifan ***@***.***>; Author ***@***.***> Subject: Re: [firedrakeproject/firedrake] Can I fd.tripcolor a part of the data in pressure field function? (Discussion #3720) This email from ***@***.*** originates from outside Imperial. Do not click on links and attachments unless you recognise the sender. If you trust the sender, add them to your safe senders list<https://spam.ic.ac.uk/SpamConsole/Senders.aspx> to disable email stamping for this address. Firedrake plotting uses matplotlib, so you can do everything that matplotlib allows you to do. Something like the following should work: import matplotlib.pyplot as plt fig, axes = plt.subplots() tripcolor(p_field_00[1000], axes=axes) axes.set_xlim(left=-1.0, right=0.0) You can find many use cases in the notebooks: notebooks<https://www.firedrakeproject.org/notebooks.html> — Reply to this email directly, view it on GitHub<#3720 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/BDENNF4FJIW26FKBV4LNYC3ZPOJXTAVCNFSM6AAAAABL4P3GAGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAMRSGQ2TMNQ>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[acse-yw11823]

I found that the output of print(p_field_1[1000].dat.data) is a 1D array. However, we can only plot the graph using fd.tripcolor. Is there any way to convert it to a 2D array so that I can plot the graph using imshow or other methods?

Many thanks!

[ksagiyam]

If I assume that p_field_1[1000] is a CGp function,

V = VectorFunctionSpace(mesh, "CG", p)
x = SpatialCoordinate(mesh)
coord = Function(V).interpolate(x)
print(coord.dat.data)

shows you the ordering of the coordinates matching p_field_1[1000].dat.data, so you can rearrange as you wish.