Homec4science
Diffusion Homework_sp4e_Peruzzo_SáezUribe (master) Browse

R9490/homework4/exercise_7.py
09c5280ab889master

Last Change
Raw File

exercise_7.py
View Options

#!/bin/env python3
import main
import numpy as np
import argparse
from scipy import optimize
def readPositions(planet_name, directory, nb_steps ):
"""
This function reads a csv file with the state of the planet system and it create a matrix with the positions of a given planet at different times
:param planet_name: input the planet name with a NON capital first letter
:param directory: input the directory
:return: a numpy array of sizes nb_steps X 3
"""
m = nb_steps
n = 3
planet_position = np.zeros((m, n), dtype=np.float64)
# we need this if condition because the name of the files in the Trajectories folder is different from the name
# of the files that we get out of the code and are stored in the direcotry dumps
if directory == "Trajectories/":
last = 11
elif directory == "dumps/":
last = 10
else:
print("wrong directory for reference or solution files")
for i in range(0, m):
if directory == "Trajectories/":
filename = directory + "step-" + '{:04}'.format(i) + ".csv"
elif directory == "dumps/":
filename = directory + "step-" + '{:05}'.format(i) + ".csv"
else:
print("wrong directory for reference or solution files")
with open(filename, 'r') as f:
data = f.readlines()
for line in data:
words = line.split()
if words[last] == planet_name:
planet_position[i, 0:n] = np.array(words[0:n], dtype=np.float64)
return planet_position
def computeError(positions_sol, positions_ref):
"""
This function takes two matrices with the positions of a given planet at different times and it computes the integral error.
:param positions_sol: a numpy array of sizes 365 X 3 that contains the computed coordinates XYZ of the planets at different times
:param positions_ref: a numpy array of sizes 365 X 3 that contains the reference coordinates XYZ of the planets at different times
:return: type float64 containing the integral error of the simulation
"""
# todo: check if the numer of rows of positions_sol is equal to the one of positions_ref
error = np.sqrt(np.sum(np.power(np.subtract(positions_sol, positions_ref), 2)))
return error
def generateInput(scale,planet_name,input_filename,output_filename):
"""
This function generates an input file from a given input file but by scaling velocity of a given planet.
:param scale: parameter to scale the velocity
:param planet_name: string with the name of the planet
:param input_filename: name of the file containing the input for the particle code to be scaled
:param output_filename: name of the file where to write the scaled input for the particle code
:return: file containing the input for the particle code for the scaled velocity
"""
out_file = open(output_filename, 'w')
separator = ' '
with open(input_filename, 'r') as f:
data = f.readlines()
for line in data:
words = line.split()
if words[10] == planet_name:
scale_velocity = np.multiply(scale, np.array(words[3:6], dtype=np.float64))
string_velocity = np.array2string(scale_velocity, separator=' ') # There is a loss of precision here
out_file.write(separator.join(words[0:3]) + string_velocity[1:-2] + " " + separator.join(words[6:11]) + '\n')
else:
out_file.write(separator.join(words) + '\n')
out_file.close()
return out_file
def launchParticles(input,nb_steps,freq):
"""
This function launches the particle code for the specific case of planets and 1 day as time step, given an input
file, the number of steps and the frequency to dump
:param input: file containing the input for the particle code
:param nb_steps: number of steps
:param freq: frequency to dump
:return: nothing
"""
main.main(nb_steps,freq,input.name,"planet",1)
def runAndComputeError(scale, planet_name, input, nb_steps, freq):
"""
This function takes an input file containing the initial conditions of the planet system, and return the integral
error for a given scale of the velocity, a given planet, and a given number of steps
:param scale: parameter to scale the velocity
:param planet_name: string with the name of the planet
:param input: file containing the input for the particle code
:param nb_steps: number of steps
:param freq: frequency to dump
:return: type float64 containing the integral error of the simulation
"""
# creating the system data for the scaled velocity
input_scale = generateInput(scale, planet_name, input.name, "Init_scaled_onthefly.csv")
# launching the particle code for the scaled velocity
launchParticles(input_scale, nb_steps, freq)
# getting the reference positions of the planet at different times
positions_ref = readPositions(planet_name, "Trajectories/", nb_steps)
# getting the computed positions of the planet at different times
positions_sol = readPositions(planet_name, "dumps/", nb_steps)
# computing the error
error = computeError(positions_sol, positions_ref)
return error
if __name__ == "__main__":
parser = argparse.ArgumentParser(description='Compute_error')
parser.add_argument('nsteps', type=int,
help='specify the number of steps to perform')
parser.add_argument('freq', type=int,
help='specify the frequency for dumps')
parser.add_argument('planet_name', type=str,
help='specify the name of the planet')
parser.add_argument('input_filename', type=str,
help='input filename or the complete path to the file if it is not in the same folder as the script.')
args = parser.parse_args()
# getting the parameters
planet_name = args.planet_name
scale = args.scale
nb_steps = args.nsteps
freq = args.freq
# getting the name of the file containing the input for the particle code to be scaled
input_filename = args.input_filename
# generate the input
input = open(input_filename, 'r')
# compute the factor to scale the velocity of Mercury
minimum = optimize.fmin(runAndComputeError, 1, args=(planet_name, input, nb_steps, freq))
# computing the error for a given scale
error = runAndComputeError(minimum, planet_name, input, nb_steps, freq)
print("The integral error between the reference and computed positions of the planet "+planet_name+
" for a scale factor equal to ", scale, " is: ")
print(error)
c4science · Help