Page MenuHomec4science
Files F83496016

converted_temperature_xdmf.py
No OneTemporary
Actions

File Metadata

Created
Tue, Sep 17, 11:20

converted_temperature_xdmf.py
View Options

#!/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')
args = parser.parse_args()
return args
def main():
args = parse_args()
layer_no = args.layer_no
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')
#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.0005, 0.0005, 0.00003
z_lim = 0.00801
x_lim = 0.055
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)
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
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
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")
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])
#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
#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.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])
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]
## 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 = "converted_temperature_"+str(layer_no).zfill(3)+".nc"
file_io_object_read = muGrid.FileIONetCDF(output_filename,
muGrid.FileIONetCDF.OpenMode.Read,
comm)
cell.get_fields().set_nb_sub_pts("quad", nb_quad)
material_phase = cell.get_fields().register_int_field("material2", 1, 'quad')
#coordinates = cell.get_fields().register_real_field("coordinates", 3, 'quad')
#temperature = 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_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()
# 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=["coordinates", "flux",
# "grad",
# "incrF",
# "material2",
# "rhs",
# "tangent", "temperature2", "vacuum_strain2"])
#file_io_object_read.register_field_collection(field_collection=material.collection, field_names=["epsilon_bar"])
#mat_fields = material.list_fields()
#file_io_object_read.register_field_collection(field_collection=material.collection)
# sub_domain_loc = cell.subdomain_locations
# #material_phase = cell.get_fields().register_int_field("material", 1, 'quad')
# material_phase_alias = np.array(material_phase, copy=False)
# coordinates_alias = np.array(coordinates, copy=False)
# for pixel_id, pixel_coord in cell.pixels.enumerate():
# i, j, k = pixel_coord[0], pixel_coord[1], pixel_coord[2]
# pix_mat = int(bulk[i, j, k])
# 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[a, b, c] = pix_mat
# coordinates_alias[0, 0, a, b, c] = i
# coordinates_alias[1, 0, a, b, c] = j
# coordinates_alias[2, 0, a, b, c] = k
timestamp_arr = file_io_object_read.read_global_attribute('timestamp')
#print(timestamp_arr)
nb_steps = len(timestamp_arr)
#nb_steps = 100
start = timer.time()
import meshio
points = msp.linear_finite_elements.get_position_3d_helper(cell, solver=solver)
nx, ny, nz = cell.nb_domain_grid_pts
xc, yc, zc = np.mgrid[:nx, :ny, :nz]
# Global node indices
def c2i(xp, yp, zp):
return xp + (nx+1) * (yp + (ny+1) * zp)
cells = cells = np.swapaxes(
[[c2i(xc, yc, zc), c2i(xc+1, yc, zc),
c2i(xc+1, yc+1, zc), c2i(xc+1, yc+1, zc+1)],
[c2i(xc, yc, zc), c2i(xc+1, yc, zc),
c2i(xc+1, yc, zc+1), c2i(xc+1, yc+1, zc+1)],
[c2i(xc, yc, zc), c2i(xc, yc+1, zc),
c2i(xc+1, yc+1, zc), c2i(xc+1, yc+1, zc+1)],
[c2i(xc, yc, zc), c2i(xc, yc+1, zc),
c2i(xc, yc+1, zc+1), c2i(xc+1, yc+1, zc+1)],
[c2i(xc, yc, zc), c2i(xc, yc, zc+1),
c2i(xc+1, yc, zc+1), c2i(xc+1, yc+1, zc+1)],
[c2i(xc, yc, zc), c2i(xc, yc, zc+1),
c2i(xc, yc+1, zc+1), c2i(xc+1, yc+1, zc+1)]],
0, 1)
cells = cells.reshape((4, -1), order='F').T
with meshio.xdmf.TimeSeriesWriter('converted_layer_'+str(layer_no).zfill(3)+'.xdmf') as writer:
writer.write_points_cells(points, {"tetra": cells})
for i in range(nb_steps):
print('=== STEP {}/{} ==='.format(i+1, nb_steps))
file_io_object_read.read(
frame=i,
field_names=["temperature"])
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)
#print('coord array length: ', np.shape(coordinates_arr), ' temp array length: ', np.shape(temperature_arr))
c_data = {
"temperature_material": np.array([temperature_material])
}
[x_0, y_0, z_0], [x_displ, y_displ, z_displ] = get_complemented_positions_class_solver(
"0d", cell, solver, np.eye(3), periodically_complemented=True)
writer.write_data(timestamp_arr[i], cell_data=c_data)
#msp.linear_finite_elements.write_3d_class(
#"ekin_tutorial_3D_spectral_output_{}.xdmf".format(i), cell, solver, cell_data=c_data)
comm_py.barrier()
#stress = cell.get_fields().get_field('flux').array()
#strain = cell.get_fields().get_field('grad').array()
#temperature = cell.get_fields().get_field('temperature').array()
#np.savez('previous_fields.npz', stress=stress, strain= strain, temperature=temperature)
file_io_object_read.close()
if __name__=='__main__':
main()

Event Timeline

c4science · Help