convert_temperatures.py
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, Nov 20, 04:43

convert_temperatures.py
View Options

	#!/usr/bin/env python3
	# -- coding: utf-8 --
	"""
	Created on Mon Jul 25 14:38:47 2022

	@author: ekinkubilay
	"""
	import os
	import sys
	#import h5py

	sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/software/muspectre/build/language_bindings/python"))
	sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/software/muspectre/build/language_bindings/libmufft/python"))
	sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/software/muspectre/build/language_bindings/libmugrid/python"))

	import muFFT
	import muSpectre as msp

	import numpy as np
	import muGrid
	from muFFT import Stencils2D
	from muFFT import Stencils3D
	import cProfile, pstats
	import time as timer

	from mpi4py import MPI
	from dolfin import Mesh, XDMFFile, MeshTransformation, Point, HDF5File, FunctionSpace, Function

	from petsc4py import PETSc
	import argparse

	import gc
	#from transfer_previous_layer import transfer_previous_layer

	def update_temperature_field(cell, temperature_array):

	sub_domain_loc = cell.subdomain_locations

	#temperature = np.zeros((cell.nb_subdomain_grid_pts[0],cell.nb_subdomain_grid_pts[1],cell.nb_subdomain_grid_pts[2]))
	temperature = np.zeros((1,cell.nb_quad_pts, cell.nb_subdomain_grid_pts[0],cell.nb_subdomain_grid_pts[1],cell.nb_subdomain_grid_pts[2]))

	for pixel_id, pixel_coord in cell.pixels.enumerate():
	i, j, k = pixel_coord[0], pixel_coord[1], pixel_coord[2]
	a = i - sub_domain_loc[0]
	b = j - sub_domain_loc[1]
	c = k - sub_domain_loc[2]
	for q in range(cell.nb_quad_pts):
	temperature[0,q,a,b,c] = np.max([temperature_array[0,q,i,j,k],0]) + 293

	return temperature


	def parse_args():
	"""
	Function to handle arguments
	"""
	parser = argparse.ArgumentParser(description='Run given layer')
	parser.add_argument("layer_no", type=int, help='layer_number')
	parser.add_argument("mu_rank", type=int, help='number of processes for muSpectre')

	args = parser.parse_args()

	return args

	def main():

	args = parse_args()
	layer_no = args.layer_no
	os.chdir('/scratch/kubilay/mechanical_model/results')

	#comm_py is WORLD communicator (with even size), it is split into two
	#set muspectre communicator
	from mpi4py import MPI
	comm_py = MPI.COMM_WORLD
	rank_py = comm_py.rank
	size_py = comm_py.size

	#print("Hello! I'm rank %d from %d running in total..." % (comm_py.rank, comm_py.size))
	#size must be even, color is tags for the split ranks
	mu_rank = args.mu_rank
	color = np.sign(int(rank_py/(size_py-mu_rank)))


	#split communicator
	comm = comm_py.Split(color, rank_py)
	rank = comm.rank
	size = comm.size

	#print(rank_py, size_py, rank, size, color)
	#print("Total number of muSpectre tasks: %d \t Total number of fenics tasks: %d" % (comm_py.size-size, size))

	#define empty variable for broadcast buffer
	nb_steps = None
	dx, dy, dz = 0.0005, 0.0005, 0.00003
	z_lim = 0.00801
	x_lim = 0.055
	element_dimensions = [dx, dy, dz]
	#define global simulation parameters
	dim = 3
	nb_quad = 6
	#nb_domain_grid_pts = [11,11,11]

	#define constant material properties
	young = 208 #GPa
	poisson = 0.3 #unitless
	alpha_s = 0.55 #unitless
	Delta_E_bs = 1 #eV
	epsilon_sat = 0.35 # unitless
	sigma_gb = 0.150 #GPa
	sigma_s0 = 0.320 #GPa
	sigma_f = 1.240 #GPa

	## define the convergence tolerance for the Newton-Raphson increment
	newton_tol = 1e-6
	equil_tol = newton_tol
	## tolerance for the solver of the linear cell
	cg_tol = -1e-6

	## macroscopic strain
	strain_step = np.zeros((3,3)) #no applied external strain
	maxiter = 1000 # for linear cell solver

	#exact coordinates of each quadrature point relative to the local element coordinates
	quad_coords = [np.array([0.75,0.5,0.25]), np.array([0.75,0.25,0.5]), np.array([0.5,0.75,0.25]), np.array([0.25,0.75,0.5]), np.array([0.5,0.25,0.75]), np.array([0.25,0.5,0.75])]
	quad_coords = quad_coords*np.array([dx,dy,dz])

	if color == 0:
	#local size and rank on the split communicator, this communicator is called comm
	rank = comm.rank
	procs = comm.size
	other_size = size_py - size
	if rank == 0: print("Total number of muSpectre tasks: %d \t Total number of fenics tasks: %d" % (comm_py.size-size, size))
	start_time = MPI.Wtime()

	send_ranks = np.array_split(np.arange(0,other_size, 1)+size, size)[rank]
	#print(send_ranks)
	from dolfin import Mesh, XDMFFile, MeshTransformation, Point, HDF5File, FunctionSpace, Function

	#get the geometry and temperature via fenics
	mesh = Mesh(comm)

	#read the mesh and the temperature function at initial timestep
	#using Fenics within the FEM communicator

	filename = "/work/lammm/kubilay/mesh/layer_"+str(layer_no).zfill(3)+".xdmf"
	temperature_filename = "/work/lammm/kubilay/results/temperature/temperature_layer_"+str(layer_no).zfill(3)+".h5"
	f = XDMFFile(comm, filename)
	f.read(mesh)
	MeshTransformation.rescale(mesh, 1e-3, Point(0,0,0)) #in m
	f.close()
	del(f)
	gc.collect()
	bounding_box = mesh.bounding_box_tree()
	hdf5_file = HDF5File(comm, temperature_filename, "r")
	V = FunctionSpace(mesh, "CG", 1)
	temperature_field_store = Function(V)
	nb_steps = hdf5_file.attributes('Temperature')['count']-1
	timestamp_arr = np.zeros(nb_steps+1)
	for t in range(nb_steps+1):
	vec_name = "/Temperature/vector_%d"%t
	time = hdf5_file.attributes(vec_name)['timestamp']
	timestamp_arr[t] = time
	no = 0
	vec_name = "/Temperature/vector_%d"%no
	hdf5_file.read(temperature_field_store, vec_name)


	#obtain local domain dimensions
	xmax = np.max(mesh.coordinates()[:,0])
	xmin = np.min(mesh.coordinates()[:,0])
	ymax = np.max(mesh.coordinates()[:,1])
	ymin = np.min(mesh.coordinates()[:,1])
	zmax = np.max(mesh.coordinates()[:,2])
	zmin = np.min(mesh.coordinates()[:,2])
	#communicate max/min domain values within FEM communicator
	xmax = comm.allreduce(xmax, op=MPI.MAX)
	ymax = comm.allreduce(ymax, op=MPI.MAX)
	zmax = comm.allreduce(zmax, op=MPI.MAX)
	xmin = comm.allreduce(xmin, op=MPI.MIN)
	ymin = comm.allreduce(ymin, op=MPI.MIN)
	zmin = comm.allreduce(zmin, op=MPI.MIN)
	#z_lim = zmin #DELETE TO CONTROL LAYERS
	#calculate domain lengths and element dimensions
	domain_lengths = [x_lim-xmin,ymax-ymin,zmax-z_lim]
	nb_domain_grid_pts = list(np.array(domain_lengths)/np.array(element_dimensions)+1)
	nb_domain_grid_pts = [round(i) for i in nb_domain_grid_pts]

	#array to store material (and not vacuum or base plate)
	bulk = np.zeros(nb_domain_grid_pts)
	#domain discretisation TO DO:muspectre has this array already
	x_dim = np.linspace(xmin, x_lim, nb_domain_grid_pts[0])
	y_dim = np.linspace(ymin, ymax, nb_domain_grid_pts[1])
	z_dim = np.linspace(z_lim, zmax, nb_domain_grid_pts[2])

	#obtain the material (and not internal vacuum) by observing if the point
	#is in or out the submesh. This creates comm.size number of a partially
	#filled array where it is non-zero only on the sub-mesh present at that rank
	#obtain the temperature by searching every point at every sub-mesh domain
	#pass if fenics raises an error (due to point not being present on that domain)
	for i,x in enumerate(x_dim):
	for j,y in enumerate(y_dim):
	for k,z in enumerate(z_dim):
	if len(bounding_box.compute_collisions(Point(x+0.5dx,y+0.5dy,z+0.5*dz))) != 0:
	bulk[i,j,k] = 1

	del(bounding_box)
	gc.collect()

	#gather the material across all FEM ranks and broadcast to all FEM ranks
	#bulk and temperature_array are global (size corresponds to complete domain)
	bulk = comm.allreduce(bulk, op=MPI.SUM)
	bulk[bulk > 1] = 1


	#add base plate at the bottom and vacuum on top and sides
	#increase the number of grid points and total domain length in each dimension

	bulk = np.insert(bulk, 0, 2, axis=2)
	nb_domain_grid_pts[2] += 1
	domain_lengths[2] += element_dimensions[2]

	#domain discretisation TO DO:muspectre has this array already
	x_dim = np.linspace(xmin, x_lim, nb_domain_grid_pts[0])
	y_dim = np.linspace(ymin, ymax, nb_domain_grid_pts[1])
	z_dim = np.linspace(z_lim-element_dimensions[2], zmax, nb_domain_grid_pts[2])

	coords = np.zeros([nb_domain_grid_pts[0],nb_domain_grid_pts[1],nb_domain_grid_pts[2],3])
	temperature_array = np.zeros((1,nb_quad)+tuple(nb_domain_grid_pts))

	for i in range(nb_domain_grid_pts[0]):
	coords[i,:,:,0] = x_dim[i]
	for i in range(nb_domain_grid_pts[1]):
	coords[:,i,:,1] = y_dim[i]
	for i in range(nb_domain_grid_pts[2]):
	coords[:,:,i,2] = z_dim[i]

	#quadrature coordinates relative to the part
	#local bulk information on each partition
	quad_coordinates_array = np.zeros(tuple(np.shape(temperature_array))+(3,))
	#bulk_local = np.zeros_like(bulk)
	for i,x in enumerate(x_dim):
	for j,y in enumerate(y_dim):
	for k,z in enumerate(z_dim):
	if bulk[i,j,k] == 1:
	for q in range(nb_quad):
	try:
	quad_coordinates_array[0,q,i,j,k] = Point(x+quad_coords[q,0],y+quad_coords[q,1],z+quad_coords[q,2])[:]
	except RuntimeError:
	pass
	local_multi_index = []
	local_multi_index_coordinates = []
	with np.nditer([bulk], flags=['multi_index'], op_flags=['readwrite']) as it:
	for iterable in it:
	if bulk[it.multi_index] == 1:
	for q in range(nb_quad):
	final_index = tuple((0,q))+it.multi_index
	try:
	temperature_field_store(quad_coordinates_array[tuple((0,q))+it.multi_index])
	#bulk_local[it.multi_index] = bulk[it.multi_index]
	local_multi_index.append(final_index)
	local_multi_index_coordinates.append(quad_coordinates_array[final_index])
	#break
	except RuntimeError:
	pass
	# print("element_dimensions", element_dimensions)
	# print("element_dimensions", nb_domain_grid_pts)
	# print(x_dim)

	#store the temperature on the array with the right index if point is found
	#this creates comm.size number of temperature_arrays each non-zero only on
	#the sub-mesh present in that rank
	#mesh and temperature are divided differently accross communicators!! TO LOOK!
	# for i,x in enumerate(x_dim):
	# for j,y in enumerate(y_dim):
	# for k,z in enumerate(z_dim):
	# if bulk[i,j,k] == 1:
	# for q in range(nb_quad):
	# try:
	# temperature_array[0,q,i,j,k] = temperature_field_store(Point(x,y,z))
	# except RuntimeError:
	# pass

	#send global arrays (material/temperature/coordinates) from FEM communicator
	#to muSpectre communicator
	for proc_rank in send_ranks:

	comm_py.send(rank, proc_rank, tag=0)

	comm_py.send(bulk, proc_rank, tag=1)

	comm_py.send(nb_domain_grid_pts, proc_rank, tag=3)

	comm_py.send(domain_lengths, proc_rank, tag=4)

	comm_py.send(coords, proc_rank, tag=5)

	comm_py.send(timestamp_arr, proc_rank, tag=6)

	comm_py.send(temperature_array, proc_rank, tag=7)


	list_subdomain_loc = [None] * size_py
	list_nb_subdomain_grid_pts = [None] * size_py

	#print(rank_py, send_ranks)
	for f, proc_rank in enumerate(send_ranks):
	list_subdomain_loc[proc_rank] = comm_py.recv(source = proc_rank, tag=proc_rank)

	# for i,j in enumerate(list_subdomain_loc):
	# print(rank, i,j)

	for f, proc_rank in enumerate(send_ranks):
	list_nb_subdomain_grid_pts[proc_rank] = comm_py.recv(source = proc_rank, tag=proc_rank)

	# for i,j in enumerate(list_nb_subdomain_grid_pts):
	# print(rank, i,j)

	end_time = MPI.Wtime()
	#print("FENICS INIT \n Local rank = {}, took {} seconds".format(rank, end_time-start_time))

	#for muSepctre processes
	if color == 1:

	#define muSpectre communicator, this communicator is called comm_mu
	comm_mu = muGrid.Communicator(comm)
	#local size and rank on the split communicator
	rank = comm.rank
	procs = comm.size

	start_time = MPI.Wtime()

	#recieve from FEM communicator material arrays and system size parameters
	recieve_rank = comm_py.recv(source = MPI.ANY_SOURCE, tag=0)

	bulk = comm_py.recv(source = recieve_rank, tag=1)

	nb_domain_grid_pts = comm_py.recv(source = recieve_rank, tag=3)

	domain_lengths = comm_py.recv(source = recieve_rank, tag=4)

	coords = comm_py.recv(source = recieve_rank, tag=5)

	#recieve empty timestamp array from FEM communicator
	timestamp_arr = comm_py.recv(source = recieve_rank, tag=6)
	delta_t_vector = timestamp_arr[1:]-timestamp_arr[:-1]
	temperature_array = comm_py.recv(source= recieve_rank, tag=7)

	#create cell using muSpectre communicator
	cell = msp.cell.CellData.make_parallel(list(nb_domain_grid_pts), list(domain_lengths), comm_mu)
	cell.nb_quad_pts = nb_quad
	cell.get_fields().set_nb_sub_pts("quad", nb_quad)

	#define materials
	material = msp.material.MaterialHEA_3d.make(cell, "material", young, poisson, alpha_s, Delta_E_bs, epsilon_sat, sigma_gb, sigma_s0, sigma_f)
	vacuum = msp.material.MaterialLinearElastic1_3d.make(cell, "vacuum", 0.0, poisson)
	base = msp.material.MaterialLinearElastic1_3d.make(cell, "base", 20*young, poisson)

	#register fields and field collection
	material_phase = cell.get_fields().register_int_field("material2", 1, 'quad')
	bulk_array = np.zeros(np.shape(material_phase.array()), np.dtype(float))
	melt_bool = -1*np.ones(np.shape(bulk_array), np.dtype(np.int32))
	sub_domain_loc = cell.subdomain_locations
	nb_subdomain_grid_pts = cell.nb_subdomain_grid_pts
	#assign materials
	for pixel_index, pixel in enumerate(cell.pixels):
	i, j, k = pixel[0], pixel[1], pixel[2]
	a = i - sub_domain_loc[0]
	b = j - sub_domain_loc[1]
	c = k - sub_domain_loc[2]
	if bulk[tuple(pixel)]==1:
	material.add_pixel(pixel_index)
	bulk_array[0,:,a,b,c] = 1
	elif bulk[tuple(pixel)]==0:
	vacuum.add_pixel(pixel_index)
	elif bulk[tuple(pixel)]==2:
	base.add_pixel(pixel_index)

	#define solver parameters
	gradient = Stencils3D.linear_finite_elements
	weights = 6 * [1]

	material.identify_temperature_loaded_quad_pts(cell, bulk_array.reshape(-1,1, order="F"))
	#material.dt = 2.5e-5
	material.T_m = 2000
	material.alpha_thermal = 20e-6


	#local_bulk_array = bulk[sub_domain_loc[0]:sub_domain_loc[0]+nb_subdomain_grid_pts[0], sub_domain_loc[1]:sub_domain_loc[1]+nb_subdomain_grid_pts[1], sub_domain_loc[2]:sub_domain_loc[2]+nb_subdomain_grid_pts[2]]

	## use solver
	krylov_solver = msp.solvers.KrylovSolverCG(cg_tol, maxiter)
	solver = msp.solvers.SolverNewtonCG(cell, krylov_solver, msp.Verbosity.Silent, newton_tol, equil_tol, maxiter, gradient, weights)
	solver.formulation = msp.Formulation.small_strain
	solver.initialise_cell()


	#define file input/output parameters for muSpectre output
	# output_filename = "layer_"+str(layer_no-1).zfill(3)+".nc"
	# if comm_mu.rank == 0 and os.path.exists(output_filename):
	# os.remove(output_filename)
	output_filename = "converted_temperature_"+str(layer_no).zfill(3)+".nc"
	if comm_mu.rank == 0 and os.path.exists(output_filename):
	os.remove(output_filename)
	comm_mu.barrier()

	file_io_object_write = muGrid.FileIONetCDF(output_filename,
	muGrid.FileIONetCDF.OpenMode.Write,
	comm_mu)

	#coordinates = cell.get_fields().register_real_field("coordinates", 3, 'quad')
	#temperature_field = cell.get_fields().register_real_field("temperature2", 1, 'quad')
	#vacuum_strain_field = cell.get_fields().register_real_field("vacuum_strain2", (3,3), "quad")

	file_io_object_write.register_field_collection(field_collection=cell.fields)
	file_io_object_write.register_field_collection(field_collection=material.collection)
	#print("cell field names:\n", cell.get_field_names())
	#initialize global time attribute (to be updated once filled)

	file_io_object_write.write_global_attribute('timestamp', timestamp_arr)


	#create pointer-like copies for fields
	material_phase_alias = np.array(material_phase, copy=False)
	#coordinates_alias = np.array(coordinates, copy=False)
	#temperature_alias = np.array(temperature_field, copy=False)
	#vacuum_strain_alias = np.array(vacuum_strain_field, copy=False)
	#melt_bool = -1*np.ones(np.shape(material_phase), np.dtype(np.int32))
	#assign values to the above created aliases
	for pixel_id, pixel_coord in cell.pixels.enumerate():
	i, j, k = pixel_coord[0], pixel_coord[1], pixel_coord[2]
	for q in range(nb_quad):
	# local coordinates on the processor
	a = i - sub_domain_loc[0]
	b = j - sub_domain_loc[1]
	c = k - sub_domain_loc[2]

	material_phase_alias[0, q, a, b, c] = int(bulk[i, j, k])
	#coordinates_alias[0, q, a, b, c] = coords[i,j,k][0]+quad_coords[q,0]
	#coordinates_alias[1, q, a, b, c] = coords[i,j,k][1]+quad_coords[q,1]
	#coordinates_alias[2, q, a, b, c] = coords[i,j,k][2]+quad_coords[q,2]
	#if k == nb_domain_grid_pts[2]-2:
	#melt_bool[0,q,a,b,c] = 0
	#use the temperature array to change the values using temperature
	#array pointer (alias)
	#temperature_alias[:,:,:,:,:] = update_temperature_field(cell, temperature_array)[:,:,:,:,:]

	comm_mu.barrier()
	if layer_no not in [10,30,50,70,60,90,130,170,210,250,270,290,330]:

	from transfer_previous_layer import transfer_previous_layer

	strain_field = cell.get_fields().get_field("grad")
	strain_field_alias = np.array(strain_field, copy=False)
	stress_field = cell.get_fields().get_field("flux")
	stress_field_alias = np.array(stress_field, copy=False)

	#previous_cell = transfer_previous_layer(layer_no-1)
	[previous_stress, previous_strain, previous_epsilon_bar, previous_melt_pool, previous_material_temperature, previous_epsilon_plastic, previous_material_vacuum_strain] = \
	transfer_previous_layer(layer_no-1, sub_domain_loc, nb_subdomain_grid_pts)

	#assign values to the above created aliases
	for pixel_id, pixel_coord in cell.pixels.enumerate():
	i, j, k = pixel_coord[0], pixel_coord[1], pixel_coord[2]
	a = i - sub_domain_loc[0]
	b = j - sub_domain_loc[1]
	c = k - sub_domain_loc[2]
	if k < nb_domain_grid_pts[2]-1:
	for q in range(nb_quad):
	# local coordinates on the processor

	strain_field_alias[:, :, q, a, b, c] = previous_strain[:, :, q, i, j, k]
	stress_field_alias[:, :, q, a, b, c] = previous_stress[:, :, q, i, j, k]
	#temperature_alias[0,q,a,b,c] = previous_temperature[0,q,i,j,k]
	#vacuum_strain_alias[:, :, q, a, b, c] = previous_vacuum_strain[:, :, q, i, j, k]
	#if k == nb_domain_grid_pts[2]-2:
	#vacuum_strain_alias[:, :, q, a, b, c] = previous_strain[:, :, q, i, j, k]


	material.update_temperature_field(cell, previous_material_temperature.reshape(-1,1, order="F"))
	material.update_epsilon_bar_field(cell, previous_epsilon_bar.reshape(-1,1, order="F"))
	material.update_melt_bool_field(cell, previous_melt_pool.reshape(-1,1, order="F"))
	material.update_epsilon_plastic_field(cell, previous_epsilon_plastic.reshape(-1,1, order="F"))
	material.update_vacuum_strain_field(cell, previous_material_vacuum_strain.reshape(-1,1, order="F"))
	del(previous_stress, previous_strain, previous_epsilon_bar, previous_melt_pool, previous_material_temperature, previous_epsilon_plastic, previous_material_vacuum_strain)
	gc.collect()

	else:
	#use the temperature array to change the values using temperature
	#array pointer (alias)
	size_array = cell.get_fields().get_field("flux").array()
	material.update_temperature_field(cell, update_temperature_field(cell, temperature_array).reshape(-1,1, order="F"))
	material.update_epsilon_bar_field(cell, np.zeros_like(bulk_array).reshape(-1,1, order="F"))
	material.update_epsilon_plastic_field(cell, np.zeros_like(size_array).reshape(-1,1, order="F"))
	material.update_melt_bool_field(cell, melt_bool.reshape(-1,1, order="F"))



	comm_mu.barrier()

	comm_py.send([cell.subdomain_locations[0],cell.subdomain_locations[1],cell.subdomain_locations[2]], recieve_rank, tag=rank+size_py-procs)

	comm_py.send([cell.nb_subdomain_grid_pts[0],cell.nb_subdomain_grid_pts[1],cell.nb_subdomain_grid_pts[2]], recieve_rank, tag=rank+size_py-procs)

	end_time = MPI.Wtime()
	#print("MUSPECTRE INIT \n Local rank = {}, took {} seconds".format(rank, end_time-start_time))

	#make nb_Steps global by broadcasting
	nb_steps = comm_py.bcast(nb_steps, root = 0)

	start = timer.time()
	comm_py.barrier()
	#seperate operations for FEM and muSpectre communicators depending on ranks
	if color == 1:
	#temperature_array = np.zeros(nb_domain_grid_pts)

	start_time = MPI.Wtime()

	dummy_temperature = np.zeros((1,cell.nb_quad_pts, cell.nb_subdomain_grid_pts[0],cell.nb_subdomain_grid_pts[1],cell.nb_subdomain_grid_pts[2]))

	#for all steps
	for s in range(nb_steps):
	step_start = MPI.Wtime()
	if rank == 0: print('=== STEP {}/{} ==='.format(s+1, nb_steps))
	file_io_object_write.append_frame().write(["temperature"])
	#recieve temp_array from constant sender request
	#temperature_array = np.zeros(nb_domain_grid_pts)
	#print("MUSPECTRE Step {} --- Recieve request created: {}".format(s, MPI.Wtime()- step_start))
	#material.dt = delta_t_vector[s]
	req = comm_py.Irecv(dummy_temperature, source = recieve_rank, tag=s)
	#print("MUSPECTRE Step {} --- Recieve request done: {}".format(s, MPI.Wtime()- step_start))

	#print(np.max(cell.get_fields().get_field("grad").array()))

	#solve_start = MPI.Wtime()
	#solve machanics problem
	#res = solver.solve_load_increment(strain_step)
	comm_mu.barrier()
	#print("MUSPECTRE Rank {} - Step {} --- Total Solve time: {}".format(rank, s, MPI.Wtime()- solve_start) )

	#write all fields in .nc format TO DO:not all fields
	#file_io_object_write.append_frame().write(["flux","grad", "epsilon_bar", "epsilon_plastic","temperature", "vacuum_strain", "melt_pool", "material2"])
	#if (s+1)%5 == 0:
	#print('=== STEP {}/{} ==='.format(s+1, nb_steps))
	#write all fields in .nc format TO DO:not all fields
	#file_io_object_write.append_frame().write(["coordinates", "flux", "grad", "material2","temperature2", "epsilon_bar", "temperature", "vacuum_strain", "epsilon_plastic"])
	#if (s+1) == 100: break
	#print("MUSPECTRE Step {} Rank {} --- Step done: {}".format(s, rank, MPI.Wtime()- step_start))

	#make sure temperature array is recieved
	req.wait()
	#print("MUSPECTRE Step {} --- Recieve request Wait done: {}".format(s, MPI.Wtime()- step_start))

	#use the temperature array to change the values using temperature
	#array pointer (alias)
	#temperature_alias[:,:,:,:,:] = dummy_temperature
	#dummy_temperature = np.ones_like(dummy_temperature)2005 - 4s
	material.update_temperature_field(cell, dummy_temperature.reshape(-1,1, order="F"))
	#if (s+1)%100 == 0:
	#if rank == 0: print('=== STEP {}/{} ==='.format(s, nb_steps))
	#break
	#res = solver.solve_load_increment(strain_step)
	#write all fields in .nc format TO DO:not all fields
	file_io_object_write.append_frame().write(["temperature"])

	end_time = MPI.Wtime()
	if rank == 0: print("Total MUSPECTRE time: {}".format(end_time-start_time))

	comm_mu.barrier()

	if color == 0:

	start_time = MPI.Wtime()
	dummy_array = np.zeros((1,nb_quad)+tuple(nb_domain_grid_pts))
	for s in range(nb_steps):

	#start constant sending request
	#print("FENICS Step {} --- Send Requests created: {}".format(s, MPI.Wtime()-start_time))
	#reqs = [None] * len(send_ranks)
	#for p, proc_rank in enumerate(send_ranks):
	# reqs[p] = comm_py.Isend([temperature_array, MPI.DOUBLE], proc_rank, tag=s)
	#print("FENICS Step {} --- Send Requests done: {}".format(s, MPI.Wtime()-start_time))
	#re-read the temperature for the next step
	no += 1
	vec_name = "/Temperature/vector_%d"%no
	timestamp = hdf5_file.attributes(vec_name)['timestamp']
	timestamp_arr[s] = timestamp
	hdf5_file.read(temperature_field_store, vec_name)
	#store it on a dummy array since we don't want to change
	#temperature array before the sending/recieving is done
	#which will be ensured by MPI.wait
	#dummy_array = np.zeros((1,nb_quad)+tuple(nb_domain_grid_pts))

	#store the temperature on the array with the right index if point is found
	#this creates comm.size number of temperature_arrays each non-zero only on
	#the sub-mesh present in that rank
	#mesh and temperature are divided differently accross communicators!! TO LOOK!
	#with np.nditer([bulk_local, dummy_array[0,0,:,:,:]], flags=['multi_index'], op_flags=[['readonly', 'arraymask'], ['writeonly', 'writemasked']]) as it:
	# for iterable in it:
	# if bulk_local[it.multi_index] == 1:
	# for q in range(nb_quad):
	# try:
	# var = temperature_field_store(quad_coordinates_array[tuple((0,q))+it.multi_index])
	# if var < 0:
	# pass
	# else:
	# dummy_array[tuple((0,q))+it.multi_index] = var
	# except RuntimeError:
	# pass
	for lmi_index, lmi in enumerate(local_multi_index):
	val = temperature_field_store(quad_coordinates_array[lmi])
	if val <= 0:
	dummy_array[lmi] = 0
	else:
	dummy_array[lmi] = val
	#make sure temperature sending is complete and create the new
	#temperature array
	if s != 0:
	for p, proc_rank in enumerate(send_ranks):
	reqs[p].wait()
	#print("FENICS Step {} --- Send Requests Wait done: {}".format(s, MPI.Wtime()-start_time))

	temperature_array = comm.allreduce(dummy_array, op=MPI.SUM)
	temperature_array = temperature_array + 293

	#start constant sending request
	reqs = [None] * len(send_ranks)
	for p, proc_rank in enumerate(send_ranks):
	send_temp = np.array(temperature_array[:,:,list_subdomain_loc[proc_rank][0]:list_subdomain_loc[proc_rank][0]+list_nb_subdomain_grid_pts[proc_rank][0],list_subdomain_loc[proc_rank][1]:list_subdomain_loc[proc_rank][1]+list_nb_subdomain_grid_pts[proc_rank][1],list_subdomain_loc[proc_rank][2]:list_subdomain_loc[proc_rank][2]+list_nb_subdomain_grid_pts[proc_rank][2]])
	#print(p,proc_rank, list_subdomain_loc[proc_rank], list_nb_subdomain_grid_pts[proc_rank])
	#print('send_shape', np.shape(send_temp))
	reqs[p] = comm_py.Isend([send_temp, MPI.DOUBLE], proc_rank, tag=s)

	#if (s+1)%100 == 0:
	#print('=== STEP {}/{} ==='.format(s+1, nb_steps))
	#break
	#print("FENICS Step {} --- Step done: {}".format(s, MPI.Wtime()-start_time))

	comm.barrier()

	end_time = MPI.Wtime()
	if rank == 0: print("Total FENICS time: {}".format(end_time-start_time))

	#make timestamp array global by broadcasting
	timestamp_arr = comm_py.bcast(timestamp_arr, root = 0)

	#if on the muSpectre communicator, update the global timestamp
	if color == 1:
	file_io_object_write.update_global_attribute('timestamp', 'timestamp', timestamp_arr)
	file_io_object_write.close()
	comm_py.barrier()
	if rank_py == 0: print("Total run time for 100 steps:", timer.time()-start)




	if __name__=='__main__':


	main()

convert_temperatures.pyNo OneTemporaryActions

File Metadata

convert_temperatures.pyView Options

Event Timeline

convert_temperatures.py
No OneTemporary
Actions

convert_temperatures.py
View Options