I’m Guru Prasad Savandaiah, pursuing my master’s thesis titled “Comparative Analysis of Computational Modeling Approaches for Flow Field Reconstruction”. In my thesis, I have to compare the models of PODGalerkin ROM, PINNs . I’m utilizing FEniCS/RBniCS for the POD Galerkin ROM implementation, specifically working on the tutorial of RBniCS/tutorials/19_navier_stokes_unsteady/tutorial_navier_stokes_unsteady_exact_1.py at master · RBniCS/RBniCS · GitHub
I was able to train the model in the offline and online modes computing the reduced solution, and then reconstruct and visualize the full order solution based on the reduced_basis functions. The code artifact is attached along with this mail.
The logic of reconstruction is
The full-order solution ( u_{\text{full}} ) is reconstructed from the reduced-order solution by combining both velocity and pressure components, with additional supremizer enrichment for stability. The The full-order solution ( u_{\text{full}} ) is reconstructed from the reduced-order solution by combining both velocity and pressure components, with additional supremizer enrichment for stability. The reconstruction process is defined as follows:
-
Full-Order Solution: The full-order solution ( u_{\text{full}} ) consists of both velocity ( u ) and pressure ( p ) components:
[
u_{\text{full}} = (u, p)
] -
Velocity Component: The velocity component ( u ) is reconstructed using the velocity POD modes ( \Phi_{\text{u}, i} ), scaled by their respective coefficients ( c_{\text{u}, i} ), along with supremizer modes ( \Phi_{\text{s}, k} ) scaled by the pressure coefficients ( c_{\text{p}, k} ):
[
u = \sum_{i} c_{\text{u}, i} \Phi_{\text{u}, i} + \sum_{k} c_{\text{p}, k} \Phi_{\text{s}, k}
] -
Pressure Component: The pressure component ( p ) is reconstructed as a sum of pressure POD modes ( \Phi_{\text{p}, j} ), each scaled by the corresponding pressure coefficient ( c_{\text{p}, j} ):
[
p = \sum_{j} c_{\text{p}, j} \Phi_{\text{p}, j}
] -
Complete Full-Order Solution: Combining both components, the final expression for the full-order solution is:
[
u_{\text{full}} = \left( \sum_{i} c_{\text{u}, i} \Phi_{\text{u}, i} + \sum_{k} c_{\text{p}, k} \Phi_{\text{s}, k}, \sum_{j} c_{\text{p}, j} \Phi_{\text{p}, j} \right)
]
The challenge is when performing POD Basis are plotted and when the reconstrcuted solution is plotted there is a unexpected blob like structure at the start of the domain, even though a parabolic profile is set for inlet with conditions of
-
self.inlet = Expression((“1./0.042025x[1](0.41 - x[1])”, “0.”), degree=2) 2) self.inlet = Expression((“6.0/0.042025x[1](0.41 - x[1])”, “0.”), degree=2)
-
self.inlet = Expression((“6.0/0.042025x[1](0.41 - x[1])”, “0.”), degree=3)
I have experimented with varying mesh = generate_mesh(domain, 64), with 64, 100, 128 mesh sizes, In all the cases the issue of blob like structure appears at the start or end when the cylinder is moved to start of the domain. The inlet condition is as default in the tutorial implementation.
Code for reconstruction:
def calculate_full_order_solution_with_supremizer_enrichment(reduced_solution, V, reduced_problem, components, debug=False):
"""
Reconstruct the full-order solution from the reduced-order solution, including supremizer enrichment.
Args:
reduced_solution: Reduced solution at a given time step.
V: Function space for the full-order solution.
reduced_problem: Reduced problem object that contains basis functions.
components: Components of the POD basis functions ('u' for velocity, 'p' for pressure, 's' for supremizers).
debug: If True, prints detailed debug information.
Returns:
full_solution: The full-order solution function containing velocity and pressure.
"""
# Define the debug print function
def debug_print(message):
if debug:
print(message)
# Initialize the full solution (velocity and pressure)
full_solution = Function(V)
print(f"Initialized full solution Function with total dimension: {V.dim()}")
# Initialize full solution vector as zero
full_solution.vector().zero()
print("Set full_solution vector to zero.")
# Get the dimension sizes for velocity and pressure components
velocity_length = V.sub(0).dim()
pressure_length = V.sub(1).dim()
total_length = velocity_length + pressure_length
print(f"Velocity component dimension: {velocity_length}, Pressure component dimension: {pressure_length}")
debug_print(f"Total full solution dimension: {total_length}")
# Create separate vectors for velocity and pressure
velocity_vector = PETScVector(MPI.comm_world, velocity_length)
pressure_vector = PETScVector(MPI.comm_world, pressure_length)
# Counter to track the index in the reduced solution vector
coefficient_index = 0
# Temporary storage for pressure coefficients for supremizer scaling
pressure_coefficients = []
# First, process the pressure components to ensure coefficients are stored before supremizers are processed
if 'p' in components:
print("--------------------------------------------------------------------------------")
print(f"Processing component: p (pressure)")
subspace = V.sub(1).collapse() # Pressure subspace
subspace_length = pressure_length
debug_print(f"Using pressure subspace with dimension {subspace.dim()}")
# Loop through the pressure POD modes
for i, pod_mode in enumerate(components['p']):
coefficient = reduced_solution.vector()[coefficient_index]
coefficient_index += 1
print(f"Coefficient for pressure mode {i}: {coefficient}")
# Store pressure coefficients for supremizer scaling
pressure_coefficients.append(coefficient)
# Create a function to store the current pressure mode
pod_function = Function(subspace)
pod_function.vector()[:] = pod_mode.vector()[velocity_length:velocity_length + pressure_length]
debug_print(f"Assigned pressure POD mode {i} to the vector with size {len(pod_function.vector())}.")
# Add to the pressure part of the full solution
pressure_vector.axpy(coefficient, pod_function.vector())
# Process the velocity components (this does not change)
if 'u' in components:
print("--------------------------------------------------------------------------------")
print(f"Processing component: u (velocity)")
subspace = V.sub(0).collapse() # Velocity subspace
subspace_length = velocity_length
debug_print(f"Using velocity subspace with dimension {subspace.dim()}")
# Loop through the velocity POD modes
for i, pod_mode in enumerate(components['u']):
coefficient = reduced_solution.vector()[coefficient_index]
coefficient_index += 1
print(f"Coefficient for velocity mode {i}: {coefficient}")
# Create a function to store the current velocity mode
pod_function = Function(subspace)
pod_function.vector()[:] = pod_mode.vector()[:velocity_length]
debug_print(f"Assigned velocity POD mode {i} to the vector with size {len(pod_function.vector())}.")
# Add to the velocity part of the full solution
velocity_vector.axpy(coefficient, pod_function.vector())
# Process the supremizer modes, ensuring the pressure coefficients are used
if 's' in components:
print("--------------------------------------------------------------------------------")
print(f"Processing component: s (supremizer)")
subspace = V.sub(0).collapse() # Supremizers affect the velocity subspace
subspace_length = velocity_length
debug_print(f"Using velocity subspace for supremizer with dimension {subspace.dim()}")
# Loop through the supremizer modes
for i, pod_mode in enumerate(components['s']):
coefficient = reduced_solution.vector()[coefficient_index]
coefficient_index += 1
print(f"Coefficient for supremizer mode {i}: {coefficient}")
# Ensure we have a corresponding pressure coefficient
if len(pressure_coefficients) > i:
pressure_coefficient = pressure_coefficients[i]
print(f"Supremizer mode {i} will be scaled by pressure coefficient: {pressure_coefficient}")
# Create a function to store the current supremizer mode
pod_function = Function(subspace)
pod_function.vector()[:] = pod_mode.vector()[:velocity_length]
debug_print(f"Assigned supremizer mode {i} to the vector with size {len(pod_function.vector())}.")
# Add to the velocity part of the full solution, scaled by the pressure coefficient
velocity_vector.axpy(pressure_coefficient, pod_function.vector())
else:
print(f"Warning: No matching pressure coefficient found for supremizer mode {i}. Coefficients list: {pressure_coefficients}")
# After summing all contributions, copy velocity and pressure vectors into full_solution
full_solution.vector().set_local(np.concatenate([velocity_vector.get_local(), pressure_vector.get_local()]))
print("Reconstructed full_solution by adding weighted POD modes and supremizers (scaled by pressure coefficients).")
print("====================================================================================")
return full_solution
# Step 1: Calculate full-order solutions for all time steps
def calculate_all_full_solutions_with_supremizer_enrichment(solution_over_time, V, reduced_problem, debug=False):
"""Calculate full-order solutions for all time steps."""
return [calculate_full_order_solution_with_supremizer_enrichment(reduced_solution, V, reduced_problem, reduced_problem.basis_functions._components, debug)
for reduced_solution in solution_over_time]
# Calculate full-order solutions for all time steps
full_solutions = calculate_all_full_solutions_with_supremizer_enrichment(solution_over_time, V, reduced_problem, debug=True)
def plot_full_order_solution(full_solution, t_idx, reduced_problem, C_x, C_y, r, debug=False):
"""
Plot the velocity magnitude and pressure field from the full-order solution
in a side-by-side layout with the cylinder filled with white.
Args:
full_solution: The full-order solution function containing velocity and pressure.
t_idx: Time index to label the plots.
reduced_problem: Reduced problem object that contains the time step information.
C_x: x-coordinate of the cylinder center.
C_y: y-coordinate of the cylinder center.
r: Radius of the cylinder.
debug: If True, prints debug information.
"""
# Define the debug print function
def debug_print(message):
if debug:
print(message)
# Split the full solution into velocity and pressure components
velocity, pressure = full_solution.split()
debug_print(f"Split full solution into velocity and pressure components.")
# Compute the velocity magnitude using UFL (Unified Form Language) operations
velocity_magnitude = sqrt(inner(velocity, velocity))
debug_print(f"Computed velocity magnitude for time step {t_idx}.")
# Project the velocity magnitude and pressure to a CG1 function space for plotting
V_mag = Function(FunctionSpace(full_solution.function_space().mesh(), 'CG', 1))
V_mag.assign(project(velocity_magnitude, V_mag.function_space()))
P_proj = Function(FunctionSpace(full_solution.function_space().mesh(), 'CG', 1))
P_proj.assign(project(pressure, P_proj.function_space()))
# Convert the FEniCS functions to NumPy arrays for plotting with Matplotlib
velocity_array = V_mag.vector().get_local()
pressure_array = P_proj.vector().get_local()
# Get the coordinates of the mesh for plotting
mesh_coordinates = V_mag.function_space().tabulate_dof_coordinates().reshape((-1, 2))
# Create a figure with side-by-side subplots (1 row, 2 columns)
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
# Plot the velocity magnitude on the left subplot
velocity_plot = ax[0].tricontourf(mesh_coordinates[:, 0], mesh_coordinates[:, 1], velocity_array, levels=50, cmap='coolwarm')
ax[0].set_title(f'Velocity Magnitude (t = {t_idx * reduced_problem.dt:.4f})')
ax[0].set_aspect('equal')
fig.colorbar(velocity_plot, ax=ax[0])
# Plot the pressure field on the right subplot
pressure_plot = ax[1].tricontourf(mesh_coordinates[:, 0], mesh_coordinates[:, 1], pressure_array, levels=50, cmap='coolwarm')
ax[1].set_title(f'Pressure (t = {t_idx * reduced_problem.dt:.4f})')
ax[1].set_aspect('equal')
fig.colorbar(pressure_plot, ax=ax[1])
# Create a white patch (filled circle) to represent the cylinder
cylinder_patch = plt.Circle((C_x, C_y), r, color='white', zorder=10)
# Add the cylinder patch (fully white) on both subplots
ax[0].add_patch(cylinder_patch) # Add to the velocity plot
ax[1].add_patch(plt.Circle((C_x, C_y), r, color='white', zorder=10)) # Add to the pressure plot
debug_print(f"Added white-filled cylinder at ({C_x}, {C_y}) with radius {r} on both plots.")
# Add a main title including the time index and time step value
plt.suptitle(f'Full-Order Solution at t = {t_idx * reduced_problem.dt:.4f} (Time Step {t_idx})', y=0.95)
# Adjust layout and show the plot
plt.tight_layout()
plt.show()
debug_print(f"Plotted velocity magnitude and pressure for time step {t_idx}.")
# Step 2: Plot full-order solutions at specified time steps
def plot_full_solutions_at_intervals(full_solutions, time_indices, reduced_problem, C_x, C_y, r, mesh=None, debug=False):
"""
Plot full-order solutions at the given time steps. If mesh is provided,
it uses a specialized plotting method with mesh details; otherwise, it uses a standard plot.
Args:
full_solutions: List of full-order solutions to be plotted.
time_indices: Indices of the time steps to be plotted (can be a list or range).
reduced_problem: The reduced problem object containing time step and other information.
C_x, C_y: Coordinates of the cylinder center for masking.
r: Radius of the cylinder.
mesh: (Optional) Mesh information. If provided, it uses mesh-aware plotting.
"""
for t_idx in time_indices:
full_solution = full_solutions[t_idx]
if mesh is not None:
# print(f"Plotting full solution at time step {t_idx} with mesh information.")
plot_full_order_solution_with_mesh(full_solution, t_idx, reduced_problem, mesh, C_x, C_y, r, debug)
else:
# print(f"Plotting full solution at time step {t_idx}")
plot_full_order_solution(full_solution, t_idx, reduced_problem, C_x, C_y, r, debug)
time_indices = range(0, len(full_solutions), 10) # Plot every 10th time step
plot_full_solutions_at_intervals(full_solutions, time_indices, reduced_problem, C_x, C_y, r)