Dear users
I usually write my perturbations functions in python to see the result that i want and then i write them in the ASPECT muparser format.
For example, i wanted to edit this CITCOM function perturbation that i found in a paper:
subsection Function
set Coordinate system = spherical
set Variable names = r,theta
set Function constants = r_inner=3481000, r_outer=6371000, p=300.0, frequency=6, w=2
set Function expression = p*sin(pi*(r-r_inner)/(r_outer-r_inner))*sin(frequency*theta/w)
end
And then adding small and large perturbations. So i wrote this python code:
import numpy as np
import matplotlib.pyplot as plt
# Parameters
r_inner = 3481000 # Inner radius in meters
r_outer = 6371000 # Outer radius in meters
p_large = 700.0 # Amplitude of large perturbations
p_small = 100.0 # Amplitude of small perturbations
frequency_large = 3 # Frequency of large perturbations
frequency_small = 30 # Frequency of small perturbations
spread_angle_large_degrees = 90 # Spread angle of large perturbations in degrees
start_angle_degrees = 270 # Starting angle for the perturbation in degrees
# Convert angles from degrees to radians
spread_angle_large = np.radians(spread_angle_large_degrees)
start_angle = np.radians(start_angle_degrees)
# Temperature function
def temperature_function(r, theta):
# Adjust theta to shift the start angle of the perturbation
adjusted_theta = (theta + start_angle) % (2 * np.pi)
# Determine the large perturbation by using a step function approach
if (adjusted_theta < spread_angle_large) or (adjusted_theta > 2 * np.pi - spread_angle_large):
# Inside the large perturbation window
temp = p_large * np.sin(np.pi * (r - r_inner) / (r_outer - r_inner)) * np.sin(frequency_large * theta)
else:
# Outside the large perturbation window, apply small perturbations
temp = p_small * np.sin(np.pi * (r - r_inner) / (r_outer - r_inner)) * np.sin(frequency_small * theta)
return temp
# Grid for r and theta values
r_values = np.linspace(r_inner, r_outer, 100)
theta_values = np.linspace(0, 2 * np.pi, 100)
r_grid, theta_grid = np.meshgrid(r_values, theta_values)
# Calculate temperature values on the grid
temperature_grid = np.vectorize(temperature_function)(r_grid, theta_grid)
# Plotting
plt.figure(figsize=(12, 6))
X, Y = r_grid * np.cos(theta_grid), r_grid * np.sin(theta_grid)
plt.contourf(X, Y, temperature_grid, 100, cmap='viridis')
plt.colorbar()
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('Adjusted Temperature Distribution in a Spherical Shell (2D) with Starting Angle')
plt.show()
To get this:
So i get 3 large perturbations and a lot of small perturbations. The problem is that I tried the whole day to put this expression into ASPECT with poor results, this is the best i could do because of the complexity:
#subsection Function
set Coordinate system = spherical
#set Variable names = r,theta
set Function constants = r_inner=3481000, r_outer=6371000, p_large=700.0, p_small=100.0, frequency_large=3*2*pi, frequency_small=30*2*pi, spread_angle_large=180*pi/180, start_angle=270*pi/180, pi=3.14159265359
set Function expression = (0.5*(1+sign(sin(frequency_large*(theta+start_angle)/2)))*sign(sin(spread_angle_large*frequency_large/2-sin(frequency_large*(theta+start_angle)/2))))*p_large*sin(pi*(r-r_inner)/(r_outer-r_inner))*sin(frequency_large*theta) + (1-(0.5*(1+sign(sin(frequency_large*(theta+start_angle)/2)))*sign(sin(spread_angle_large*frequency_large/2-sin(frequency_large*(theta+start_angle)/2)))))*p_small*sin(pi*(r-r_inner)/(r_outer-r_inner))*sin(frequency_small*theta)
end
The function is very bugged and is not clean like my python code.
Can someone please help me to convert this function?
Thank you in advance for the help
Francesco
And thatâs exactly why it is hard to debug. My suggestion would be to (i) make use of the fact that you can break lines in input files by putting a backslash at the end, and then to write the formula in such a way that you can align parentheses, say in the following way:
