Python Computer Language

tahanson43206 · 2024-07-27 17:57:27

It's been a while since I posted anything to this topic....

today I am exploring with ChatGPT4o how we might automate one of GW Johnson's many images he has created for his course on basic orbital mechanics. My reasoning is that if we automate the images (or one of them in this case) we might be able to reach readers who are not math experts, but who remember enough from education to be able to understand and relate to animation.

ChatGPT4o has provided a set of instructions that I will attempt to follow on one of my computers that has Python installed.

import matplotlib.pyplot as plt
import matplotlib.animation as animation
from PIL import Image
# Load the image
image_path = 'your_image_file.png'
image = Image.open(image_path)
# Convert the image to an array
image_array = np.array(image)
# Set up the figure and axis
fig, ax = plt.subplots()
im = ax.imshow(image_array)
# Function to update the animation
def update(frame):
# For now, we'll just move the image slightly to the right
ax.set_xlim(frame, frame + image_array.shape[1])
return im,
# Create the animation
ani = animation.FuncAnimation(fig, update, frames=np.arange(0, 100), blit=True)
# Show the animation
plt.show()

1) Install library

pip install matplotlib pillow

2) Create animation script

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np
from PIL import Image
# Load the image
image_path = 'your_image_file.png' # Replace with your actual image file name
image = Image.open(image_path)
# Convert the image to an array
image_array = np.array(image)
# Set up the figure and axis
fig, ax = plt.subplots()
im = ax.imshow(image_array)
# Function to update the animation
def update(frame):
# For now, we'll just move the image slightly to the right
ax.set_xlim(frame, frame + image_array.shape[1])
return im,
# Create the animation
ani = animation.FuncAnimation(fig, update, frames=np.arange(0, 100), blit=True)
# Show the animation
plt.show()

4) Run the script

python animate_image.py

The image should shift to the right in small steps

(th)

tahanson43206 · 2024-08-17 19:07:39

Python provides modules to support a web site. I'm interested in this capability because I'd like to (at least try to) provide online animation of the course material GW Johnson has developed to cover Basic Orbital Mechanics.

My initial hope was to run the software on our Azure Test site, but I've run into a snag there... the Python software wants to run on port 5000, and so far the Azure environment has steadfastly refused to allow access to anything on that port. To be more accurate, the Azure environment allows access to port 5000 ** inside ** the environment at IP address 10.0.0.4, but it refuses to allow access to the port from the outside on 40.75.112.55.

After many hours of struggle, and after enlisting Microsoft Co-Pilot as well as ChatGPT4o (which may be the same thing for all I know) I finally threw in the towel and installed the code on my Linux workstation. The installation went smoothly, and the application ran immediately. What is ** more ** the website at port 5000 was accessible to my browser, and the run produced a flood of errors. While I would have preferred to see the program work as intended, at least we have ** something ** to look at, instead of the blank screen that showed up on Azure.

It occurs to me to wonder if all those errors might have occurred on Azure, but (perhaps?) it had no way to tell me about them? That remains to be seen.

In the mean time, there are a flood of errors to correct.

If anyone is interested, here is the Python code that produced the flood of errors:

# Test 01 Python Flask web site on port 5000
from flask import Flask, render_template
import numpy as np
import plotly.graph_objs as go
from astropy.constants import G
from astropy import units as u
app = Flask(__name__)
@app.route('/orbit')
def orbit():
def get_orbit(t, r0, v0, mu):
r = np.zeros((len(t), 3))
v = np.zeros((len(t), 3))
r[0] = r0
v[0] = v0
for i in range(1, len(t)):
dt = t{i} - t[i-1]
r_mag = np.linalg.norm(r[i-1])
a = -mu / r_mag**3 * r[i-1]
v{i} = v[i-1] + a * dt
r{i} = r[i-1] + v[i-1] * dt
return r, v
r0 = np.array([7000.0, 0.0, 0.0])
v0 = np.array([0.0, 7.12, 0.0])
mu = G.to(u.km**3 / u.s**2).value * 5.972e24
t = np.linspace(0, 5400, 1000)
r, v = get_orbit(t, r0, v0, mu)
fig = go.Figure(data=[go.Scatter3d(x=r[:,0], y=r[:,1], z=r[:,2], mode='lines')])
fig.update_layout(scene=dict(
xaxis_title='X (km)',
yaxis_title='Y (km)',
zaxis_title='Z (km)'
))
return fig.to_html(full_html=False)
@app.route('/')
def index():
return render_template('index.html', orbit_plot=orbit())
if __name__ == '__main__':
app.run(debug=True)

I had to replace brackets with braces to persuade FluxBB to accept the code.

(th)

tahanson43206 · 2024-08-18 10:56:59

And here is the code that worked. The key problem had to do with how mu is calculated. The solution was to break the computation into two parts...

# Test 04 Python Flask web site on port 5000
# In Version 02 a units error is corrected
# In Version 03 we make an additional correction for mu
# In Version 04 we split calculation of mu into two steps
from flask import Flask, render_template
import numpy as np
import plotly.graph_objs as go
from astropy.constants import G
from astropy import units as u
app = Flask(__name__)
@app.route('/orbit')
def orbit():
def get_orbit(t, r0, v0, mu):
r = np.zeros((len(t), 3))
v = np.zeros((len(t), 3))
r[0] = r0
v[0] = v0
for i in range(1, len(t)):
dt = t{i} - t[i-1]
r_mag = np.linalg.norm(r[i-1])
a = -mu / r_mag**3 * r[i-1]
v{i} = v[i-1] + a * dt
r{i} = r[i-1] + v[i-1] * dt
return r, v
r0 = np.array([7000.0, 0.0, 0.0])
v0 = np.array([0.0, 7.12, 0.0])
# mu = G.to(u.km**3 / u.s**2).value * 5.972e24 wrong units
# mu = (G * (u.km**3 / u.m**3)).to(u.km**3 / u.s**2).value * 5.972e24
# mu = (G * 5.972e24 * (u.km**3 / u.m**3)).to(u.km**3 / u.s**2).value
mu = G.value * 5.972e24 # G in m³ / (kg s²) and M in kg
mu = (mu * (u.m**3 / u.s**2)).to(u.km**3 / u.s**2).value # Convert to km³ / s²
t = np.linspace(0, 5400, 1000)
r, v = get_orbit(t, r0, v0, mu)
fig = go.Figure(data=[go.Scatter3d(x=r[:,0], y=r[:,1], z=r[:,2], mode='lines')])
fig.update_layout(scene=dict(
xaxis_title='X (km)',
yaxis_title='Y (km)',
zaxis_title='Z (km)'
))
return fig.to_html(full_html=False)
@app.route('/')
def index():
return render_template('index.html', orbit_plot=orbit())
if __name__ == '__main__':
app.run(debug=True)

As before, the bracket i bracket configuration was changed to {i}

I'd like to emphasize that practically anyone should be able to run this code. Python runs on Windows, Linux, Apple and (I'm told) on Android (with Linux installed as an app).

(th)

tahanson43206 · 2024-09-07 17:49:04

We don't seem to have anyone interested in working with Python at this point, but that could change at any time.

Here is a little Python program that ChatGPT4o designed to work on the little energy storage problem under consideration in another topic...

blocked by AISE ... our webmaster (kbd512) is working on an upgrade to our FluxBB system that will not have this problem.

(th)

tahanson43206 · 2024-09-07 18:14:43

The program ran without error.

I'm unable to verify the results, but here they are:

(Note that 1 kilowatt hour is 3,600,000 joules)

Mass of water input: 1000 kg
Air pressure in Compression tank: 1000000 pascal (1 bar)
Calculate:
Total energy stored: 198100 joules or about .05 kwh
Total height: 12.55 meters
Vacuum chamber height: 1.27 meters
Compression tank height: 1.27 meters

Doubling the mass of water doubles the energy stored.
Doubling the air pressure had a lesser effect: 298100 joules, up from 198100 joules

It occurs to me the mass of the water in the vertical pipe is not being taken into account. That mass has to be lifted in order to deliver the ton of water into the receiver.

It would be interesting to study the source code, but that will not be possible until the new version of FluxBB is up and running.

(th)

tahanson43206 · 2024-09-07 19:03:54

After thinking about the behavior of the little test program a bit more, I realized that this problem is one that is a good candidate for calculus, because the force required to move an increment of water from the vacuum tank to the compression tank is not constant. The amount of force required to deal with gravity is nearly constant, but not completely, because the top of the water column rises as water flows into the receiver tank. However, the main changes are at the ends. The force required to pull an increment of water out of the vacuum tank may increase as the volume of vacuum increases although I'm not sure about that. What I ** am ** sure about is that the force required to push an increment of water into the compression tank will increase. Thus, I expect that the theoretical performance of the system will improve as the code is refined. That said, I'm not expecting dramatic improvements. We are now looking at .05 kwh to move a ton of water, and that seems likely to increase slightly with the changes proposed.

(th)

tahanson43206 · 2024-10-07 20:03:58

This post will attempt to upload a revised version of the Energy Storage calculator program that was investigated last month.

This version eliminates the gravity storage component, and concentrates on gas compression storage.

The specific scenario of business interest is oil or gas wells that are in limbo, and thus potentially available for energy storage.

Because we have not yet completed upgrade of our software, the AISE error may block this attempt.

Unbelievable! The code was accepted. Anyone with a modern computer can install Python.

This program will run under Windows. Linux (where i work) and Apple.

I'm told it will run on Android devices (such as Chromebooks) if the Linux subsystem is installed.

This preliminary version should be considered an experiment to create the framework for the study.

However, that said, it does run without obvious errors, so the output may be correct. I certainly don't guarantee it.

# Compressed Gas Energy Storage Calculator - Version 3
# Focus: Evaluating energy storage potential in oil and gas wells
import math
import tkinter as tk
from tkinter import messagebox
# Function to calculate work of isothermal and adiabatic compression
def work_compression(p_initial, v_initial, v_final, gamma=1.4, isothermal=True):
if isothermal:
# Isothermal compression work: W = P_initial * V_initial * ln(V_initial / V_final)
return p_initial * v_initial * math.log(v_initial / v_final)
else:
# Adiabatic compression work: W = (P_initial * V_initial * [(V_final/V_initial)^(1-gamma) - 1]) / (gamma - 1)
return (p_initial * v_initial * ((v_initial / v_final) ** (gamma - 1) - 1)) / (1 - gamma)
# Main function to perform the calculation
# Main function to perform the calculation
def calculate():
try:
# Retrieve input values from the GUI
P_initial = float(entry_P_initial.get()) * 1000 # Initial pressure (kPa to Pa)
P_final = float(entry_P_final.get()) * 1000 # Final pressure (kPa to Pa)
pipe_diameter = float(entry_pipe_diameter.get()) # Diameter of riser pipe (m)
pipe_length = float(entry_pipe_length.get()) # Height of riser pipe (m)
compression_type = compression_mode.get() # Compression type: Isothermal or Adiabatic
# Derived values
pipe_area = math.pi * (pipe_diameter / 2) ** 2 # Cross-sectional area of the pipe (m²)
pipe_volume = pipe_area * pipe_length # Volume of the pipe (m³)
# Choose the appropriate calculation method
isothermal = (compression_type == "Isothermal")
gamma = 1.4 # Specific heat ratio (default for air)
# Calculate the compression work based on selected mode
compression_work = work_compression(P_initial, pipe_volume, pipe_volume / (P_final / P_initial), gamma, isothermal)
# Convert compression work to kilojoules (kJ)
compression_work_kilojoules = compression_work / 1000
# Display results in the output labels
output_compression_kj.config(text=f"{compression_work_kilojoules:.2f} kJ")
except ValueError:
messagebox.showerror("Input Error", "Please enter valid numeric values.")
# Set up the Tkinter GUI# Main function to perform the calculation
# revised calculate function for Version 03
def calculate():
try:
# Retrieve input values from the GUI
P_initial = float(entry_P_initial.get()) * 1000 # Initial pressure (kPa to Pa)
P_final = float(entry_P_final.get()) * 1000 # Final pressure (kPa to Pa)
pipe_diameter = float(entry_pipe_diameter.get()) # Diameter of riser pipe (m)
pipe_length = float(entry_pipe_length.get()) # Height of riser pipe (m)
compression_type = compression_mode.get() # Compression type: Isothermal or Adiabatic
# Derived values
pipe_area = math.pi * (pipe_diameter / 2) ** 2 # Cross-sectional area of the pipe (m²)
pipe_volume = pipe_area * pipe_length # Volume of the pipe (m³)
# Choose the appropriate calculation method
isothermal = (compression_type == "Isothermal")
gamma = 1.4 # Specific heat ratio (default for air)
# Calculate the compression work based on selected mode
compression_work = work_compression(P_initial, pipe_volume, pipe_volume / (P_final / P_initial), gamma, isothermal)
# Convert compression work to kilojoules (kJ)
compression_work_kilojoules = compression_work / 1000
# Display results in the output labels
output_compression_kj.config(text=f"{compression_work_kilojoules:.2f} kJ")
except ValueError:
messagebox.showerror("Input Error", "Please enter valid numeric values.")
root = tk.Tk()
root.title("Compressed Gas Energy Storage Calculator - Version 3")
# Labels and input fields
tk.Label(root, text="Initial Pressure (kPa):").grid(row=0, column=0, padx=10, pady=5)
entry_P_initial = tk.Entry(root)
entry_P_initial.grid(row=0, column=1, padx=10, pady=5)
entry_P_initial.insert(0, "101.325")
tk.Label(root, text="Final Pressure (kPa):").grid(row=1, column=0, padx=10, pady=5)
entry_P_final = tk.Entry(root)
entry_P_final.grid(row=1, column=1, padx=10, pady=5)
entry_P_final.insert(0, "202.65")
tk.Label(root, text="Riser Pipe Diameter (m):").grid(row=2, column=0, padx=10, pady=5)
entry_pipe_diameter = tk.Entry(root)
entry_pipe_diameter.grid(row=2, column=1, padx=10, pady=5)
entry_pipe_diameter.insert(0, "0.04")
tk.Label(root, text="Riser Pipe Height (m):").grid(row=3, column=0, padx=10, pady=5)
entry_pipe_length = tk.Entry(root)
entry_pipe_length.grid(row=3, column=1, padx=10, pady=5)
entry_pipe_length.insert(0, "10")
# Compression type selection
compression_mode = tk.StringVar(value="Isothermal")
tk.Label(root, text="Compression Type:").grid(row=4, column=0, padx=10, pady=5)
tk.Radiobutton(root, text="Isothermal", variable=compression_mode, value="Isothermal").grid(row=4, column=1, padx=10, pady=5)
tk.Radiobutton(root, text="Adiabatic", variable=compression_mode, value="Adiabatic").grid(row=4, column=2, padx=10, pady=5)
# above is first half
# Output labels for results
tk.Label(root, text="Compression Work (kJ):").grid(row=5, column=0, padx=10, pady=5)
output_compression_kj = tk.Label(root, text="0.00 kJ")
output_compression_kj.grid(row=5, column=1, padx=10, pady=5)
# Button to trigger calculation
calc_button = tk.Button(root, text="Calculate", command=calculate)
calc_button.grid(row=6, column=1, columnspan=2, pady=10)
# Start the Tkinter event loop
root.mainloop()

(th)

tahanson43206 · 2024-11-09 14:45:24

Work on developing a browser accessible Python application has been ongoing for some time. There should be reports somewhere in the archive, but I can't remember where I put them. This topic is as good as any for a short progress report.

Today I worked with ChatGPT4o to update the previous version of pythontest05.py, which is available for viewing at:
http://40.75.112.55/flask or http://40.75.112.55/flask2

The new version is running on my development PC as pythontest06.py. It uses an html file: index06.html.

The new version has input fields for X, Y and Z, and a recalculate button.

I don't think the new version is ready for public evaluation.

A detail that may be of interest to someone in the readership is that the Linux "touch" command turned out to be useful for activating new code without restarting the Apache server. In the case of a Flask application such as this one, the file that we "touch" is flask.wsgi.

There may be something comparable in a php application, but I don't know that right now.

Today's concentration was getting the new updates to the Flask application to work.

(th)

tahanson43206 · 2024-11-13 23:17:46

ChatGPT4o and I are working on pythontest08.py. The program visible on the Azure web site is still the simple little starter program we used to get Flask running as an online service from Azure. Development is progressing on the local development PC.

This evening we hammered out a series of small steps to advance toward the goal of supporting GW Johnson's course on Basic Orbital mechanics.

For one thing, we now have a little blue Earth in the center of the 3D graph, and the orbit now closes properly. We've added a button to allow the student to try different starting positions for the satellite.

Enhancements to come include advancing from a simple circle to a true ellipse, which (after all) is the focus of the first lesson in GW's course.

(th)

tahanson43206 · 2024-12-14 16:28:57

It's been a while since we had something to offer in this topic.

Today, ChatGPT4o showed off some new features, and "we" created a program to compute color changes for a 3D printer project.
The program is not yet perfect, but it does work, and it will run on any computer that has Python3 and tkinter for the GUI.

I'll give the new FluxBB software a severe test now, by attempting to paste the source code.
Amazing !!!! It worked !!! Bravo to kbd512 and everyone else who worked on the upgrade!

#20241214-3DcolorsV06.py
import tkinter as tk
from tkinter import filedialog, messagebox
import itertools
# Function to generate combinations
def generate_combinations(colors, timing, model_count):
schedule = []
for model in range(1, model_count + 1):
stage = 1
for combo in itertools.permutations(colors):
schedule.append((f"Model {model}: Stage {stage}", " -> ".join(combo), timing))
stage += 1
return schedule
# Function to save schedule to file
def save_to_file(schedule):
file_path = filedialog.asksaveasfilename(defaultextension=".txt",
filetypes=[("Text Files", "*.txt"), ("All Files", "*.*")])
if not file_path:
return
with open(file_path, "w") as file:
for model_stage, combination, time in schedule:
file.write(f"{model_stage}: {combination} - {time} minutes\n")
messagebox.showinfo("Success", "Schedule saved successfully!")
# Function to handle generate button click
def on_generate():
global current_schedule # Make schedule globally accessible
try:
colors = [color.strip() for color in color_entry.get().split(",")]
if not colors or any(c == "" for c in colors):
raise ValueError("Colors: Please enter a comma-separated list of color names.")
try:
timing = int(timing_var.get())
if timing <= 0:
raise ValueError("Timing: Please enter a positive integer.")
except ValueError:
raise ValueError("Timing: Please enter a positive integer.")
try:
model_count = int(model_var.get())
if model_count <= 0:
raise ValueError("Models: Please enter a positive integer.")
except ValueError:
raise ValueError("Models: Please enter a positive integer.")
current_schedule = generate_combinations(colors, timing, model_count)
# Display schedule
output_text.delete("1.0", tk.END)
for model_stage, combination, time in current_schedule:
output_text.insert(tk.END, f"{model_stage}: {combination} - {time} minutes\n")
# Enable save button
save_button.config(state=tk.NORMAL)
except ValueError as e:
messagebox.showerror("Input Error", str(e))
# Create main Tkinter window
root = tk.Tk()
root.title("3D Printing Schedule Planner")
# Input frame
input_frame = tk.Frame(root)
input_frame.pack(pady=10)
# Colors input
color_label = tk.Label(input_frame, text="Enter Colors (comma-separated):")
color_label.grid(row=0, column=0, padx=5, pady=5)
color_entry = tk.Entry(input_frame, width=30)
color_entry.insert(0, "Red,Green,Blue,Yellow,White")
color_entry.grid(row=0, column=1, padx=5, pady=5)
# Timing input
timing_label = tk.Label(input_frame, text="Time per Stage (minutes):")
timing_label.grid(row=1, column=0, padx=5, pady=5)
timing_var = tk.StringVar(value="10")
timing_entry = tk.Entry(input_frame, textvariable=timing_var, width=5)
timing_entry.grid(row=1, column=1, padx=5, pady=5)
# Model count input
model_label = tk.Label(input_frame, text="Number of Models:")
model_label.grid(row=2, column=0, padx=5, pady=5)
model_var = tk.StringVar(value="1")
model_entry = tk.Entry(input_frame, textvariable=model_var, width=5)
model_entry.grid(row=2, column=1, padx=5, pady=5)
# Buttons
button_frame = tk.Frame(root)
button_frame.pack(pady=10)
generate_button = tk.Button(button_frame, text="Generate Schedule", command=on_generate)
generate_button.grid(row=0, column=0, padx=10)
save_button = tk.Button(button_frame, text="Save to File", command=lambda: save_to_file(current_schedule), state=tk.DISABLED)
save_button.grid(row=0, column=1, padx=10)
# Output display
output_frame = tk.Frame(root)
output_frame.pack(pady=10)
output_label = tk.Label(output_frame, text="Generated Schedule:")
output_label.pack(anchor="w")
output_text = tk.Text(output_frame, width=50, height=15)
output_text.pack()
# Run the Tkinter event loop
root.mainloop()

(th)

tahanson43206 · 2025-02-25 20:46:05

This will be a test of the updated FluxBB software we are running...

The quote below will be the Version #3 of a program to show animation of a payload on a tether from a rotating object.

The program should run on any modern computer with Windows, Linux, Android (running Linux), or Apple.

# Version 03 - Asteroid Tether Animation
import tkinter as tk
import math
# Constants
CANVAS_SIZE = 600 # Size of Tkinter canvas
ASTEROID_RADIUS = 80 # Radius of the asteroid (scaled for visualization)
TETHER_LENGTH = 300 # Max tether length (scaled down for visualization)
FRAME_RATE = 10 # Milliseconds per frame update
def reset_simulation():
global payload_x, payload_y, payload_vx, payload_vy, tether_angle, running, rotation_speed
payload_x = CANVAS_SIZE // 2 + ASTEROID_RADIUS # Start on asteroid surface
payload_y = CANVAS_SIZE // 2
payload_vx = 1.5 # Initial outward velocity (scaled for visualization)
payload_vy = 0
tether_angle = 0
running = False # Control simulation start/stop
try:
rpm = float(rotation_speed_entry.get())
rotation_speed = (rpm * 2 * math.pi) / 60 # Convert RPM to radians per second
except ValueError:
rotation_speed = (1 * 2 * math.pi) / 60 # Default to 1 RPM if input is invalid
def update():
global payload_x, payload_y, payload_vx, payload_vy, tether_angle, running

if not running:
return

# Compute new position (simple outward motion for now)
payload_x += payload_vx
payload_y += payload_vy

# Compute new tether angle as asteroid rotates
tether_angle += rotation_speed * (FRAME_RATE / 1000) # Adjust for frame timing

# Check if tether reaches max length
dx = payload_x - CANVAS_SIZE // 2
dy = payload_y - CANVAS_SIZE // 2
distance = math.sqrt(dx**2 + dy**2)

if distance >= TETHER_LENGTH:
# Constrain payload to max tether length
angle = math.atan2(dy, dx)
payload_x = CANVAS_SIZE // 2 + TETHER_LENGTH * math.cos(angle)
payload_y = CANVAS_SIZE // 2 + TETHER_LENGTH * math.sin(angle)

# Stop outward velocity when tether goes taut
payload_vx = 0
payload_vy = 0

# Update canvas
canvas.delete("all")

# Draw asteroid
canvas.create_oval(CANVAS_SIZE//2 - ASTEROID_RADIUS, CANVAS_SIZE//2 - ASTEROID_RADIUS,
CANVAS_SIZE//2 + ASTEROID_RADIUS, CANVAS_SIZE//2 + ASTEROID_RADIUS,
fill='gray', outline='black')

# Draw payload
canvas.create_oval(payload_x - 5, payload_y - 5, payload_x + 5, payload_y + 5, fill='red')

# Draw tether
canvas.create_line(CANVAS_SIZE // 2, CANVAS_SIZE // 2, payload_x, payload_y, fill='black')

root.after(FRAME_RATE, update) # Schedule next frame
def start_simulation():
global running
reset_simulation()
running = True
update()
# Create Tkinter window
root = tk.Tk()
root.title("Asteroid Tether Animation")
canvas = tk.Canvas(root, width=CANVAS_SIZE, height=CANVAS_SIZE, bg="white")
canvas.pack()
# Input for Rotation Speed
rotation_speed_label = tk.Label(root, text="Rotation Speed RPM")
rotation_speed_label.pack(anchor="w")
rotation_speed_entry = tk.Entry(root)
rotation_speed_entry.insert(0, "1") # Default value 1 RPM
rotation_speed_entry.pack(anchor="w")
# Start/Restart button
start_button = tk.Button(root, text="Start/Restart Simulation", command=start_simulation)
start_button.pack()
reset_simulation()
root.mainloop()

(th)

tahanson43206 · 2025-02-26 10:44:10

This is an experimental post...

Apparently ChatGPT4o can share a program in development:

https://chatgpt.com/canvas/shared/67bf4 … 2309dc3322

The version in progress at the time I created the link above is #13 ... we are working on animation of the tether which is particularly difficult.

Update: The link definitely points to version 13. I have completed a number of updates since, and the link above still pulls up 13.

We decided to implement a physics engine called "pymunk" and the latest version has plenty of dynamic activity. The payload is bouncing wildly against the asteroid which means the physics engine is treating the asteroid as a solid and implementing collision detection.

The tether is being rendered as a kinky spring instead of a straight line as I'd expect.

Update: the "kinky spring" rendering was due to a setting inside the program. Segments of the tether had been set to 20, to reduce processing time. By setting the segments count to 50, the behavior of the tether improved.

Once we got the physics engine going, I discovered that the scenario I'd set up was going to behave differently than I'd expected.

The payload is bouncing against the asteroid. Since there are currently no provisions for damping, and since the tether is considered as elastic, the wild banging around of the payload makes sense. In the Real Universe, if an astronaut tossed a payload out away from the surface and trailed it with a cord, I'd (NOW) expect to see something like this simulation.

(th)

tahanson43206 · 2025-04-10 07:38:21

This post is reserved for a display of Python to compute temperature rise in a (hypothetical) Solar Powered Space Vessel)

It should be noted that the design in study is a direct result of an idea originated with kbd512, for a way to update 1950's concepts using modern technology.

The program listed below should run on any modern desktop, under Windows, Linux, Apple or Android (with the Linux extension)

# HeatPipeSimulator01 2025/04/09 by ChatGPT4o
import tkinter as tk
from tkinter import ttk, scrolledtext
# Constants for hydrogen
DEFAULT_CP = 14300 # J/kg-K, specific heat capacity at high temp
class HeatPipeSimulator:
def __init__(self, root):
self.root = root
self.root.title("Heat Pipe Simulator")
self.setup_widgets()
def setup_widgets(self):
# Frame for inputs
input_frame = ttk.Frame(self.root, padding="10")
input_frame.grid(row=0, column=0, sticky="w")
# Input labels and entries
labels = ["Pipe Length (m)", "Segment Length (m)", "Power per Segment (W)",
"Mass Flow Rate (kg/s)", "Specific Heat Capacity (J/kg-K)",
"Initial Temperature (K)", "Efficiency (%)"]
defaults = [240, 0.038, 6325, 1.0, DEFAULT_CP, 300, 100]
self.entries = []
for i, (label, default) in enumerate(zip(labels, defaults)):
ttk.Label(input_frame, text=label).grid(row=i, column=0, sticky="w")
entry = ttk.Entry(input_frame)
entry.insert(0, str(default))
entry.grid(row=i, column=1, sticky="w")
self.entries.append(entry)
# Run button
ttk.Button(input_frame, text="Run Simulation", command=self.run_simulation).grid(row=len(labels), column=0, columnspan=2, pady=10)
# Output area
self.output = scrolledtext.ScrolledText(self.root, width=80, height=25, wrap=tk.WORD)
self.output.grid(row=1, column=0, padx=10, pady=10)
def run_simulation(self):
# Get input values
try:
pipe_length = float(self.entries[0].get())
segment_length = float(self.entries[1].get())
power_per_segment = float(self.entries[2].get())
mass_flow = float(self.entries[3].get())
cp = float(self.entries[4].get())
initial_temp = float(self.entries[5].get())
efficiency = float(self.entries[6].get()) / 100.0
except ValueError:
self.output.insert(tk.END, "\nError: Please enter valid numeric values.\n")
return
# Simulation
num_segments = int(pipe_length / segment_length)
temperature = initial_temp
total_energy = 0
self.output.delete(1.0, tk.END)
self.output.insert(tk.END, f"Running simulation with {num_segments} segments...\n\n")
for i in range(num_segments):
energy = power_per_segment * efficiency
delta_T = energy / (mass_flow * cp)
temperature += delta_T
total_energy += energy
if i % 100 == 0 or i == num_segments - 1:
position = (i + 1) * segment_length
self.output.insert(tk.END, f"Position: {position:.2f} m, Temperature: {temperature:.2f} K\n")
self.output.insert(tk.END, f"\nFinal Temperature: {temperature:.2f} K\n")
self.output.insert(tk.END, f"Total Energy Input: {total_energy / 1e6:.2f} MJ\n")
if temperature > 4000:
self.output.insert(tk.END, "\nWARNING: Final temperature exceeds 4000 K — material or hydrogen dissociation limits may be exceeded!\n")
if __name__ == '__main__':
root = tk.Tk()
app = HeatPipeSimulator(root)
root.mainloop()

Note: to change insolation, you can update the watts per segment value. The default value is 1 watt per optical device.

This value is based upon a conservative estimate of insolation of 1000 watts per square meter.

To change the value to show performance at various other levels of insolation, just multiply the default number by a factor.

For example, the factor for LEO would be 1.361 (1 point 361)

The factor for Mars would be: .5862 (point 5862)

In contrast, the factor for Venus would be: 2600

Please note that the system as designed cannot handle the insolation available at Venus. One solution is to simply cover optical receivers as needed to keep power levels in reasonable limits.

Reminder: to run this program on a desktop, save the file as (your-preferred-name).py

Then enter python3 your-preffered-name.py

If you have never run Python on your machine, you'll need to install some components, but that should go quickly.

NewMars members who would like to try this are welcome to report any difficulties.

Hopefully assistance will be forthcoming.

(th)

tahanson43206 · 2025-04-10 20:50:19

This post is intended to provide the source code for Version 7 of a Python program to compute heating using solar power.

This version includes an option to introduce hydrogen directly as liquid, or to run the liquid through the nozzle and up a long pipe to the top of the heat pipe.

I held the program back until I had a chance to run some trials. It seems to be working but I'd appreciate someone else helping with evaluation of the performance of the program.

# Heat Pipe Simulator Version 07 2025/04/10 by ChatGPT4o
import tkinter as tk
from tkinter import ttk, scrolledtext
# Constants for hydrogen
DEFAULT_CP = 14300 # J/kg-K, specific heat capacity at high temp
LATENT_HEAT_VAP_H2 = 445_000 # J/kg
TEMP_BOIL_H2 = 20.3 # K
TEMP_IN_LH2 = 20 # K
TEMP_IN_GAS = TEMP_BOIL_H2 # K
class HeatPipeSimulator:
def __init__(self, root):
self.root = root
self.root.title("Heat Pipe Simulator Version 7")
self.setup_widgets()
def setup_widgets(self):
# Frame for inputs
input_frame = ttk.Frame(self.root, padding="10")
input_frame.grid(row=0, column=0, sticky="w")
# Input labels and entries
labels = ["Pipe Length (m)", "Segment Length (m)", "Power per Segment (W)",
"Mass Flow Rate (kg/s)", "Specific Heat Capacity (J/kg-K)",
"Initial Temperature (K)", "Efficiency (%)"]
defaults = [240, 0.038, 6325, 1.0, DEFAULT_CP, 20.3, 100]
self.entries = []
for i, (label, default) in enumerate(zip(labels, defaults)):
ttk.Label(input_frame, text=label).grid(row=i, column=0, sticky="w")
entry = ttk.Entry(input_frame)
entry.insert(0, str(default))
entry.grid(row=i, column=1, sticky="w")
self.entries.append(entry)
# Mode selector
ttk.Label(input_frame, text="Hydrogen Entry Mode:").grid(row=len(labels), column=0, sticky="w", pady=(10, 0))
self.mode_var = tk.StringVar(value="LH2")
ttk.Radiobutton(input_frame, text="LH2 (Liquid Hydrogen)", variable=self.mode_var, value="LH2", command=self.update_initial_temp).grid(row=len(labels)+1, column=0, sticky="w")
ttk.Radiobutton(input_frame, text="Gas (Pre-heated)", variable=self.mode_var, value="Gas", command=self.update_initial_temp).grid(row=len(labels)+2, column=0, sticky="w")
# Run button
ttk.Button(input_frame, text="Run Simulation", command=self.run_simulation).grid(row=len(labels)+3, column=0, columnspan=2, pady=10)
# Output area
self.output = scrolledtext.ScrolledText(self.root, width=80, height=25, wrap=tk.WORD)
self.output.grid(row=1, column=0, padx=10, pady=10)
def update_initial_temp(self):
mode = self.mode_var.get()
initial_temp_entry = self.entries[5]
if mode == "LH2":
initial_temp_entry.delete(0, tk.END)
initial_temp_entry.insert(0, str(TEMP_IN_LH2))
else:
initial_temp_entry.delete(0, tk.END)
initial_temp_entry.insert(0, str(TEMP_IN_GAS))
def run_simulation(self):
# Get input values
try:
pipe_length = float(self.entries[0].get())
segment_length = float(self.entries[1].get())
power_per_segment = float(self.entries[2].get())
mass_flow = float(self.entries[3].get())
cp = float(self.entries[4].get())
initial_temp_input = float(self.entries[5].get())
efficiency = float(self.entries[6].get()) / 100.0
mode = self.mode_var.get()
except ValueError:
self.output.insert(tk.END, "\nError: Please enter valid numeric values.\n")
return
# Adjust initial temperature and energy budget if LH2 selected
if mode == "LH2":
energy_vaporize = mass_flow * LATENT_HEAT_VAP_H2
energy_preheat = energy_vaporize
temperature = TEMP_BOIL_H2
else:
energy_preheat = 0
temperature = initial_temp_input
# Simulation
num_segments = int(pipe_length / segment_length)
total_energy = -energy_preheat
self.output.delete(1.0, tk.END)
self.output.insert(tk.END, f"Running simulation with {num_segments} segments...\n")
self.output.insert(tk.END, f"Mode: {mode} (Hydrogen enters as {'Liquid' if mode == 'LH2' else 'Gas'})\n\n")
if mode == "LH2":
self.output.insert(tk.END, f"Energy to Vaporize: {energy_vaporize / 1e6:.2f} MJ\n")
self.output.insert(tk.END, f"Total Preheat Energy: {energy_preheat / 1e6:.2f} MJ\n")
for i in range(num_segments):
energy = power_per_segment * efficiency
delta_T = energy / (mass_flow * cp)
temperature += delta_T
total_energy += energy
if i % 100 == 0 or i == num_segments - 1:
position = (i + 1) * segment_length
self.output.insert(tk.END, f"Position: {position:.2f} m, Temperature: {temperature:.2f} K\n")
self.output.insert(tk.END, f"\nFinal Temperature: {temperature:.2f} K\n")
self.output.insert(tk.END, f"Total Net Energy Input: {total_energy / 1e6:.2f} MJ\n")
if temperature > 4000:
self.output.insert(tk.END, "\nWARNING: Final temperature exceeds 4000 K — material or hydrogen dissociation limits may be exceeded!\n")
if __name__ == '__main__':
root = tk.Tk()
app = HeatPipeSimulator(root)
root.mainloop()

A detail that i don't fully understand is how to estimate the temperature of the gas after it is vaporized by the preheating operation.

(th)

tahanson43206 · 2025-04-23 20:19:16

This post provides the Python code for a utility to read plot data points and show the points in the X-Y format.

The program expects a simple CSV file format with three columns of data

# 20250423OpenFOAMPlot Revised program to display points X Y and Z on a plot
# Revision 06 to close cleanly 07 to fix time slider
import tkinter as tk
from tkinter import ttk, filedialog
import matplotlib.pyplot as plt
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
import csv
class MeshViewer:
def __init__(self, root):
self.root = root
self.root.title("MeshViewer: Merlin Edition V07")
self.points = []
self.index = 0
self.auto_play_id = None
self.create_widgets()
# Handle window close event
self.root.protocol("WM_DELETE_WINDOW", self.close)
def create_widgets(self):
# Control panel
control_frame = ttk.Frame(self.root)
control_frame.pack(side=tk.TOP, fill=tk.X)
ttk.Button(control_frame, text="Load CSV", command=self.load_csv).pack(side=tk.LEFT)
ttk.Button(control_frame, text="Next Point", command=self.plot_next_point).pack(side=tk.LEFT)
ttk.Button(control_frame, text="Auto Play", command=self.auto_play).pack(side=tk.LEFT)
ttk.Button(control_frame, text="Reset", command=self.reset).pack(side=tk.LEFT)
self.status_label = ttk.Label(control_frame, text="No file loaded.")
self.status_label.pack(side=tk.LEFT, padx=10)
# Speed control slider
speed_frame = ttk.Frame(self.root)
speed_frame.pack(side=tk.TOP, fill=tk.X)
self.speed_label = ttk.Label(speed_frame, text="Speed: 3.57 ms/point") # Default
self.speed_label.pack(side=tk.LEFT)
self.speed_slider = ttk.Scale(speed_frame, from_=2, to=100, orient=tk.HORIZONTAL, command=self.update_speed_label)
self.speed_slider.set(3.57) # default ~28,000 points / 60 sec
self.speed_slider.pack(side=tk.LEFT, fill=tk.X, expand=True, padx=10)
# Plotting area
self.fig, self.ax = plt.subplots()
self.ax.set_title("Mesh Point Animation")
self.ax.set_xlabel("X")
self.ax.set_ylabel("Y")
self.canvas = FigureCanvasTkAgg(self.fig, master=self.root)
self.canvas.get_tk_widget().pack(fill=tk.BOTH, expand=True)
ttk.Button(control_frame, text="Close", command=self.close).pack(side=tk.RIGHT)
def update_speed_label(self, val):
ms = float(val)
self.speed_label.config(text=f"Speed: {ms:.2f} ms/point")
def close(self):
if self.auto_play_id is not None:
try:
self.root.after_cancel(self.auto_play_id)
except Exception:
pass
self.auto_play_id = None
# Explicitly tell tkinter to break out of the mainloop
try:
self.root.quit()
except Exception:
pass
self.root.destroy()
def load_csv(self):
file_path = filedialog.askopenfilename(filetypes=[("CSV Files", "*.csv")])
if not file_path:
return
with open(file_path, newline='') as csvfile:
reader = csv.reader(csvfile, delimiter=' ', skipinitialspace=True)
self.points = [(float(row[0]), float(row[1])) for row in reader if len(row) >= 2]
self.status_label.config(text=f"{len(self.points)} points loaded.")
self.reset()
def plot_next_point(self):
if self.index < len(self.points):
x, y = self.points[self.index]
self.ax.plot(x, y, 'bo', markersize=2)
self.canvas.draw()
self.index += 1
def auto_play(self):
if self.index < len(self.points):
self.plot_next_point()
delay = max(1, int(float(self.speed_slider.get())))
self.auto_play_id = self.root.after(delay, self.auto_play)
else:
self.auto_play_id = None
def reset(self):
if self.auto_play_id:
self.root.after_cancel(self.auto_play_id)
self.auto_play_id = None
self.index = 0
self.ax.clear()
self.ax.set_title("Mesh Point Animation")
self.ax.set_xlabel("X")
self.ax.set_ylabel("Y")
self.canvas.draw()
# Run the application
root = tk.Tk()
app = MeshViewer(root)
root.mainloop()
print("Application closed cleanly.")

And here is some data for testing:

-0.3770523164 0.1838259122 0.06303653966
-0.3541049515 0.1836518267 0.06303653966
-0.3311599763 0.1833424661 0.06303653966
-0.3082172508 0.182989529 0.06303653966

Update 2025/04/24 ... I just installed the program and ran it on a Windows 10 machine.

Python3 was already installed and some features added. However, matplotlib was missing.

The command to add it is: pip install matplotlib

If you want to run this program and don't have Python installed, ask Gemini, or look at the instructions in this post:
https://newmars.com/forums/viewtopic.ph … 92#p212392

And! When you save the program with Wordpad, use Text MSDOS format, and use the file type of "py"

When you save the data using Wordpad, use Text MSDOS format, and use the file type of "csv"

The program works well, and the data set show four points. You can load them a point at a time, or all at once.

To create your own set of points to plot, just use Excel to generate a set of points.

If anyone is interested in seeing the complete set of 28,000+ points generated for the Merlin Engine OpenFOAM model, it is up on Dropbox. The display shows clearly how the flow lines were generated. See below image for the link.

Update: This is a view of the points in the Merlin nozzle simulation as the right side of the field is being added.

Update:

https://www.dropbox.com/scl/fi/ekvgtj0i … uowjl&dl=0

The file at the link above contains points that define the structure created by our anonymous donor, of a working model of Merlin Engine for flow analysis.

(th)

tahanson43206 · 2025-04-27 07:21:53

This post provides the Python code that generates a plot of Courant Number changes during a run of OpenFOAM against a model of Merlin Engine created by an anonymous third party.

# 2025/04/26 Version 02 of OpenFOAM log viewer
# This viewer reads a log, extracts Courant Number data and plots it
import re
import tkinter as tk
from tkinter import filedialog
import matplotlib.pyplot as plt

def main():
    # Start Tkinter root (hide the main window)
    root = tk.Tk()
    root.withdraw()  # Hide the blank main window

    # Open file dialog
    file_path = filedialog.askopenfilename(
        title="Select OpenFOAM log file",
        filetypes=[("Log files", "*.log"), ("All files", "*.*")]
    )

    if not file_path:
        print("No file selected. Exiting.")
        return

    # Prepare lists
    times = []
    courant_means = []
    courant_maxes = []

    current_time = None  # Track latest time

    # Read and parse the selected log file
    with open(file_path, "r") as file:
        for line in file:
            time_match = re.search(r'Time = ([\deE\+\-\.]+)', line)
            if time_match:
                current_time = float(time_match.group(1))

            courant_match = re.search(r'Courant Number mean: ([\d\.eE\+\-]+) max: ([\d\.eE\+\-]+)', line)
            if courant_match and current_time is not None:
                courant_means.append(float(courant_match.group(1)))
                courant_maxes.append(float(courant_match.group(2)))
                times.append(current_time)

    if not times:
        print("No Courant Number entries found in the file!")
        return

    # Plotting
    plt.figure(figsize=(10,5))
    plt.plot(times, courant_means, label="Mean Courant Number", marker="o")
    plt.plot(times, courant_maxes, label="Max Courant Number", marker="x")
    plt.xlabel("Time [s]")
    plt.ylabel("Courant Number")
    plt.title("Courant Number vs Time")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    main()

Here is the output:

(th)

tahanson43206 · 2025-04-28 13:27:01

This program is the second of a series I am hoping to create for analysis of an OpenFOAM run against a model of the Merlin Engine created several years ago by an anonymous contributor. This on examines the Velocity of incoming gas in the X direction. It shows a very stable value across the 281 time steps recorded before the simulation crashed at the throat of the nozzle.

# 2025/04/28 Python program to read OpenFOAM run log, to study Ux residuals
import re
import tkinter as tk
from tkinter import filedialog
import matplotlib.pyplot as plt

def main():
    # Start Tkinter root (hide the main window)
    root = tk.Tk()
    root.withdraw()  # Hide the blank main window

    # Open file dialog
    file_path = filedialog.askopenfilename(
        title="Select OpenFOAM log file",
        filetypes=[("Log files", "*.log"), ("All files", "*.*")]
    )

    if not file_path:
        print("No file selected. Exiting.")
        return

    # Prepare lists
    times = []
    ux_initial_residuals = []

    current_time = None  # Track latest time

    # Read and parse the selected log file
    with open(file_path, "r") as file:
        for line in file:
            # Match Time lines
            time_match = re.search(r'Time = ([\deE\+\-\.]+)', line)
            if time_match:
                current_time = float(time_match.group(1))

            # Match Ux residual lines
            ux_match = re.search(r'smoothSolver: *Solving for Ux, *Initial residual = ([\d\.eE\+\-]+)', line)
            if ux_match and current_time is not None:
                ux_initial_residuals.append(float(ux_match.group(1)))
                times.append(current_time)

    if not times:
        print("No Ux residual entries found in the file!")
        return

    # Plotting
    plt.figure(figsize=(10,5))
    plt.semilogy(times, ux_initial_residuals, label="Initial Residual (Ux)", marker="o")
    plt.xlabel("Time [s]")
    plt.ylabel("Initial Residual (log scale)")
    plt.title("Ux Initial Residual vs Time")
    plt.legend()
    plt.grid(True, which="both", linestyle="--", linewidth=0.5)

    # Optional Save Button (Save as PNG)
    def save_plot():
        save_path = filedialog.asksaveasfilename(defaultextension=".png",
                                                  filetypes=[("PNG files", "*.png"), ("All files", "*.*")])
        if save_path:
            plt.savefig(save_path)
            print(f"Plot saved as {save_path}")

    # Add Save Button
    save_button_ax = plt.axes([0.81, 0.01, 0.1, 0.05])
    save_button = plt.Button(save_button_ax, "Save Plot")
    save_button.on_clicked(lambda event: save_plot())

    plt.tight_layout(rect=[0, 0.05, 1, 1])  # Leave room for save button
    plt.show()

if __name__ == "__main__":
    main()

fGI1vvg_d.webp?maxwidth=256&shape=square
(th)

tahanson43206 · 2025-04-28 19:47:34

This post is intended to hold the code for a log viewer to report on Uy residuals.

# 20250428LogViewerUy.py
import re
import tkinter as tk
from tkinter import filedialog
import matplotlib.pyplot as plt

def main():
    # Start Tkinter root (hide the main window)
    root = tk.Tk()
    root.withdraw()  # Hide the blank main window

    # Open file dialog
    file_path = filedialog.askopenfilename(
        title="Select OpenFOAM log file",
        filetypes=[("Log files", "*.log"), ("All files", "*.*")]
    )

    if not file_path:
        print("No file selected. Exiting.")
        return

    # Prepare lists
    times = []
    uy_initial_residuals = []

    current_time = None  # Track latest time

    # Read and parse the selected log file
    with open(file_path, "r") as file:
        for line in file:
            # Match Time lines
            time_match = re.search(r'Time = ([\deE\+\-\.]+)', line)
            if time_match:
                current_time = float(time_match.group(1))

            # Match Uy residual lines
            uy_match = re.search(r'smoothSolver: *Solving for Uy, *Initial residual = ([\d\.eE\+\-]+)', line)
            if uy_match and current_time is not None:
                uy_initial_residuals.append(float(uy_match.group(1)))
                times.append(current_time)

    if not times:
        print("No Uy residual entries found in the file!")
        return

    # Plotting
    plt.figure(figsize=(10,5))
    plt.semilogy(times, uy_initial_residuals, label="Initial Residual (Uy)", marker="o")
    plt.xlabel("Time [s]")
    plt.ylabel("Initial Residual (log scale)")
    plt.title("Uy Initial Residual vs Time")
    plt.legend()
    plt.grid(True, which="both", linestyle="--", linewidth=0.5)

    # Optional Save Button (Save as PNG)
    def save_plot():
        save_path = filedialog.asksaveasfilename(defaultextension=".png",
                                                  filetypes=[("PNG files", "*.png"), ("All files", "*.*")])
        if save_path:
            plt.savefig(save_path)
            print(f"Plot saved as {save_path}")

    # Add Save Button
    save_button_ax = plt.axes([0.81, 0.01, 0.1, 0.05])
    save_button = plt.Button(save_button_ax, "Save Plot")
    save_button.on_clicked(lambda event: save_plot())

    plt.tight_layout(rect=[0, 0.05, 1, 1])  # Leave room for save button
    plt.show()

if __name__ == "__main__":
    main()

Here is a copy of the output:

(th)

tahanson43206 · 2025-04-29 09:18:10

This post is reserved for a Python program to deliver a plot of "p" residuals from an OpenFOAM log. This program ** should ** work for anyone running OpenFOAM. The value is to see the behavior of the software as it solves differential equations for every cell in the model, at each time step.

# Log viewer for "p"
import re
import tkinter as tk
from tkinter import filedialog
import matplotlib.pyplot as plt

def main():
    # Start Tkinter root (hide the main window)
    root = tk.Tk()
    root.withdraw()  # Hide the blank main window

    # Open file dialog
    file_path = filedialog.askopenfilename(
        title="Select OpenFOAM log file",
        filetypes=[("Log files", "*.log"), ("All files", "*.*")]
    )

    if not file_path:
        print("No file selected. Exiting.")
        return

    # Prepare lists
    times = []
    p_initial_residuals = []

    current_time = None  # Track latest time

    # Read and parse the selected log file
    with open(file_path, "r") as file:
        for line in file:
            # Match Time lines
            time_match = re.search(r'Time = ([\deE\+\-\.]+)', line)
            if time_match:
                current_time = float(time_match.group(1))

            # Match p residual lines
            p_match = re.search(r'Solving for p, *Initial residual = ([\d\.eE\+\-]+)', line)
            if p_match and current_time is not None:
                p_initial_residuals.append(float(p_match.group(1)))
                times.append(current_time)

    if not times:
        print("No p residual entries found in the file!")
        return

    # Plotting
    plt.figure(figsize=(10,5))
    plt.semilogy(times, p_initial_residuals, label="Initial Residual (p)", marker="o")
    plt.xlabel("Time [s]")
    plt.ylabel("Initial Residual (log scale)")
    plt.title("p Initial Residual vs Time")
    plt.legend()
    plt.grid(True, which="both", linestyle="--", linewidth=0.5)

    # Optional Save Button (Save as PNG)
    def save_plot():
        save_path = filedialog.asksaveasfilename(defaultextension=".png",
                                                  filetypes=[("PNG files", "*.png"), ("All files", "*.*")])
        if save_path:
            plt.savefig(save_path)
            print(f"Plot saved as {save_path}")

    # Add Save Button
    save_button_ax = plt.axes([0.81, 0.01, 0.1, 0.05])
    save_button = plt.Button(save_button_ax, "Save Plot")
    save_button.on_clicked(lambda event: save_plot())

    plt.tight_layout(rect=[0, 0.05, 1, 1])  # Leave room for save button
    plt.show()

if __name__ == "__main__":
    main()

Output:

(th)

tahanson43206 · 2025-04-29 09:50:34

This post is for the Python program to show changes in residuals for 'e' (energy) during the Merlin Engine test using foamRun_Vacuum model.

# 20250429PythonLogViewerFor-e-V01.py
import re
import tkinter as tk
from tkinter import filedialog
import matplotlib.pyplot as plt

def main():
    # Start Tkinter root (hide the main window)
    root = tk.Tk()
    root.withdraw()  # Hide the blank main window

    # Open file dialog
    file_path = filedialog.askopenfilename(
        title="Select OpenFOAM log file",
        filetypes=[("Log files", "*.log"), ("All files", "*.*")]
    )

    if not file_path:
        print("No file selected. Exiting.")
        return

    # Prepare lists
    times = []
    e_initial_residuals = []

    current_time = None  # Track latest time

    # Read and parse the selected log file
    with open(file_path, "r") as file:
        for line in file:
            # Match Time lines
            time_match = re.search(r'Time = ([\deE\+\-\.]+)', line)
            if time_match:
                current_time = float(time_match.group(1))

            # Match e residual lines
            e_match = re.search(r'smoothSolver: *Solving for e, *Initial residual = ([\d\.eE\+\-]+)', line)
            if e_match and current_time is not None:
                e_initial_residuals.append(float(e_match.group(1)))
                times.append(current_time)

    if not times:
        print("No e residual entries found in the file!")
        return

    # Plotting
    plt.figure(figsize=(10,5))
    plt.semilogy(times, e_initial_residuals, label="Initial Residual (e)", marker="o")
    plt.xlabel("Time [s]")
    plt.ylabel("Initial Residual (log scale)")
    plt.title("e Initial Residual vs Time")
    plt.legend()
    plt.grid(True, which="both", linestyle="--", linewidth=0.5)

    # Optional Save Button (Save as PNG)
    def save_plot():
        save_path = filedialog.asksaveasfilename(defaultextension=".png",
                                                  filetypes=[("PNG files", "*.png"), ("All files", "*.*")])
        if save_path:
            plt.savefig(save_path)
            print(f"Plot saved as {save_path}")

    # Add Save Button
    save_button_ax = plt.axes([0.81, 0.01, 0.1, 0.05])
    save_button = plt.Button(save_button_ax, "Save Plot")
    save_button.on_clicked(lambda event: save_plot())

    plt.tight_layout(rect=[0, 0.05, 1, 1])  # Leave room for save button
    plt.show()

if __name__ == "__main__":
    main()

Result goes here >>

(th)

tahanson43206 · 2025-04-30 17:25:47

This post is about OpenFOAM's recording of pressure in cells in models created to study Computational Fluid Dynamics. The model currently under study my NewMars members is of a Merlin Rocket Engine originally designed for kerosene and oxygen propellants. It is crashing shortly after hydrogen gas exits the throat of the nozzle. The current program is one of several intended to look closely at runtime parameters, in hopes of discovering possible causes for the crash.

The program below will examine a file of pressure readings and create a histogram.

# 20250430PythonReadPressureV04.py - Pressure Histogram with Bar Labels (Safe Version)
import tkinter as tk
from tkinter import filedialog
import matplotlib.pyplot as plt
import numpy as np
import re

def main():
    # Create file picker
    root = tk.Tk()
    root.withdraw()
    file_path = filedialog.askopenfilename(title="Select OpenFOAM pressure (p) file")

    if not file_path:
        print("No file selected.")
        return

    # Check filename
    filename = file_path.split('/')[-1]
    if filename != "p":
        print("⚠️  This tool is designed for OpenFOAM 'p' field files only.")
        return

    # Read pressure values
    pressures = []
    with open(file_path, 'r') as f:
        lines = f.readlines()
        if not any("internalField" in line and "nonuniform" in line for line in lines):
            print("⚠️  File structure doesn't match OpenFOAM pressure format.")
            return

        reading = False
        for line in lines:
            if reading:
                if ')' in line:
                    break
                try:
                    value = float(line.strip())
                    pressures.append(value)
                except ValueError:
                    continue
            elif re.search(r'internalField\s+nonuniform\s+List<scalar>', line):
                reading = True

    if not pressures:
        print("No pressure values found.")
        return

    # Convert and analyze
    pressure_array = np.array(pressures)
    min_p = np.min(pressure_array)
    max_p = np.max(pressure_array)
    min_bar = min_p / 100000
    max_bar = max_p / 100000
    print(f"✅ Pressure range: {min_p:.2f} Pa ({min_bar:.3f} bar) to {max_p:.2f} Pa ({max_bar:.3f} bar)")

    # Plot histogram
    plt.figure(figsize=(10, 6))
    bins = np.linspace(min_p, max_p, 100)
    plt.hist(pressure_array, bins=bins, color='blue', edgecolor='black')
    plt.title("Pressure Field Histogram")
    plt.xlabel("Pressure (Pa)")
    plt.ylabel("Number of Cells (log scale)")
    plt.yscale('log')
    plt.grid(True, which='both', linestyle='--', alpha=0.6)

    # Annotate pressure range on the plot
    textstr = f"Min: {min_p:.0f} Pa ({min_bar:.2f} bar)\nMax: {max_p:.0f} Pa ({max_bar:.2f} bar)"
    props = dict(boxstyle='round', facecolor='white', alpha=0.9)
    plt.text(0.98, 0.95, textstr, transform=plt.gca().transAxes,
             fontsize=10, verticalalignment='top', horizontalalignment='right', bbox=props)

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    main()

This space is reserved for a link to an image of output.

(th)

tahanson43206 · 2025-05-01 11:20:52

This post is to show Python code to extract sections of a polyMesh so they can be given meaningful names. The program will read the points and faces files of the polyMesh structure.

import numpy as np
import matplotlib.pyplot as plt
import csv
import re

# Define your face and point file locations
faces_file = "faces"
points_file = "points"

# Wall face range from boundary file
wall_start = 27436
wall_end = 27889

# Step 1: Read points
with open(points_file, 'r') as f:
    point_lines = f.readlines()

# Clean points and extract data
points = []
for line in point_lines:
    match = re.search(r'\(([-0-9eE+.]+) ([-0-9eE+.]+) ([-0-9eE+.]+)\)', line)
    if match:
        points.append(tuple(map(float, match.groups())))

points = np.array(points)

# Step 2: Read faces
with open(faces_file, 'r') as f:
    face_lines = f.readlines()

# Extract only faces in the wall range
wall_faces = []
count = 0
for line in face_lines:
    if wall_start <= count <= wall_end:
        match = re.findall(r'\d+', line)
        if len(match) > 1:  # Skip lone counts like "4(" lines
            indices = list(map(int, match[1:]))  # skip the face vertex count
            wall_faces.append(indices)
    count += 1

# Step 3: Calculate centroids
centroids = []
for face in wall_faces:
    coords = points[face]
    centroid = np.mean(coords, axis=0)
    centroids.append(centroid)

# Step 4: Save to CSV
with open("wall_centroids.csv", "w", newline='') as csvfile:
    writer = csv.writer(csvfile)
    writer.writerow(["FaceIndex", "X", "Y", "Z"])
    for i, c in enumerate(centroids):
        writer.writerow([i + wall_start, *c])

# Step 5: Plot
y_coords = [c[1] for c in centroids]
plt.figure(figsize=(10, 6))
plt.plot(range(wall_start, wall_end + 1), y_coords, 'bo-', markersize=2)
plt.xlabel("Face Index")
plt.ylabel("Y Centroid Position (m)")
plt.title("Y Centroids of 'wall' Faces")
plt.grid(True)
plt.tight_layout()
plt.savefig("wall_centroids_plot.png")
plt.show()

The work it will be doing is suggested by this image.

The problem is that the highlighted sections on the chart do not have names, and names are needed to resolve error messages.

(th)

tahanson43206 · 2025-05-04 12:24:24

This post is reserved for an addition to the library of utility programs to understand OpenFOAM mesh.

This version is able to find faces that include points from the points file. In this case, I used a simple spreadsheet sort to find the leftmost points in the points file. I gave that list to ChatGPT4o, and it did the massive search operation to come up with 177 faces that include at least one of the points. I am requesting a slight enhancement and hopefully that will be forthcoming. In the meantime, here is fully functional utility:

# 20250504_FindLeftFaces_V03.py
# Finds all faces that reference points from left377only.csv in the OpenFOAM mesh

def read_points_file(filepath):
    with open(filepath, "r") as f:
        lines = f.readlines()
    return [tuple(map(float, line.strip().split())) for line in lines]

def read_openfoam_points(filepath):
    with open(filepath, "r") as f:
        lines = f.readlines()
    start = lines.index("(\n") + 1
    end = lines.index(")\n")
    point_lines = lines[start:end]
    points = []
    for line in point_lines:
        coords = line.strip().strip("()").split()
        points.append(tuple(map(float, coords)))
    return points

def read_faces_file(filepath):
    with open(filepath, "r") as f:
        lines = f.readlines()
    start = lines.index("(\n") + 1
    end = lines.index(")\n")
    face_lines = lines[start:end]
    faces = []
    for line in face_lines:
        stripped = line.strip().strip("()")
        if stripped.startswith("4"):  # Match format like: 4(0 1 2 3)
            face_data = stripped[1:].strip("()").split()
            faces.append(list(map(int, face_data)))
    return faces

# Load all required data
left_points = read_points_file("left377only.csv")
all_points = read_openfoam_points("points")
faces = read_faces_file("faces")

# Build mapping from point coordinates to indices
point_index_map = {pt: i for i, pt in enumerate(all_points)}

# Identify point indices that match left_points
matched_indices = []
for pt in left_points:
    idx = point_index_map.get(pt)
    if idx is not None:
        matched_indices.append(idx)
print(f"✅ Found {len(matched_indices)} matching point indices in points file.")

# Find ALL face indices that reference any of those matched points
all_matched_faces = []
for i, face in enumerate(faces):
    for pt_idx in face:
        if pt_idx in matched_indices:
            all_matched_faces.append((i, face))
            break

print(f"✅ Found {len(all_matched_faces)} total face references using left points.\n")

# Print all results (with duplicates if any)
print("? Full face list (with possible duplicates):")
for i, face in all_matched_faces:
    print(f"  Face {i}: {face}")

Output consists of the id of each face, followed by the details for that face.

(th)

New Mars Forums

Announcement

#26 2024-07-27 17:57:27

Re: Python Computer Language

#27 2024-08-17 19:07:39

Re: Python Computer Language

#28 2024-08-18 10:56:59

Re: Python Computer Language

#29 2024-09-07 17:49:04

Re: Python Computer Language

#30 2024-09-07 18:14:43

Re: Python Computer Language

#31 2024-09-07 19:03:54

Re: Python Computer Language

#32 2024-10-07 20:03:58

Re: Python Computer Language

#33 2024-11-09 14:45:24

Re: Python Computer Language

#34 2024-11-13 23:17:46

Re: Python Computer Language

#35 2024-12-14 16:28:57

Re: Python Computer Language

#36 2025-02-25 20:46:05

Re: Python Computer Language

#37 2025-02-26 10:44:10

Re: Python Computer Language

#38 2025-04-10 07:38:21

Re: Python Computer Language

#39 2025-04-10 20:50:19

Re: Python Computer Language

#40 2025-04-23 20:19:16

Re: Python Computer Language

#41 2025-04-27 07:21:53

Re: Python Computer Language

#42 2025-04-28 13:27:01

Re: Python Computer Language

#43 2025-04-28 19:47:34

Re: Python Computer Language

#44 2025-04-29 09:18:10

Re: Python Computer Language

#45 2025-04-29 09:50:34

Re: Python Computer Language

#46 2025-04-30 17:25:47

Re: Python Computer Language

#47 2025-05-01 11:20:52

Re: Python Computer Language

#48 2025-05-04 12:24:24

Re: Python Computer Language

Board footer