Homec4science
Diffusion Additive Manufacturing Work (optimize_step_cost) Browse

R11910/mechanical_model/solver_trial.py
7f4f712f8fc1optimize_step_cost

Last Change
Raw File

solver_trial.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 13 11:03:05 2023
@author: ekinkubilay
"""
print('kin')
import os
import sys
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"))
print("current working directory: '{}'".format(os.getcwd()))
import muFFT
import muSpectre as msp
import numpy as np
import muGrid
from muFFT import Stencils3D
from mpi4py import MPI
import time as timer
class EigenStrain3:
"""
Class whose eigen_strain_func function is used to apply eigen strains
(gel expansion)
"""
def __init__(self,eigenstrain_shape, thermal_expansion,cell):
self.cell = cell
self.pixels = self.cell.pixels
self.nb_sub_grid = self.cell.nb_subdomain_grid_pts
self.sub_loc = self.cell.subdomain_locations
self.eigenstrain_shape = eigenstrain_shape
self.alpha = thermal_expansion
self.nb_quad_pts = self.cell.nb_quad_pts
self.step = 0
def __call__(self, strain_field):
self.alpha *= 1.1
self.eigen_strain_func(strain_field)
def eigen_strain_func(self, strain_field):
#dummy_time = MPI.Wtime()
# Looping over pixels locally and apply eigen strain based on the
# global boolean field eigen_pixels_in_structure_eigen
for i in range(self.nb_sub_grid[0]):
for j in range(self.nb_sub_grid[1]):
for k in range(self.nb_sub_grid[2]):
for q in range(self.nb_quad_pts):
strain_field[:, :, q, i, j, k] -= self.alpha*self.eigenstrain_shape
#print("EIGENSTRAIN CLASS --- Total EIGENSTRAIN time: {}".format(MPI.Wtime()- dummy_time) )
comm = MPI.COMM_WORLD
comm = muGrid.Communicator(comm)
nb_steps = 5
dx, dy, dz = 0.05, 0.05, 0.05
#z_lim = 0.00216
element_dimensions = [dx, dy, dz]
domain_lengths = [1,1,1]
nb_domain_grid_pts = list(np.array(domain_lengths)/np.array(element_dimensions)+1)
nb_domain_grid_pts = [int(i) for i in nb_domain_grid_pts]
#define global simulation parameters
dim = 3
nb_quad = 8
thermal_expansion = 1e-5
#nb_domain_grid_pts = [11,11,11]
#define constant material properties
Youngs_modulus = 100e9
Poisson_ratio = 1/3
## define the convergence tolerance for the Newton-Raphson increment
newton_tol = 1e-4
equil_tol = newton_tol
## tolerance for the solver of the linear cell
cg_tol = -1e-3
## 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])
#create cell using muSpectre communicator
cell = msp.cell.CellData.make(list(nb_domain_grid_pts), list(domain_lengths))
cell.nb_quad_pts = nb_quad
#define materials
material = msp.material.MaterialLinearElastic1_3d.make(cell, "material", Youngs_modulus , Poisson_ratio)
base = msp.material.MaterialLinearElastic1_3d.make(cell, "base", 2*Youngs_modulus, Poisson_ratio)
#assign materials
for pixel_index, pixel in enumerate(cell.pixels):
if np.random.rand()>0.75:
base.add_pixel(pixel_index)
else:
material.add_pixel(pixel_index)
#define solver parameters
eigenstrain_shape = np.identity(3)
gradient = Stencils3D.linear_finite_elements
weights = 8 * [1]
## use solver
krylov_solver = msp.solvers.KrylovSolverPCG(cg_tol, maxiter)
fem_stencil = msp.FEMStencil.trilinear_hexahedron(cell)
discretisation = msp.Discretisation(fem_stencil)
solver = msp.solvers.SolverFEMNewtonPCG(discretisation, krylov_solver, msp.Verbosity.Full, newton_tol, equil_tol, maxiter)
solver.formulation = msp.Formulation.small_strain
solver.initialise_cell()
solver.evaluate_stress_tangent()
reference_material = solver.tangent.map.mean()
solver.set_reference_material(reference_material)
eigen_class = EigenStrain3(eigenstrain_shape, thermal_expansion, cell)
# output_filename = "solver_trial.nc"
# comm = muGrid.Communicator(MPI.COMM_SELF)
# file_io_object_write = muGrid.FileIONetCDF(output_filename, muGrid.FileIONetCDF.OpenMode.Write, comm)
# #register fields and field collection
# cell.get_fields().set_nb_sub_pts("quad", nb_quad)
# material_phase = cell.get_fields().register_int_field("material", 1, 'quad')
# coordinates = cell.get_fields().register_real_field("coordinates", 3, 'quad')
# temperature_field = cell.get_fields().register_real_field("temperature", 1, 'quad')
# vacuum_strain_field = cell.get_fields().register_real_field("vacuum_strain", (3,3), "quad")
# file_io_object_write.register_field_collection(field_collection=cell.get_fields(),field_names=cell.get_field_names())
start = timer.time()
for i in range(nb_steps):
res = solver.solve_load_increment(strain_step, eigen_class.eigen_strain_func)
#krylov_solver.initialise()
#write all fields in .nc format TO DO:not all fields
# file_io_object_write.append_frame().write(["flux","grad"])
print("Total time:", timer.time()-start)
c4science · Help