where:
Note: depending on the choice of $k$ the other terms of the equation must be dimensionaly consistent.
Here we consider 2D models in the region, $0 \leqslant x \leqslant 1 $ and $ y_{0}\leqslant y \leqslant y_{1}$
with no variation in the x direction, i.e. $ \frac{\partial T}{\partial x} = 0 $ and a constant value $ k $ across the domain.
Leading the 1D general solution:
$ T = -\frac{h}{2 k}y^{2} + c_{0}y + c_{1} $
where $c_{0}, c_{1}$ are arbitrary constants found by applying each model's boundary conditions
Three models are presented below, each with an analytic solution that the numerical results are tested against.
# analytic solution definitions
def analyticTemperature(y, h, k, c0, c1):
return -h/(2.*k)*y**2 + c0*y + c1
def exactDeriv(y, h, k, c0):
return -h/k*y + c0
import underworld as uw
import glucifer
import numpy as np
uw.matplotlib_inline()
import matplotlib.pyplot as pyplot
import matplotlib.pylab as pylab
pyplot.ion() # needed to ensure pure python jobs do now hang on show()
rank = uw.rank()
# for machines without matplotlib #
make_graphs = True
try:
import matplotlib
except ImportError:
make_graphs=False
# depth range
y0 = -.60
y1 = 1.3
# build mesh and fields
mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q1"),
elementRes = (10, 20),
minCoord = (0., y0),
maxCoord = (1., y1))
tField = mesh.add_variable( nodeDofCount=1, dataType='double')
topWall = mesh.specialSets['MaxJ_VertexSet']
bottomWall = mesh.specialSets['MinJ_VertexSet']
DirichletCondition
- topWall$ T(x,y_{1}) = T_{1} $
NeumannCondition
- bottomWall.$ q \cdot n_{b} = (\,0.0\,,\, f\,) \cdot (\,0.0\,,\,-1.0\,) = - f$
**Note** The heat flow is calculated using the heat flux vector $q$ multiplied by the outward surface normal, $n$.
The bottom surface outward normal $n_{b}$ point along the negative j-axis
When the NeumannCondition
object is associated with the SteadyStateHeat
system it defines a flux along a boundary such that:
$ q \cdot n = \phi $ on $ \Gamma_{\phi} $
where:
An example implementation
nbc
nodeIndexSet=mesh.specialSets["MinJ_VertexSet"] )
Applies phi
as a flow to the tField
over the boundary vertices in the set nodeIndexSet
.
Here phi
can be an underworld.Function
or underworld.MeshVariable
type.
Arbitrary constants are:
$c_{0} = -\frac{1}{k} (\,f - hy_{0}\,) $
$c_{1} = T_{1} + \frac{h}{2 k}y_{1}^2 - c_{0}y_{1} $
T1 = 8.0 # surface temperature
k = 6.7 # diffusivity
h = 8.0 # heat production, source term
f = 2.
# analytic solution definitions
# 1 dirichlet conditions (top) + 1 neumann (bottom)
c0 = -1./k*(f-h*y0)
c1 = T1 + h/(2.0*k)*y1**2 - c0*y1
for ii in topWall:
tField.data[ii] = T1
# flag the dirichlet conditions on the topWall only
bc = uw.conditions.DirichletCondition(tField, indexSetsPerDof=(topWall) )
# define neumann condition
nbc = uw.conditions.NeumannCondition( variable=tField,
fn_flux=-f,
indexSetsPerDof=(bottomWall))
# flag the dirichlet conditions on the topWall only
bc = uw.conditions.DirichletCondition(tField, indexSetsPerDof=(topWall) )
# define heat eq. system
ss = uw.systems.SteadyStateHeat( temperatureField = tField,
fn_diffusivity = k,
fn_heating = h,
conditions = [bc, nbc] )
solver = uw.systems.Solver(ss)
solver.solve()
# create numpy arrays for analytics, parallel check first
yvals = np.zeros(mesh.specialSets['MinI_VertexSet'].count)
ycoord = np.zeros_like(yvals)
analytic = np.zeros_like(yvals)
ids = mesh.specialSets['MinI_VertexSet']
yvals[:] = tField.evaluate(ids).T
ycoord = tField.mesh.data[ids.data,[1]]
analytic = analyticTemperature(ycoord, h, k, c0, c1)
abserr = uw.utils._nps_2norm(analytic - yvals)
mag = uw.utils._nps_2norm(analytic)
relerr = abserr / mag
# measure border flux, analytic is easy, parallel check needed for numeric result
yspot = y0
ana_flux = exactDeriv(yspot,h,k,c0)
tmp = tField.fn_gradient.evaluate_global([0.2,yspot])
if tmp is not None: num_flux = tmp[0][1]
else: num_flux = 0.
from mpi4py import MPI
comm = MPI.COMM_WORLD
# assuming order in the allgather is the same
coords = comm.allgather(ycoord)
numerical = comm.allgather(yvals)
if rank == 0:
# build matplot lib graph of result only on proc 0
# 1st build exact solution hiRes
big = np.linspace(y0,y1)
cool = analyticTemperature(big, h, k, c0, c1)
pylab.rcParams[ 'figure.figsize'] = 12, 6
pyplot.plot(coords, numerical, 'o', color = 'black', label='numerical')
pyplot.plot(big, cool, color = 'red', label="exact")
pyplot.xlabel('Y coords')
pyplot.ylabel('Temperature')
pyplot.show()
if rank == 0:
threshold = 1.0e-4
yspot = y0
print "Numerical flux at y = " ,yspot,"is", num_flux
print "Exact flux at y = " ,yspot,"is", ana_flux
print "\nAbs. error = {0:.3e}".format(abserr)
print "Rel. error = {0:.3e}".format(relerr)
if relerr > threshold:
raise RuntimeError("The numerical solution is outside the error threshold of the analytic solution." \
"The Relative error was ", relerr," the threshold is ", threshold)
Numerical flux at y = -0.6 is -0.355222812219 Exact flux at y = -0.6 is -0.298507462687 Abs. error = 3.009e-04 Rel. error = 6.742e-06
$ T(x,y_{0}) = T_{0} $
$ q \cdot n_{u} = (\,0.0\,,\, f\,) \cdot (\,0.0\,,\,1.0\,) = f$
**Note** The top surface outward normal $n_{u}$ point along the j-axis
Arbitrary constants are:
$c_{0} = -\frac{1}{k} (\,f - hy_{1}\,) $
$c_{1} = T_{0} + \frac{h}{2 k}y_{0}^2 - c_{0}y_{0} $
T0 = 8.0 # surface temperature
k = 2.2 # diffusivity
h = -7.4 # heat production, source term
f = 4.0 # temperature flow, implies negative gradient
# analytic solution definitions
# 1 dirichlet conditions (top) + 1 neumann (bottom)
c0 = -1.0*(f-h*y1)/k
c1 = T0 + h/(2.0*k)*y0**2 - c0*y0
for ii in bottomWall:
tField.data[ii] = T0
# define neumann condition
nbc = uw.conditions.NeumannCondition( fn_flux=f,
variable=tField,
indexSetsPerDof=(topWall) )
# flag the dirichlet conditions on the topWall only
bc = uw.conditions.DirichletCondition(tField, indexSetsPerDof=(bottomWall) )
# define heat eq. system
ss = uw.systems.SteadyStateHeat( temperatureField = tField,
fn_diffusivity = k,
fn_heating = h,
conditions = [bc, nbc] )
solver = uw.systems.Solver(ss)
solver.solve()
# create numpy arrays for analytics
yvals = np.zeros(mesh.specialSets['MinI_VertexSet'].count)
ycoord = np.zeros_like(yvals)
analytic = np.zeros_like(yvals)
ids = mesh.specialSets['MinI_VertexSet']
yvals[:] = tField.evaluate(ids).T
ycoord = tField.mesh.data[ids.data,[1]]
analytic = analyticTemperature(ycoord, h, k, c0, c1)
abserr = uw.utils._nps_2norm(analytic - yvals)
mag = uw.utils._nps_2norm(analytic)
relerr = abserr / mag
# measure border flux, analytic is easy, parallel check needed for numeric result
yspot = y0
ana_flux = exactDeriv(yspot,h,k,c0)
tmp = tField.fn_gradient.evaluate_global([0.2,yspot])
if tmp is not None: num_flux = tmp[0][1]
else: num_flux = 0.
from mpi4py import MPI
comm = MPI.COMM_WORLD
# assuming order in the allgather is the same
coords = comm.allgather(ycoord)
numerical = comm.allgather(yvals)
if rank == 0:
# build matplot lib graph of result only on proc 0
# 1st build exact solution hiRes
big = np.linspace(y0,y1)
cool = analyticTemperature(big, h, k, c0, c1)
pylab.rcParams[ 'figure.figsize'] = 12, 6
pyplot.plot(coords, numerical, 'o', color = 'black', label='numerical')
pyplot.plot(big, cool, color = 'red', label="exact")
pyplot.xlabel('Y coords')
pyplot.ylabel('Temperature')
pyplot.show()
if rank == 0:
threshold = 1.0e-4
yspot = y1
print "Numerical flux at y = " ,yspot,"is", num_flux
print "Exact flux at y = " ,yspot,"is", ana_flux
print "\nAbs. error = {0:.3e}".format(abserr)
print "Rel. error = {0:.3e}".format(relerr)
if relerr > threshold:
raise RuntimeError("The numerical solution is outside the error threshold of the analytic solution." \
"The Relative error was ", relerr," the threshold is ", threshold)
Numerical flux at y = 1.3 is -8.04936746828 Exact flux at y = 1.3 is -8.20909090909 Abs. error = 4.497e-04 Rel. error = 2.647e-05
2D, Steady State Heat Equation with Dirichlet BC at the top and bottom surfaces.
$T(x,y_{1}) = T_{1}$
$ T(x,y_{0}) = T_{0} $
arbitrary constants are:
$ c_{0} = \frac{1}{y_{1}-y_{0}} \left[ T_{1}-T_{0}+\frac{h} {2k}(y_{1}^2-y_{0}^2) \right] $
$c_{1} = T_{1} + \frac{h}{2k}y_{1}^2 - c_{0}y_{1}$
# Model parameters
T1 = 8.0 # top surface temperature
T0 = 4.0 # bottom surface temperature
k = 0.50 # diffusivity
h = 10 # heat production, source term
# arbitrary constant given the 2 dirichlet conditions
c0 = (T1-T0+h/(2*k)*(y1**2-y0**2)) / (y1-y0)
c1 = T1 + h/(2*k)*y1**2 - c0*y1
# set boundary conditions
for ii in topWall:
tField.data[ii] = T1
for ii in bottomWall:
tField.data[ii] = T0
# flag boundary conditions
bc = uw.conditions.DirichletCondition(tField, indexSetsPerDof=(topWall+bottomWall) )
# define heat eq. system
ss = uw.systems.SteadyStateHeat( temperatureField = tField,
fn_diffusivity = k,
fn_heating = h,
conditions = [bc] )
solver = uw.systems.Solver(ss)
solver.solve()
# create numpy arrays for analytics
yvals = np.zeros(mesh.specialSets['MinI_VertexSet'].count)
ycoord = np.zeros_like(yvals)
analytic = np.zeros_like(yvals)
ids = mesh.specialSets['MinI_VertexSet']
yvals[:] = tField.evaluate(ids).T
ycoord = tField.mesh.data[ids.data,[1]]
analytic = analyticTemperature(ycoord, h, k, c0, c1)
abserr = uw.utils._nps_2norm(analytic - yvals)
mag = uw.utils._nps_2norm(analytic)
relerr = abserr / mag
# measure border flux, analytic is easy, parallel check needed for numeric result
yspot = y0
ana_flux = exactDeriv(yspot,h,k,c0)
tmp = tField.fn_gradient.evaluate_global([0.2,yspot])
if tmp is not None: num_flux = tmp[0][1]
else: num_flux = 0.
from mpi4py import MPI
comm = MPI.COMM_WORLD
# assuming order in the allgather is the same
coords = comm.allgather(ycoord)
numerical = comm.allgather(yvals)
if rank == 0:
threshold = 1.0e-4
yspot = y0
print "Numerical flux at y = " ,yspot,"is", num_flux
print "Exact flux at y=" ,yspot,"is", ana_flux
print "\nAbs. error = {0:.3e}".format(abserr)
print "Rel. error = {0:.3e}".format(relerr)
if relerr > threshold:
raise RuntimeError("The numerical solution is outside the error threshold of the analytic solution." \
"The Relative error was ", relerr," the threshold is ", threshold)
Numerical flux at y = -0.6 is 20.1552435024 Exact flux at y= -0.6 is 21.1052631579 Abs. error = 1.988e-04 Rel. error = 3.573e-06
if rank == 0:
# build matplot lib graph of result only on proc 0
# 1st build exact solution hiRes
big = np.linspace(y0,y1)
cool = analyticTemperature(big, h, k, c0, c1)
pylab.rcParams[ 'figure.figsize'] = 12, 6
pyplot.plot(coords, numerical, 'o', color = 'black', label='numerical')
pyplot.plot(big, cool, color = 'red', label="exact")
pyplot.xlabel('Y coords')
pyplot.ylabel('Temperature')
pyplot.show()