# Impossible to write this expression in ASPECT muparser format

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
start_angle_degrees = 270  # Starting angle for the perturbation in degrees

# Convert angles from degrees to radians

# 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
# 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
end

``````

The function is very bugged and is not clean like my python code.

Thank you in advance for the help
Francesco

Better result with:

``````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=160*pi/180, start_angle=0*pi/180, pi=3.14159265359
end
``````

But not what i hoped, changing angles doesnât help.

Francesco:
The formula is of course unreadable 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:

``````  if (x > 0, \
y, \
z \
)
``````

or some such style so that it is easy to see which parentheses match, which expressions you have, etc. It will be substantially easier to debug the formula this way.

The alternative is to write your own plugin. The `Function` plugin is simply a convenient way to evaluate a formula that you provide in an input file. But there are of course cases where the formula becomes so complicated that it is easier to write it in C++ than as a single formula in the `.prm` file. In that case, creating your own plugin is the way to go.

Best
Wolfgang

Hi @bangerth , thank you for the suggestion of the "", it will surely help!
Regarding the plugin, i will use it in case as a last resource, because i didnât think that a temperature perturbation would cost so a lot of time, also iâm not so good at C++ now.

Here is the temperature expression in a readable format, i hope this can help to solve the problem:

expression:

conditions for plarge and psmall:

costants:

I hope that helps

Francesco

Hi Francesco,

A third option to add Wolfgangâs suggestion.

You can also provide initial conditions (temperature, composition) to ASPECT through an ascii data file. It simply requires writing the data out from your python code in a specific order.

An example of how to do this for 2D spherical model can be found here:
chunk_ascii_initial_temperature.prm
ascii initial temperature file

Cheers,
John

I can surely try. This means that i need to select a list of model names, the first is adiabatic, and the second is the ascii. In this way my ascii file will add or subtract the existent temperature like i want, right?

I can surely try. This means that i need to select a list of model names, the first is adiabatic, and the second is the ascii. In this way my ascii file will add or subtract the existent temperature like i want, right?

Yes, if what you generate in the python script is a perturbation from the background profile (adiabatic), you would want to use both initial temperature plugins and either add or subtract the ascii values depending on how you define the sign of the perturbation.

1 Like

Thank you, Iâll try as soon as possible and iâll let you know

I take the opportunity to ask another question. Is it possible to do the same thing for a compositional field? I wanted do a simulation of earth core segregation and replicate the results of a paper. I did a function in python to plot all the little spheres of iron as a compositional field, obviously aspect refused to do that as a function. is it possible to do the same in ASCII?

Thank you for the availability

Is it possible to do the same thing for a compositional field?

Yes, you can use the exact same approach with an ascii data file. I commonly do this to prescribed initial randomized fields that are not directly to the points on the FEM grid.

Here is an example from the test suite for using this approach in a 2D box, and as you will see it will be straightforward to apply to the 2D shell case.

Hi @jbnaliboff
I regret to inform you that the idea of using an ASCII file, as promising as it might be, unfortunately did not work out.

In my opinionâŚ

âŚthe reason is as follows: my temperature function is too complex, consequently too many points are represented in the ASCII file, including radii, angles, and temperature values. This causes ASPECT to go into a loop. I started my simulation 3 hours ago, and still nothingâŚ

Unfortunately, I canât simplify it further. According to Python, I would lose too much key information and with it the resolution. Unfortunately, only the function may help me.

IDK if this is an ASPECT problem or a problem in my file, so I attach it here, but i found similar problems using also ascii files for depth dependent simulations, specially when i wanted a viscosity profile like the Steinberger one using depth dependent.
I hope we can Solve
Francesco
temperatura_data.txt (26.2 KB)

@Francyrad How large is this file? And what happens if you run ASPECT in debug mode? Oftentimes, if a file has the wrong format, you might get a meaningful error message in debug mode, but not in release mode.
Best
W.

Dear @bangerth
The file is 27 KB, but with thousands of lines.

When i run in DEBUG, still nothing happens like in RELEASE mode, at least after 30 minutesâŚ

From my observations, it appears that ASPECT may encounter difficulties when attempting to process extensive ASCII datasets. This is intriguing because there doesnât seem to be a similar issue when handling Steinberger-related data, such as viscosity, tables, and material properties, which could potentially be due to the difference in data size or the underlying programming mechanisms.

Moreover, iâve noted that examples like the âjellyfishâ run without hitches, leading to some confusion about the challenges faced with importing complex temperature perturbations in ASCII format.

Iâm hopeful that with further investigation and collaboration, we can resolve this matter effectively.

Thank you for your attention to this issue.

Best regards,
Francesco

@Francyrad Yes, 27 kB is a rather small data set It cannot have more than 27,000 lines, and thatâs a very small number. I am fairly sure that ASPECT has successfully read files that are on the order of 100s of MB.

See if you can make your model (i.e., the .prm file) as simple as possible and upload it for us to check here!

Best
W.

@bangerth here is a minimum reproducible example. I deleted a lot of stuff, but it still doesnât work.
ascii.prm (5.4 KB)
temperatura_data.txt (26.2 KB)

i hope you can fix

Francesco

@Francyrad I must admit that I do not know why ASPECT doesnât print this error for you (perhaps try running with just one process?), but here is the error I get:

``````   Loading Ascii data initial file ./temperatura_data.txt.

---------------------------------------------------------
TimerOutput objects finalize timed values printed to the
screen by communicating over MPI in their destructors.
Since an exception is currently uncaught, this
synchronization (and subsequent output) will be skipped
---------------------------------------------------------

----------------------------------------------------
Exception 'ExcMessage("Invalid coordinate in column " + Utilities::int_to_string(column_num) + " in row " + Utilities::int_to_string(row_num) + " in file " + filename + "\nThis class expects the coordinates to be structured, meaning " "the coordinate values in each coordinate direction repeat exactly " "each time. This also means each row in the data file has to have " "the same number of columns as the first row containing data.")' on rank 0 on processing:

--------------------------------------------------------
An error occurred in line <511> of file </home/fac/g/bangerth/p/deal.II/1/projects/aspect/source/structured_data.cc> in function
void aspect::Utilities::StructuredDataLookup<dim>::load_ascii(const string&, MPI_Comm) [with int dim = 2; std::string = std::__cxx11::basic_string<char>; MPI_Comm = ompi_communicator_t*]
The violated condition was:
old_value == 0. || (std::abs(old_value-temp_data) < 1e-8*std::abs(old_value))
Invalid coordinate in column 1 in row 2 in file ./temperatura_data.txt
This class expects the coordinates to be structured, meaning the
coordinate values in each coordinate direction repeat exactly each
time. This also means each row in the data file has to have the same
number of columns as the first row containing data.
``````

The header of this file is

``````# POINTS: 3 3 3
``````

I think that you need to only specify two numbers (because youâre in two dimensions), but in any case you have more than 3x3 data points in the table below. So the header does not match the table.

Best
W.

@bangerth I finally solved a problem: now ASPECT can read my ascii file

It was due to column format, which had to be perfectly aligned, and that there were NaN values, for example T=100.2.
I added that â.â due to a misunderstanding of the example file, where the T values had a . at the end.

There is another problem now. Iâm pretty sure that my ASCII should be correct and have correct datas.

I even plotted it, it should contain my perturbations to add to a background adiabatic profile:

Did i format something wrong? I attach here the scripts and the ascii file
temperature_data_ascii.txt (1.9 MB)

``````import numpy as np
import matplotlib.pyplot as plt

# Defining the parameters and the temperature function as provided in the original script
r_inner = 3481000  # Inner radius in meters
r_outer = 6371000  # Outer radius in meters
p_large = 400.0  # Amplitude of large perturbations
p_small = 30.0  # Amplitude of small perturbations
frequency_large = 3  # Frequency of large perturbations
frequency_small = 30  # Frequency of small perturbations
start_angle_degrees = 330  # Starting angle for perturbation in degrees

# Temperature function to use degrees and handle arrays
def temperature_function_degrees_array(r, theta_degrees):

temp = np.zeros_like(r)

return temp

# Generating data for visualization
theta = np.linspace(0, 360, 360)  # Angles from 0 to 360 degrees

# Creating a grid for plotting
Temp = temperature_function_degrees_array(R, Theta)

# Plotting the original temperature distribution
plt.figure(figsize=(8, 6))
plt.colorbar(contour, label='Temperature')
plt.title('Original Temperature Distribution in Polar Coordinates')
plt.xlabel('X Coordinate (m)')
plt.ylabel('Y Coordinate (m)')
plt.axis('equal')
plt.show()

# Function to write ASCII data
with open(filename, 'w') as file:
file.write("# Test data for ascii data initial conditions.\n")
file.write("# Only next line is parsed in format: [nx] [ny] [nz] because of keyword \"POINTS:\"\n")
file.write("# POINTS: 100 360 1\n")
file.write("# Columns: r phi      temperature [K]\n")

for j in range(len(theta_degrees)):
phi = theta_degrees[j]
temp = temperatures[j, i]
file.write(f"{r:<10} {phi:<10} {temp}\n")

# Saving the ASCII data to a file
ascii_filename = 'temperature_data_ascii.txt'

# Reading the ASCII data for plotting

# Plotting the ASCII data
plt.figure(figsize=(8, 6))
plt.scatter(ascii_data[:, 1], ascii_data[:, 0], c=ascii_data[:, 2], cmap='viridis', label='ASCII Data')
plt.colorbar(label='Temperature')
plt.title('Temperature Distribution from ASCII Data')
plt.xlabel('Phi (degrees)')
plt.legend()
plt.show()

ascii_filename

``````

From the look of it, I bet you have depth and angle the wrong way around in your file.
Best
W.

@bangerth I attached the wrong file:
30-3-p30-P400.txt (969.5 KB)

The format should be exactly the same as the example:

The radius increase while the angle is constant. When all the angle is covered (0 for example), the radius begin again to the lowest point.

I really donât understand why it doesnât work