Homec4science
Diffusion Additive Manufacturing Work (temperature_mesh_coarsening) Browse

R11910/Thermal_model/fem_layer_scan_final.py
3e9a1d30b01ftemperature_mesh_coarse…

Last Change
Raw File

fem_layer_scan_final.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Mar 14 12:30:17 2022
@author: ekinkubilay
"""
from fenics import *
import os
import numpy as np
import time as timeit
import matplotlib.pyplot as plt
import argparse
from fem_scan_strategy import *
def parse_args():
"""
Function to handle arguments
"""
parser = argparse.ArgumentParser(description='Run FEM of single scan on a simple box')
parser.add_argument("layer_no", type=int, help='layer number')
parser.add_argument("rho", type=float, help='density in kg/m^3')
parser.add_argument("c_p", type=float, help='specific heat in J/kgK')
parser.add_argument("k", type=float, help='thermal conductivity in W/Km')
parser.add_argument("power", type=float, help='laser power in W')
parser.add_argument("vel_mag", type=float, help='velocity magnitude in m/s')
parser.add_argument("time_step", type=float, help='time step for simulation (s)')
args = parser.parse_args()
return args
def main():
#get the arguments
args = parse_args()
#define simulation parameters
layer_no = args.layer_no
dt = args.time_step
store = fem_scan_strategy(dt, layer_no)
time = store[-1,-1]
num_steps = len(store)
print('Simulation will be {} steps'.format(num_steps))
#num_steps = 100
#define simulation parameters
rho = args.rho #kg/m^3
c_p = args.c_p #J/kgK
k = args.k #W/mK
alpha = k/(rho*c_p) #m^2/s
power = args.power #W
tol = 1e-14
#read the mesh
mesh = Mesh()
mesh_dir = '../Part_geometry/mesh/'
temp_data_dir = 'results/temperature/'
mesh_name = mesh_dir+'layer_{}.xdmf'.format(str(layer_no).zfill(3))
f = XDMFFile(MPI.comm_world, mesh_name)
f.read(mesh)
MeshTransformation.rescale(mesh, 1e-3, Point(0,0,0)) #in m
f.close()
#determine the limits of the domain in m
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])
#define markers for the boundary faces
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 4)
#define boundary subclasses for bottom, top and side surfaces, mark the boundaries
class BoundaryBottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], z_min, tol)
b_bottom = BoundaryBottom()
b_bottom.mark(boundary_markers, 0)
class BoundaryTop(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], z_max, tol)
b_top = BoundaryTop()
b_top.mark(boundary_markers, 1)
class BoundarySides(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[0], x_min, tol) or near(x[0], x_max, tol) or near(x[1], y_min, tol) or near(x[1], y_max, tol))
b_sides = BoundarySides()
b_sides.mark(boundary_markers, 2)
#define function space for v
V = FunctionSpace(mesh, "CG", 1)
#redefine d in therm of the boundary markers
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
#define Dirichlet initial condition
T_d = Constant(0)
#define Neumann boundary conditions
g_top = Constant(0.0)
g_sides = Constant(0)
#list of boundary conditions for convenience
boundary_conditions = {0: {'Dirichlet': T_d},
1: {'Neumann': g_top},
2: {'Neumann': g_sides}}
#interpolate the initial condition
#T_0 = interpolate(Constant(0), V)
T_0 = Function(V)
try:
print(temp_data_dir+'adjusted_temperature_layer_{}.h5'.format(str(layer_no).zfill(3)))
hdf5_file = HDF5File(MPI.comm_world, temp_data_dir+'adjusted_temperature_layer_{}.h5'.format(str(layer_no).zfill(3)), "r")
#define the last time step
#t = hdf5_file.attributes('Temperature')['count']-1
#define name of the vector
#vec_name = "/Temperature/vector_%d"%t
#read the temperature of the last time step of previous layer
hdf5_file.read(T_0, 'Temperature')
#timestamp = hdf5_file.attributes(vec_name)["timestamp"]
hdf5_file.close()
except RuntimeError:
T_0 = interpolate(Constant(0), V)
#define test and trial function and the source term
T = TrialFunction(V)
v = TestFunction(V)
#gather all Dirichlet boundary conditions in one list
bcs = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
bcs.append(bc)
#gather all Neumann boundary conditions in weak form
integrals_N = []
for i in boundary_conditions:
if 'Neumann' in boundary_conditions[i]:
g = boundary_conditions[i]['Neumann']
integrals_N.append(g*v*ds(i))
#write the weak form and seperate into right- and left-hand sides
F = T*v*dx + alpha*dt*dot(grad(T), grad(v))*dx - T_0*v*dx - alpha*dt*sum(integrals_N)
a, L = lhs(F), rhs(F)
#define the function and time
T = Function(V)
t = 0
#create directory to save .xdmf temperature data files
visuals_dir = 'results/visualization'
temp_data_dir = 'results/temperature'
filename = '/temperature_layer_{}'.format(str(layer_no).zfill(3))
try:
os.mkdir(visuals_dir)
except FileExistsError:
pass
#create directory to .h5 temperature data files
try:
os.mkdir(temp_data_dir)
except FileExistsError:
pass
xdmf_file = XDMFFile(MPI.comm_world, visuals_dir+filename+'.xdmf')
xdmf_file.parameters["flush_output"] = True
xdmf_file.parameters["functions_share_mesh"] = True
xdmf_file.parameters["rewrite_function_mesh"] = False
hdf5_file = HDF5File(MPI.comm_world, temp_data_dir+filename+'.h5', "w")
#calculate time
start = timeit.time()
A, b = assemble_system(a, L, bcs)
#solve the weak form for each timestep
for i in range(num_steps):
#assemble the equations
b = assemble(L, tensor=b)
#determine the location of the source in given time
#apply the source if it is on
#b = M*T.vector()
try:
delta = PointSource(V, Point(store[i][0], store[i][1], z_max-tol), dt*power/(rho*c_p))
delta.apply(b)
except RuntimeError:
pass
#solve the system
solve(A, T.vector(), b, 'cg', 'sor')
#assign the temperature as the initial condition for the next step
T_0.assign(T)
#write temperature data
#if i%10==0:
T.rename("Temperature", "Temperature")
#vtkfile_fem << (T,t)
#xdmf_file.write_checkpoint(T, "FEM Temperature", t, XDMFFile.Encoding.HDF5)
xdmf_file.write(T, t)
hdf5_file.write(T, "Temperature", t)
#increment time and recalculate time dependent functions
t += dt
dt = 0.05
cooling_time = 3*time
cooling_num_steps = int(cooling_time/0.05)
#solve the weak form for each timestep
for i in range(cooling_num_steps):
b = assemble(L, tensor=b)
#solve the system
solve(A, T.vector(), b, 'cg', 'sor')
#assign the temperature as the initial condition for the next step
T_0.assign(T)
#store the temperature at points of interests
"""
T_start[i] = T(start_point)
T_middle[i] = T(mid_point)
T_end[i] = T(end_point)
"""
#write temperature data
T.rename("Temperature", "Temperature")
if i%20==0:
xdmf_file.write(T, t)
hdf5_file.write(T, "Temperature", t)
#increment time and recalculate time dependent functions
t += dt
#output total simulation time
print('it took:', timeit.time()-start)
hdf5_file.close()
xdmf_file.close()
if __name__=='__main__':
main()
c4science · Help