diff --git a/mechanical_model/convert_array.sh b/mechanical_model/convert_array.sh index 5d046ef..8359b2d 100644 --- a/mechanical_model/convert_array.sh +++ b/mechanical_model/convert_array.sh @@ -1,79 +1,78 @@ #!/bin/bash #SBATCH --qos serial #SBATCH --nodes 1 #SBATCH --ntasks 1 #SBATCH --cpus-per-task 1 #SBATCH --mem 230G #SBATCH --time 1-01:01:01 #SBATCH --error=slurm-%j.stderr #SBATCH --output=slurm-%j.stdout #SBATCH --job-name=convert_temperature #SBATCH --mail-type=NONE #SBATCH --mail-user=recep.kubilay@epfl.ch -#SBATCH --array=268-280 - +#SBATCH --array=281-333 module purge source /home/kubilay/anaconda3/bin/activate conda activate muspectre_env2 export LD_LIBRARY_PATH=/home/kubilay/pnetcdf/lib/:$LD_LIBRARY_PATH export OPENBLAS_NUM_THREADS=1 LOG_FILE="convert_temperature_mesh_"${SLURM_ARRAY_TASK_ID}".log" echo $z >$LOG_FILE python3 adjust_temperature_mesh.py ${SLURM_ARRAY_TASK_ID} >> $LOG_FILE conda deactivate echo "*** JOB FINISHED ***" >>$LOG_FILE #sbatch ./convert_layer.sh ${SLURM_ARRAY_TASK_ID} #LOG_FILE="convert_temperature_mesh.log" #START_LAYER=268 #LAYERS=268 #echo "*** STARTING JOB ***" >$LOG_FILE #echo "START PART GEOMETRY" #module purge #source /home/kubilay/anaconda3/bin/activate #conda activate freecad_env #for z in $(seq $START_LAYER 1 $LAYERS) #do #echo $z >>$LOG_FILE #python3 meshing.py $z >> $LOG_FILE #done #echo "END PART GEOMETRY" #conda deactivate #echo "START XDMF GEOMETRY" #module purge #source /home/kubilay/anaconda3/bin/activate #conda activate muspectre_env2 #for z in $(seq $START_LAYER 1 $LAYERS) #do #echo $z >>$LOG_FILE #python3 convert_mesh.py $z >> $LOG_FILE #done #echo "END TEMPERATURE TRANSFER" #conda deactivate #echo "START XDMF GEOMETRY" #module purge #source /home/kubilay/anaconda3/bin/activate #conda activate muspectre_env2 #export LD_LIBRARY_PATH=/home/kubilay/pnetcdf/lib/:$LD_LIBRARY_PATH #export OPENBLAS_NUM_THREADS=1 ##for z in $(seq $START_LAYER 1 $LAYERS) ##do #echo $z >>$LOG_FILE #python3 adjust_temperature_mesh.py $z >> $LOG_FILE ##done #echo "END TEMPERATURE TRANSFER" #conda deactivate #echo "*** JOB FINISHED ***" >>$LOG_FILE #SOURCE_DIRECTORY="/scratch/kubilay/mechanical_model/results/coarse/" #DEST_DIRECTORY="/work/lammm/kubilay/mechanical_model/" #cp -r $SOURCE_DIRECTORY $DEST_DIRECTORY diff --git a/mechanical_model/netcdf2h5_hea.py b/mechanical_model/netcdf2h5_hea.py index 32f83ce..06b05f8 100644 --- a/mechanical_model/netcdf2h5_hea.py +++ b/mechanical_model/netcdf2h5_hea.py @@ -1,374 +1,486 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Tue Aug 2 13:47:04 2022 @author: ekinkubilay """ #from netCDF4 import Dataset """ nc_file_name = 'spectral_3D_output.nc' data = Dataset(nc_file_name, 'r') print(data) grad = data.variables['grad'][:,:,:,:,:,:] nx = data.dimensions['nx'] print(type(grad)) hf = h5py.File('data.h5', 'w') hf.create_dataset('dataset_1', data=grad) hf.close() """ 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")) sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/software/muspectre/build/language_bindings/python/muSpectre")) print("current working directory: '{}'".format(os.getcwd())) import muFFT import muSpectre as msp import numpy as np import matplotlib.pyplot as plt import muGrid from muFFT import Stencils2D from muFFT import Stencils3D import cProfile, pstats import time as timer import fenics from petsc4py import PETSc import argparse from muSpectre.gradient_integration import get_complemented_positions_class_solver 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("factor", type=float, help='factor for material parametric analysis') + parser.add_argument("target_directory", help='target directory for output/input files') args = parser.parse_args() return args +def read_data_of_step(file_io_object_read, cell, solver, read_netcdf, material): + + nb_quad = cell.nb_quad_pts + nb_domain_grid_pts = cell.nb_domain_grid_pts + dim = cell.dims + alpha_thermal = material.alpha_thermal + T_m = material.T_m + epsilon_sat = 0.35 + sigma_for = 1.200 #GPa + if read_netcdf: + file_io_object_read.read( + frame=-1, + field_names=["flux", "grad", "epsilon_bar", "epsilon_plastic", "temperature", "vacuum_strain", "melt_pool", "material2"]) + + material_arr = cell.get_fields().get_field('material2').array().reshape((1, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((1, 1, -1), order='F').T.swapaxes(1, 2) + flux_arr = cell.get_fields().get_field('flux').array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, dim, -1), order='F').T.swapaxes(1, 2) + grad_arr = cell.get_fields().get_field('grad').array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, dim, -1), order='F').T.swapaxes(1, 2) + epsilon_plastic = cell.get_globalised_current_real_field("epsilon_plastic").array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, dim, -1), order='F').T.swapaxes(1, 2) + epsilon_bar = cell.get_globalised_current_real_field("epsilon_bar").array().reshape((1, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((1, 1, -1), order='F').T.swapaxes(1, 2) + vacuum_strain = cell.get_globalised_internal_real_field("vacuum_strain").array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, dim, -1), order='F').T.swapaxes(1, 2) + temperature_material = cell.get_globalised_internal_real_field("temperature").array().reshape((1, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((1, 1, -1), order='F').T.swapaxes(1, 2) + melt_pool = cell.get_globalised_internal_int_field("melt_pool").array().reshape((1, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((1, 1, -1), order='F').T.swapaxes(1, 2) + thermal_strain = np.zeros_like(grad_arr) + thermal_strain_value = (temperature_material[:,0,0] - T_m)*alpha_thermal + for i in range(3): + thermal_strain[:,i,i] = np.array([thermal_strain_value]) + elastic_strain = grad_arr - epsilon_plastic - vacuum_strain - thermal_strain + total_actual_strain = grad_arr - vacuum_strain + sigma_forest = np.zeros_like(temperature_material) + sigma_forest = sigma_for*(1-np.exp(-epsilon_bar[:,0,0]/epsilon_sat))*(1- temperature_material[:,0,0]/T_m)**0.25 + stress_deviatoric = np.zeros_like(flux_arr) + trace_array = (flux_arr[:,0,0] + flux_arr[:,1,1] + flux_arr[:,2,2])/3 + stress_deviatoric += flux_arr + for i in range(3): + stress_deviatoric[:,i,i] -= trace_array + sigma_vonmises = np.sqrt((3/2)*np.sum(stress_deviatoric*stress_deviatoric, axis=(1,2))) + c_data = { + "material": np.array([material_arr]), + "stress": np.array([flux_arr]), + "strain": np.array([grad_arr]), + "epsilon_plastic": np.array([epsilon_plastic]), + "epsilon_bar": np.array([epsilon_bar]), + "vacuum_strain": np.array([vacuum_strain]), + "temperature_material": np.array([temperature_material]), + "melt_pool": np.array([melt_pool]), + "epsilon_elastic": np.array([elastic_strain]), + "sigma_forest": np.array([sigma_forest]), + "sigma_vonmises": np.array([sigma_vonmises]), + "stress_deviatoric": np.array([stress_deviatoric]), + "total_actual_strain": np.array([total_actual_strain]) + } + + strain_field = cell.get_fields().get_field("grad") + strain_field_alias = np.array(strain_field, copy=False) + vacuum_strain_field = cell.get_globalised_internal_real_field("vacuum_strain").array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f') + strain_field_alias -= vacuum_strain_field + #strain_field_alias[:,:,:,-1,:,:] = 0 + #strain_field_alias[:,:,:,:,-1,:] = 0 + #strain_field_alias[:,:,:,:,:,-1] = 0 + #strain_field_alias[:,:,:,:,:,:] = np.roll(strain_field_alias, (nb_domain_grid_pts[0]//2, nb_domain_grid_pts[1]//2, nb_domain_grid_pts[2]//2), (3,4,5)) + [x_0, y_0, z_0], [x_displ, y_displ, z_displ] = get_complemented_positions_class_solver("0n", cell, solver, np.eye(3), periodically_complemented=True) + #x_displ = np.roll(x_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2)) + #y_displ = np.roll(y_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2)) + #z_displ = np.roll(z_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2)) + actual_displacement = np.transpose([x_displ.ravel(order='F'), y_displ.ravel(order='F'), z_displ.ravel(order='F')]) + + file_io_object_read.read(frame=-1,field_names=["grad"]) + strain_field = cell.get_fields().get_field("grad") + strain_field_alias = np.array(strain_field, copy=False) + #strain_field_alias[:,:,:,-1,:,:] = 0 + #strain_field_alias[:,:,:,:,-1,:] = 0 + #strain_field_alias[:,:,:,:,:,-1] = 0 + #strain_field_alias[:,:,:,:,:,:] = np.roll(strain_field_alias, (nb_domain_grid_pts[0]//2, nb_domain_grid_pts[1]//2, nb_domain_grid_pts[2]//2), (3,4,5)) + [x_0, y_0, z_0], [x_displ, y_displ, z_displ] = get_complemented_positions_class_solver("0n", cell, solver, np.eye(3), periodically_complemented=True) + #x_displ = np.roll(x_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2)) + #y_displ = np.roll(y_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2)) + #z_displ = np.roll(z_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2)) + total_displacement = np.transpose([x_displ.ravel(order='F'), y_displ.ravel(order='F'), z_displ.ravel(order='F')]) + + strain_field = cell.get_fields().get_field("grad") + strain_field_alias = np.array(strain_field, copy=False) + #strain_field_alias[:,:,:,-1,:,:] = 0 + #strain_field_alias[:,:,:,:,-1,:] = 0 + #strain_field_alias[:,:,:,:,:,-1] = 0 + #strain_field_alias[:,:,:,:,:,:] = np.roll(strain_field_alias, (nb_domain_grid_pts[0]//2, nb_domain_grid_pts[1]//2, nb_domain_grid_pts[2]//2), (3,4,5)) + strain_field_alias[:,:,:,:,:] = cell.get_globalised_current_real_field("epsilon_plastic").array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f')[:,:,:,:,:] + [x_0, y_0, z_0], [x_displ, y_displ, z_displ] = get_complemented_positions_class_solver("0n", cell, solver, np.eye(3), periodically_complemented=True) + #x_displ = np.roll(x_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2)) + #y_displ = np.roll(y_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2)) + #z_displ = np.roll(z_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2)) + plastic_displacement = np.transpose([x_displ.ravel(order='F'), y_displ.ravel(order='F'), z_displ.ravel(order='F')]) + + + point_data = {'actual_displacement' : actual_displacement, 'total_displacement': total_displacement, 'plastic_displacement': plastic_displacement} + file_io_object_read.read(frame=-1,field_names=["grad"]) + + return c_data, point_data + def main(): args = parse_args() layer_no = args.layer_no - + factor = args.factor + target_directory = args.target_directory + from mpi4py import MPI comm_py = MPI.COMM_WORLD comm = muGrid.Communicator(comm_py) #rank = comm.rank #procs = comm.size - os.chdir('/scratch/kubilay/mechanical_model/results/mu_layer') - + #os.chdir('/scratch/kubilay/mechanical_model/results/mu_layer') + os.chdir(target_directory) #get the geometry via fenics mesh = fenics.Mesh() 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 = fenics.XDMFFile(filename) f.read(mesh) fenics.MeshTransformation.rescale(mesh, 1e-3, fenics.Point(0,0,0)) #in m f.close() bounding_box = mesh.bounding_box_tree() #nb_steps = 20 dim = 3 nb_quad = 6 x_max = np.max(mesh.coordinates()[:,0]) x_min = np.min(mesh.coordinates()[:,0]) y_max = np.max(mesh.coordinates()[:,1]) y_min = np.min(mesh.coordinates()[:,1]) z_max = np.max(mesh.coordinates()[:,2]) z_min = np.min(mesh.coordinates()[:,2]) xmax = np.array([np.NAN]) ymax = np.array([np.NAN]) zmax = np.array([np.NAN]) xmin = np.array([np.NAN]) ymin = np.array([np.NAN]) zmin = np.array([np.NAN]) - dx, dy, dz = 0.0002, 0.0003125, 0.00003 + dx, dy, dz = 0.00045, 0.0005, 0.00003 z_lim = 0.00801 - x_lim = 0.055 + #x_lim = 0.05496 + #x_lim_min = -0.00004 + x_lim = 0.05505 + x_lim_min = 0.00005 element_dimensions = [dx, dy, dz] #nb_domain_grid_pts = [11,11,11] - domain_lengths = [x_lim-x_min,y_max-y_min,z_max-z_lim] - nb_domain_grid_pts = list(np.array(domain_lengths)/np.array(element_dimensions)+1) + domain_lengths = [x_lim-x_lim_min,y_max-y_min,z_max-z_lim] + nb_domain_grid_pts = list(np.array(domain_lengths)/np.array(element_dimensions)) nb_domain_grid_pts = [round(i) for i in nb_domain_grid_pts] #domain_lengths = [x_max-x_min,y_max-y_min,z_max-z_min] #domain_lengths = [x_max-x_min,y_max-y_min,z_max-z_min] #element_dimensions = np.array(domain_lengths)/np.array(nb_domain_grid_pts-np.ones(dim)) #define constant material properties - young = 208e9 #GPa - poisson = 0.3 #unitless + young = 162 #GPa + poisson = 0.277 #unitless alpha_s = 0.55 #unitless - Delta_E_bs = 1 #eV + Delta_E_bs = 1.13 #eV epsilon_sat = 0.35 # unitless - sigma_gb = 150e6 #GPa - sigma_s0 = 320e6 #GPa - sigma_f = 1240e6 #GPa + sigma_gb = 0.150 #GPa + sigma_s0 = 0.254 #GPa + sigma_f = 1.200 #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 - maxiter = 1000 # for linear cell solver + maxiter = 10000 # for linear cell solver #bulk = np.load("bulk.npy") #coords = np.load("coords.npy") #nb_domain_grid_pts = np.load("nb_domain_grid_pts.npy") #domain_lengths = np.load("domain_lengths.npy") - bulk = np.zeros(nb_domain_grid_pts) - x_dim = np.linspace(x_min, x_lim, nb_domain_grid_pts[0]) - y_dim = np.linspace(y_min, y_max, nb_domain_grid_pts[1]) - z_dim = np.linspace(z_lim, z_max, nb_domain_grid_pts[2]) + bulk = np.ones(nb_domain_grid_pts) + #x_dim = np.linspace(x_min, x_lim, nb_domain_grid_pts[0]) + #y_dim = np.linspace(y_min, y_max, nb_domain_grid_pts[1]) + #z_dim = np.linspace(z_lim, z_max, nb_domain_grid_pts[2]) #assign material - 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(fenics.Point(x+0.5*dx,y+0.5*dy,z+0.5*dz))) != 0: - bulk[i,j,k] = 1 + #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(fenics.Point(x+0.5*dx,y+0.5*dy,z+0.5*dz))) != 0: + # bulk[i,j,k] = 1 #add base plate at the bottom and vacuum on top and sides # bulk = np.insert(bulk, nb_domain_grid_pts[2], 0, axis=2) # nb_domain_grid_pts[2] += 1 # domain_lengths[2] += element_dimensions[2] # bulk = np.insert(bulk, nb_domain_grid_pts[1], 0, axis=1) # nb_domain_grid_pts[1] += 1 # domain_lengths[1] += element_dimensions[1] # bulk = np.insert(bulk, nb_domain_grid_pts[0], 0, axis=0) # nb_domain_grid_pts[0] += 1 # domain_lengths[0] += element_dimensions[0] #bulk[:,:,0] *= 2 bulk = np.insert(bulk, 0, 2, axis=2) nb_domain_grid_pts[2] += 1 domain_lengths[2] += element_dimensions[2] - - x_dim = np.linspace(x_min, x_lim, nb_domain_grid_pts[0]) - y_dim = np.linspace(y_min, y_max, nb_domain_grid_pts[1]) - z_dim = np.linspace(z_lim-element_dimensions[2], z_max, nb_domain_grid_pts[2]) + + bulk = np.insert(bulk, nb_domain_grid_pts[2], 0, axis=2) + nb_domain_grid_pts[2] += 1 + domain_lengths[2] += element_dimensions[2] + + bulk = np.insert(bulk, nb_domain_grid_pts[1], 0, axis=1) + nb_domain_grid_pts[1] += 1 + domain_lengths[1] += element_dimensions[1] + + bulk = np.insert(bulk, nb_domain_grid_pts[0], 0, axis=0) + nb_domain_grid_pts[0] += 1 + domain_lengths[0] += element_dimensions[0] + + #domain discretisation TO DO:muspectre has this array already + x_dim = np.linspace(x_lim_min, x_lim, nb_domain_grid_pts[0])+0.5*dx + y_dim = np.linspace(ymin, ymax, nb_domain_grid_pts[1])+0.5*dy + z_dim = np.linspace(z_lim-element_dimensions[2], zmax, nb_domain_grid_pts[2])+0.5*dz 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)) + bulk_array = np.zeros_like(temperature_array) + 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] + + s_list = [(0.005, 0.055)] + support_mask = np.ones_like(x_dim) + for i,x in enumerate(x_dim): + loc = x + for j,s in enumerate(s_list): + if loc>s[0] and loc 1] = 1 + #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] + + bulk = np.insert(bulk, nb_domain_grid_pts[2], 0, axis=2) + nb_domain_grid_pts[2] += 1 + domain_lengths[2] += element_dimensions[2] + + bulk = np.insert(bulk, nb_domain_grid_pts[1], 0, axis=1) + nb_domain_grid_pts[1] += 1 + domain_lengths[1] += element_dimensions[1] + + bulk = np.insert(bulk, nb_domain_grid_pts[0], 0, axis=0) + nb_domain_grid_pts[0] += 1 + domain_lengths[0] += element_dimensions[0] #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]) + x_dim = np.linspace(x_lim_min, x_lim, nb_domain_grid_pts[0])+0.5*dx + y_dim = np.linspace(ymin, ymax, nb_domain_grid_pts[1])+0.5*dy + z_dim = np.linspace(z_lim-element_dimensions[2], zmax, nb_domain_grid_pts[2])+0.5*dz 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] - + + s_list = [(i,i+0.001) for i in np.linspace(0.005, 0.050, 26)] + #s_list = [(0.005, 0.05)] + support_mask = np.ones_like(x_dim) + for i,x in enumerate(x_dim): + loc = x + for j,s in enumerate(s_list): + if loc>s[0] and loc710: + #file_io_object_write.append_frame().write(["flux", "grad", "epsilon_bar", "epsilon_plastic", "temperature", "vacuum_strain", "melt_pool", "material2"]) + #if (s+1)%810 == 0: #if rank == 0: print('=== STEP {}/{} ==='.format(s, nb_steps)) #file_io_object_write.append_frame().write(["flux", "grad", "epsilon_bar", "epsilon_plastic", "temperature", "vacuum_strain", "melt_pool", "material2"]) - break + #break + #if rank == 0: print("MUSPECTRE Step {} COMPLETE --- in: {}".format(s, MPI.Wtime()- starting_step)) 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(["flux", "grad", "epsilon_bar", "epsilon_plastic", "temperature", "vacuum_strain", "melt_pool", "material2"]) end_time = MPI.Wtime() - if rank == 0: print("Total MUSPECTRE time: {}".format(end_time-start_time)) + if rank == 0: print("Total MUSPECTRE time: {}".format(end_time-start)) comm_mu.barrier() if color == 0: dummy_array = np.zeros((1,nb_quad)+tuple(nb_domain_grid_pts)) for s in range(nb_steps): start_time = MPI.Wtime() + starting_step = MPI.Wtime() #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)) - if rank == 0: print("FENICS Step {} --- File Reading: {}".format(s, MPI.Wtime()-start_time)) + #if rank == 0: print("FENICS Step {} --- File Reading: {}".format(s, MPI.Wtime()-start_time)) #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! loop_time = MPI.Wtime() - 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 + #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 #if rank == 0: print("FENICS Step {} --- Looping done: {}".format(s, MPI.Wtime()-loop_time)) - #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 start_time = MPI.Wtime() if s != 0: for p, proc_rank in enumerate(send_ranks): reqs[p].wait() - #print("FENICS Step {} --- Previous Requests Wait done: {}".format(s, MPI.Wtime()-start_time)) + #if rank == 0: print("FENICS Step {} --- Previous Requests Wait done: {}".format(s, MPI.Wtime()-start_time)) start_time = MPI.Wtime() temperature_array = comm.allreduce(dummy_array, op=MPI.SUM) temperature_array = temperature_array + 293 - #print("FENICS Step {} --- Temperature Reduce done: {}".format(s, MPI.Wtime()-start_time)) + #if rank == 0: print("FENICS Step {} --- Temperature Reduce done: {}".format(s, MPI.Wtime()-start_time)) #start constant sending request start_time = MPI.Wtime() 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)%90 == 0: + #if (s+1)%810 == 0: #print('=== STEP {}/{} ==='.format(s+1, nb_steps)) - break + #break #if rank == 0: print("FENICS Step {} --- Request sending done: {}".format(s, MPI.Wtime()-start_time)) + #if rank == 0: print("FENICS Step {} COMPLETE --- in: {}".format(s, MPI.Wtime()- starting_step)) #if rank == 0: 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)) + if rank == 0: print("Total FENICS time: {}".format(end_time-start)) #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 rank_py == 0: print("Total run time:", timer.time()-start) if __name__=='__main__': main() diff --git a/mechanical_model/run_hea.sh b/mechanical_model/run_hea.sh index 7e72762..7fc43c9 100644 --- a/mechanical_model/run_hea.sh +++ b/mechanical_model/run_hea.sh @@ -1,37 +1,40 @@ #!/bin/bash -#SBATCH --qos serial +#SBATCH --qos debug #SBATCH --nodes 1 #SBATCH --ntasks 1 -#SBATCH --cpus-per-task 36 -#SBATCH --mem 234G -#SBATCH --time 00:30:10 +#SBATCH --cpus-per-task 18 +#SBATCH --mem 254G +#SBATCH --time 00:59:10 #SBATCH --error=slurm-%j.stderr #SBATCH --output=slurm-%j.stdout #SBATCH --job-name=mech_model module purge source /home/kubilay/anaconda3/bin/activate conda activate muspectre_env2 export LD_LIBRARY_PATH=/home/kubilay/pnetcdf/lib/:$LD_LIBRARY_PATH LOG_FILE="mechanical_model_hea.log" echo "*** STARTING JOB ***" >$LOG_FILE -START_LAYER=268 -LAYERS=268 -PROCS=30 +START_LAYER=324 +LAYERS=325 +PROCS=27 export OPENBLAS_NUM_THREADS=1 +FACTOR=1.2 +TARGET_DIRECTORY="/scratch/kubilay/mechanical_model/results/factor"_$FACTOR +mkdir -p $TARGET_DIRECTORY for z in $(seq $START_LAYER 1 $LAYERS) do echo $z >>$LOG_FILE -mpirun -np 36 python3 -u residual_stress_analysis_hea.py $z $PROCS >>$LOG_FILE -python3 netcdf2h5_hea.py $z >> $LOG_FILE +#mpirun -np 36 python3 residual_stress_analysis_hea.py $z $PROCS $FACTOR $TARGET_DIRECTORY >>$LOG_FILE +python3 netcdf2h5_hea.py $z $FACTOR $TARGET_DIRECTORY >> $LOG_FILE done conda deactivate echo "*** JOB FINISHED ***" >>$LOG_FILE diff --git a/mechanical_model/run_time_optimization.sh b/mechanical_model/run_time_optimization.sh index 1be26a5..473627f 100644 --- a/mechanical_model/run_time_optimization.sh +++ b/mechanical_model/run_time_optimization.sh @@ -1,41 +1,41 @@ #!/bin/bash #SBATCH --qos serial #SBATCH --nodes 1 #SBATCH --ntasks 1 #SBATCH --cpus-per-task 36 #SBATCH --mem=252G -#SBATCH --time 01:30:00 +#SBATCH --time 00:30:00 #SBATCH --error=slurm-%j.stderr #SBATCH --output=slurm-%j.stdout #SBATCH --job-name=time_optimization module purge source /home/kubilay/anaconda3/bin/activate conda activate muspectre_env2 export LD_LIBRARY_PATH=/home/kubilay/pnetcdf/lib/:$LD_LIBRARY_PATH LOG_FILE="time_optimization.log" echo "*** STARTING JOB ***" >$LOG_FILE START_LAYER=268 -LAYERS=333 +LAYERS=300 export OPENBLAS_NUM_THREADS=1 #mpirun -np 36 python3 -u residual_stress_analysis_hea2.py 269 30 >>$LOG_FILE -python3 -u netcdf2h5_hea2.py 269 >>$LOG_FILE +python3 -u netcdf2h5_hea_disp.py 299 >>$LOG_FILE for z in $(seq $START_LAYER 5 $LAYERS) do for PROCS in 18 21 24 27 30 33 do echo $z >>$LOG_FILE #mpirun -np 36 python3 -u residual_stress_analysis_hea.py $z $PROCS >>$LOG_FILE done done conda deactivate echo "*** JOB FINISHED ***" >>$LOG_FILE diff --git a/mechanical_model/transfer_previous_layer.py b/mechanical_model/transfer_previous_layer.py index 389642c..27e2d05 100644 --- a/mechanical_model/transfer_previous_layer.py +++ b/mechanical_model/transfer_previous_layer.py @@ -1,204 +1,223 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Tue Aug 2 13:47:04 2022 @author: ekinkubilay """ import os import sys #import h5py sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/muspectre/build/language_bindings/python")) sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/muspectre/build/language_bindings/libmufft/python")) sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/muspectre/build/language_bindings/libmugrid/python")) print("current working directory: '{}'".format(os.getcwd())) #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 fenics import Mesh, XDMFFile, MeshTransformation, Point #from petsc4py import PETSc #import argparse #import gc -def transfer_previous_layer(layer_no, x_max, y_max, z_max, x_min, y_min, z_min): +def transfer_previous_layer(target_directory, layer_no, x_max, y_max, z_max, x_min, y_min, z_min): from mpi4py import MPI comm_py = MPI.COMM_SELF comm = muGrid.Communicator(comm_py) - os.chdir('/scratch/kubilay/mechanical_model/results') - + #os.chdir('/scratch/kubilay/mechanical_model/results/mu_layer') + os.chdir(target_directory) #nb_steps = 20 dim = 3 nb_quad = 6 - dx, dy, dz = 0.00050, 0.00050, 0.00003 + dx, dy, dz = 0.00045, 0.0005, 0.00003 + x_max -= 0.5*dx + y_max -= 0.5*dy + z_max -= 0.5*dz + x_min -= 0.5*dx + y_min -= 0.5*dy + z_min -= 0.5*dz z_lim = 0.00801 - x_lim = 0.055 + #x_lim = 0.05496 + #x_lim_min = -0.00004 + x_lim = 0.05505 + x_lim_min = 0.00005 z_max -= dz element_dimensions = [dx, dy, dz] #nb_domain_grid_pts = [11,11,11] - domain_lengths = [x_lim-x_min,y_max-y_min,z_max-z_lim] - nb_domain_grid_pts = list(np.array(domain_lengths)/np.array(element_dimensions)+1) + domain_lengths = [x_lim-x_lim_min,y_max-y_min,z_max-z_lim] + nb_domain_grid_pts = list(np.array(domain_lengths)/np.array(element_dimensions)) nb_domain_grid_pts = [round(i) for i in nb_domain_grid_pts] #domain_lengths = [x_max-x_min,y_max-y_min,z_max-z_min] #domain_lengths = [x_max-x_min,y_max-y_min,z_max-z_min] #element_dimensions = np.array(domain_lengths)/np.array(nb_domain_grid_pts-np.ones(dim)) #define constant material properties - young = 208 #GPa - poisson = 0.3 #unitless + young = 162 #GPa + poisson = 0.277 #unitless alpha_s = 0.55 #unitless - Delta_E_bs = 1 #eV + Delta_E_bs = 1.13 #eV epsilon_sat = 0.35 # unitless sigma_gb = 0.150 #GPa - sigma_s0 = 0.320 #GPa - sigma_f = 1.240 #GPa + sigma_s0 = 0.254 #GPa + sigma_f = 1.200 #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 maxiter = 1000 # for linear cell solver #bulk = np.load("bulk.npy") #coords = np.load("coords.npy") #nb_domain_grid_pts = np.load("nb_domain_grid_pts.npy") #domain_lengths = np.load("domain_lengths.npy") #add base plate at the bottom and vacuum on top and sides # bulk = np.insert(bulk, nb_domain_grid_pts[2], 0, axis=2) # nb_domain_grid_pts[2] += 1 # domain_lengths[2] += element_dimensions[2] # bulk = np.insert(bulk, nb_domain_grid_pts[1], 0, axis=1) # nb_domain_grid_pts[1] += 1 # domain_lengths[1] += element_dimensions[1] # bulk = np.insert(bulk, nb_domain_grid_pts[0], 0, axis=0) # nb_domain_grid_pts[0] += 1 # domain_lengths[0] += element_dimensions[0] bulk = np.ones(nb_domain_grid_pts) - bulk[:,:,-1] = 0 - bulk[:,-1,:] = 0 - bulk[-1,:,:] = 0 + bulk = np.insert(bulk, 0, 2, axis=2) nb_domain_grid_pts[2] += 1 domain_lengths[2] += element_dimensions[2] + + bulk = np.insert(bulk, nb_domain_grid_pts[2], 0, axis=2) + nb_domain_grid_pts[2] += 1 + domain_lengths[2] += element_dimensions[2] + + bulk = np.insert(bulk, nb_domain_grid_pts[1], 0, axis=1) + nb_domain_grid_pts[1] += 1 + domain_lengths[1] += element_dimensions[1] + + bulk = np.insert(bulk, nb_domain_grid_pts[0], 0, axis=0) + nb_domain_grid_pts[0] += 1 + domain_lengths[0] += element_dimensions[0] #x_dim = np.linspace(x_min, x_max, nb_domain_grid_pts[0]) #y_dim = np.linspace(y_min, y_max, nb_domain_grid_pts[1]) #z_dim = np.linspace(z_min-element_dimensions[2], z_max, nb_domain_grid_pts[2]) #coords = np.zeros([nb_domain_grid_pts[0],nb_domain_grid_pts[1],nb_domain_grid_pts[2],3]) #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] cell = msp.cell.CellData.make(list(nb_domain_grid_pts), list(domain_lengths)) 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) for pixel_index, pixel in enumerate(cell.pixels): if bulk[tuple(pixel)]==1: material.add_pixel(pixel_index) elif bulk[tuple(pixel)]==0: vacuum.add_pixel(pixel_index) elif bulk[tuple(pixel)]==2: base.add_pixel(pixel_index) gradient = Stencils3D.linear_finite_elements - weights = 6 * [1] + weights = nb_quad * [1] ## use solver krylov_solver = msp.solvers.KrylovSolverCG(cg_tol, maxiter) solver = msp.solvers.SolverNewtonCG(cell, krylov_solver, msp.Verbosity.Full, newton_tol, equil_tol, maxiter, gradient, weights) solver.formulation = msp.Formulation.small_strain solver.initialise_cell() output_filename = "layer_"+str(layer_no).zfill(3)+".nc" file_io_object_read = muGrid.FileIONetCDF(output_filename, muGrid.FileIONetCDF.OpenMode.Read, comm) #register fields and field collection #material_phase = cell.get_fields().register_int_field("material2", 1, 'quad') #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") # register global fields of the cell which you want to write # print("cell field names:\n", cell.get_field_names()) #cell_field_names = cell.get_field_names() file_io_object_read.register_field_collection(field_collection=cell.get_fields(),field_names=cell.get_field_names()) file_io_object_read.register_field_collection(field_collection=material.collection) # register global fields of the cell which you want to write #print("cell field names:\n", cell.get_field_names()) #cell_field_names = cell.get_field_names() #file_io_object_read.register_field_collection(field_collection=material.collection) file_io_object_read.read( frame=-1, field_names=["flux","grad", "epsilon_bar", "epsilon_plastic","temperature", "vacuum_strain", "melt_pool"]) material.save_history_variables() previous_stress = cell.get_fields().get_field('flux').array() previous_strain = cell.get_fields().get_field('grad').array() #previous_temperature = cell.get_fields().get_field('temperature2').array() #previous_vacuum_strain = cell.get_fields().get_field('vacuum_strain2').array() previous_epsilon_bar = np.insert(cell.get_globalised_current_real_field('epsilon_bar').array(), -1, 0, axis=-1) #previous_epsilon_bar = previous_epsilon_bar[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz] previous_melt_pool = np.insert(cell.get_globalised_internal_int_field('melt_pool').array(), -1, 0, axis=-1) previous_melt_pool[:,:,:,:,-2] = -1 #previous_melt_pool = previous_melt_pool[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz] previous_material_temperature = np.insert(cell.get_globalised_internal_real_field('temperature').array(), -1, 0, axis=-1) #previous_material_temperature = previous_material_temperature[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz] #previous_native_stress = np.insert(cell.get_globalised_internal_real_field('native_stress').array(), -1, 0, axis=-1) #previous_native_stress = previous_native_stress[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz] previous_epsilon_plastic = np.insert(cell.get_globalised_current_real_field('epsilon_plastic').array(), -1, 0, axis=-1) #previous_epsilon_plastic = previous_epsilon_plastic[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz] previous_material_vacuum_strain = np.insert(cell.get_globalised_internal_real_field('vacuum_strain').array(), -1, 0, axis=-1) previous_material_vacuum_strain[:,:,:,:,-2] = previous_strain.reshape((dim*dim, nb_quad) + tuple(nb_domain_grid_pts),order='f')[:,:,:,:,-1] #previous_material_vacuum_strain = previous_material_vacuum_strain[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz] file_io_object_read.close() return [previous_stress, previous_strain, previous_epsilon_bar, previous_melt_pool, previous_material_temperature, previous_epsilon_plastic, previous_material_vacuum_strain]