Page MenuHomec4science
Files F85247262

superposition.py
No OneTemporary
Actions

File Metadata

Created
Fri, Sep 27, 19:03

superposition.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 22 16:08:08 2021
@author: kubilay
"""
from fenics import *
import os
import numpy as np
from helpers import *
import time as timer
import multiprocessing as mp
#define simulation parameters
velocity = 0.8*np.array([1, 0, 0])
velocity_mag = np.linalg.norm(velocity)
source_loc = np.array([-1.95, -1.95, 4])
alpha = 3
time = 1
V = velocity/(2*alpha)
N = 20
num_steps = 40
dt = time/num_steps
tol = 1e-14
#read the mesh
mesh = Mesh('../Part_geometry/mesh/layer_004.xml')
boundary_mesh = BoundaryMesh(mesh, 'exterior')
#determine the limits of the domain
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 function space for v
V = FunctionSpace(mesh, 'P', 1)
#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
#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)
#initial source location
class BoundarySource(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[2], z_max, tol) and ((x[0]-1.8)**2 + x[1]**2 < 0.14))
b_source = BoundarySource()
b_source.mark(boundary_markers, 3)
#redefine d in therm of the boundary markers
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
#define initial condition
T_d = Constant(0)
#define Neumann boundary conditions
g_top = Constant(0)
g_top = Constant(0.0)
g_sides = Constant(0)
source = Constant(0)
#list of boundary conditions for convenience
boundary_conditions = {0: {'Dirichlet': T_d},
1: {'Neumann': g_top},
2: {'Neumann': g_sides},
3: {'Neumann': source}}
#project the initial condition
T_0 = project(Constant(0.0), V)
#define test and trial function and the source term
T = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
#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]:
if boundary_conditions[i]['Neumann'] != 0:
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 + dt*dot(grad(T), grad(v))*dx - T_0*v*dx - dt*f*v*dx + dt*sum(integrals_N)
a, L = lhs(F), rhs(F)
#define the function and time
T = Function(V)
t = 0
#create directory to save files
try:
os.mkdir('results')
except FileExistsError:
for file in os.scandir('results'):
os.remove(file.path)
#define vtk file to save the data at each iteration
vtkfile = File('results/analytic_temperature.pvd')
def run_analytic_analysis(t):
analytical_temperature = MyExpression(time=t, velocity=velocity, source_loc=source_loc+velocity*t, alpha=alpha)
T = interpolate(analytical_temperature, V)
print(np.max(T.vector()[:]))
T.rename("Temperature", "Analytic Temperature")
vtkfile << (T,t)
start = timer.time()
#T_final = Function(V)
#solve the weak form for each timestep
for n in range(num_steps):
#increment time and recalculate time dependent functions
t += dt
T_d.t = t
g_top.t = t
g_sides.t = t
run_analytic_analysis(t)
print('it took:', timer.time()-start)

Event Timeline

c4science · Help