Page MenuHomec4science
Files F77499629

adjust_temperature_mesh.py
No OneTemporary
Actions

File Metadata

Created
Wed, Aug 14, 20:27

adjust_temperature_mesh.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, LagrangeInterpolator
from petsc4py import PETSc
import argparse
import gc
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
os.chdir('/scratch/kubilay/mechanical_model/results/coarse/')
#comm_py is WORLD communicator (with even size), it is split into two
#set muspectre communicator
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
size = comm.size
#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")
#xdmf_file = XDMFFile(comm, 'relevant_temperature_'+str(layer_no).zfill(3)+'.xdmf')
#xdmf_file.parameters["flush_output"] = True
#xdmf_file.parameters["functions_share_mesh"] = True
#xdmf_file.parameters["rewrite_function_mesh"] = False
hdf5_file_write = HDF5File(comm, 'relevant_temperature_'+str(layer_no).zfill(3)+'.h5', "w")
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)
temperature_field_store.set_allow_extrapolation(True)
mesh_small = Mesh(comm)
name = 'temperature_coarse_'+str(layer_no).zfill(3)
directory = "/scratch/kubilay/mechanical_model/results/coarse/"
filename = directory+name+'.xdmf'
g = XDMFFile(comm, filename)
g.read(mesh_small)
g.close()
del(g)
gc.collect()
V_small = FunctionSpace(mesh_small, "CG", 1)
temperature_field_store_small = Function(V_small)
temperature_field_store_small.rename("Temperature", "Temperature")
for i in range(nb_steps+1):
vec_name = "/Temperature/vector_%d"%i
time = hdf5_file.attributes(vec_name)['timestamp']
hdf5_file.read(temperature_field_store, vec_name)
LagrangeInterpolator.interpolate(temperature_field_store_small, temperature_field_store)
#temperature_field_store_small.interpolate(temperature_field_store)
hdf5_file_write.write(temperature_field_store_small, "Temperature", time)
#xdmf_file.write(temperature_field_store_small, time)
print('COMPLETE: ', nb_steps)
hdf5_file.close()
#xdmf_file.close()
if __name__=='__main__':
main()

Event Timeline

c4science · Help