3.1. A Basic Use Case¶
To demonstrate the functionality of the postprocessor, consider the case of the 3D heat equation with
variable diffusivity. The full demo can be found in Basic.py
.
The general heat equation reads
where u typically denotes the temperature and \(\alpha\) denotes the material diffusivity.
Boundary conditions are in our example given as
and initial condition
We also use f=0, and solve the equations at the unit cube for \(t \in (0,3]\).
3.1.1. Setting up the problem¶
We start by defining a set of parameters for our problem:
from cbcpost import *
from cbcpost.utils import cbc_print
from dolfin import *
set_log_level(WARNING)
# Create parameters for problem
params = ParamDict(
T = 3.0, # End time
dt = 0.05, # Time step
theta = 0.5, # Time stepping scheme (0.5=Crank-Nicolson)
alpha0 = 10.0, # Outer diffusivity
alpha1 = 1e-3, # Inner diffusivity
amplitude = 3.0, # Amplitude of boundary condition
)
The parameters are created using the utility class ParamDict
, which extends the built-in python
dict with dot notation to access values. We use the parameters to set up the problem using FEniCS:
# Create mesh
mesh = UnitCubeMesh(21,21,21)
# Function spaces
V = FunctionSpace(mesh, "CG", 1)
u,v = TrialFunction(V), TestFunction(V)
# Time and time-stepping
t = 0.0
timestep = 0
dt = Constant(params.dt)
# Initial condition
U = Function(V)
# Define inner domain
def inside(x):
return (0.5 < x[0] < 0.8) and (0.3 < x[1] < 0.6) and (0.2 < x[2] < 0.7)
class Alpha(Expression):
"Variable conductivity expression"
def __init__(self, alpha0, alpha1):
self.alpha0 = alpha0
self.alpha1 = alpha1
def eval(self, value, x):
if inside(x):
value[0] = self.alpha1
else:
value[0] = self.alpha0
# Conductivity
alpha = project(Alpha(params.alpha0, params.alpha1), V)
# Boundary condition
u0 = Expression("ampl*sin(x[0]*2*pi*t)", t=t, ampl=params.amplitude)
bc = DirichletBC(V, u0, "on_boundary")
# Source term
f = Constant(0)
# Bilinear form
a = (1.0/dt*inner(u,v)*dx()
+ Constant(params.theta)*alpha*inner(grad(u), grad(v))*dx)
L = (1.0/dt*inner(U,v)*dx()
+ Constant(1-params.theta)*alpha*inner(grad(U), grad(v))*dx()
+ inner(f,v)*dx())
A = assemble(a)
b = assemble(L)
bc.apply(A)
3.1.2. Setting up the PostProcessor¶
To set up the use case, we specify the case directory, and asks to clean out the case directory if there is any data remaining from a previous simulation:
pp = PostProcessor(dict(casedir="Results", clean_casedir=True))
Since we’re solving for temperature, we add a SolutionField to the postprocessor:
pp.add_field(SolutionField("Temperature", dict(save=True,
save_as=["hdf5", "xdmf"],
plot=True,
plot_args=dict(range_min=-params.amplitude,
range_max=params.amplitude),
)))
Note that we pass parameters, specifying that the field is to be saved in hdf5 and xdmf formats. These formats are default for dolfin.Function-type objects. We also ask for the Field to be plotted, with plot_args specifying the plot window. These arguments are passed directly to the dolfin.plot-command.
3.1.2.1. Time derivatives and time integrals¶
We can compute both integrals and derivatives of other Fields. Here, we add the integral of temperature from t=1.0 to t=2.0, the time-average from t=0.0 to t=5.0 as well as the derivative of the temperature field.
pp.add_fields([
TimeIntegral("Temperature", dict(save=True, start_time=1.0,
end_time=2.0)),
TimeAverage("Temperature", dict(save=True, end_time=params.T)),
TimeDerivative("Temperature", dict(save=True)),
])
Again, we ask the fields to be saved. The storage formats are determined by the datatype returned from the compute-functions.
3.1.2.2. Inspecting parts of a solution¶
We can also define fields to inspect parts of other fields. For this, we use some utilities from
cbcpost.utils
.
For this problem, the domain of a different diffusivity lies entirely within the unit cube, and thus it may
make sense to view some of the interior. We start by creating (sub)meshes of the domains we wish to inspect:
from cbcpost.utils import create_submesh, create_slice
celldomains = CellFunction("size_t", mesh)
celldomains.set_all(0)
AutoSubDomain(inside).mark(celldomains, 1)
slicemesh = create_slice(mesh, (0.7,0.5,0.5), (0.0,0.0,1.0))
submesh = create_submesh(mesh, celldomains, 1)
We then add instances of the fields PointEval
, SubFunction
and Restrict
to the
postprocessor:
pp.add_fields([
PointEval("Temperature", [[0.7,0.5, 0.5]], dict(plot=True)),
SubFunction("Temperature", slicemesh, dict(plot=True,
plot_args=dict(range_min=-params.amplitude,
range_max=params.amplitude, mode="color"))),
Restrict("Temperature", submesh, dict(plot=True, save=True)),
])
3.1.2.3. Averages and norms¶
We can also compute scalars from other fields. DomainAvg
compute the average of a specified domain
(if not specified, the whole domain). Here, we compute the average temperature inside and outside the domain
of different diffusivity, as specified by the variable cell_domains:
pp.add_fields([
DomainAvg("Temperature", cell_domains=cell_domains,
indicator=1, label="inner"),
DomainAvg("Temperature", cell_domains=cell_domains,
indicator=0, label="outer"),
])
The added parameter label does that these fields are now identified by DomainAvg_Temperature-inner and DomainAvg_Temperature-inner, respectively.
We can also compute the norm of any field:
pp.add_field(Norm("Temperature", dict(save=True)))
If no norm is specified, the L2-norm (or l2-norm) is computed.
3.1.2.4. Custom fields¶
The user may also customize fields with custom computations. In this section we demonstrate two ways to compute the difference in average temperature between the two areas of different diffusivity at any given time. First, we take an approach based solely on accessing the Temperature-field:
class TempDiff1(Field):
def __init__(self, domains, ind1, ind2, *args, **kwargs):
Field.__init__(self, *args, **kwargs)
self.domains = domains
self.dx = Measure("dx", domain=self.domains.mesh(),
subdomain_data=self.domains)
self.ind1 = ind1
self.ind2 = ind2
def before_first_compute(self, get):
self.V1 = assemble(Constant(1)*self.dx(self.ind1))
self.V2 = assemble(Constant(1)*self.dx(self.ind2))
def compute(self, get):
u = get("Temperature")
T1 = 1.0/self.V1*assemble(u*self.dx(self.ind1))
T2 = 1.0/self.V2*assemble(u*self.dx(self.ind2))
return T1-T2
In this implementation we have to specify the domains, as well as compute the respective averages directly each time. However, since we already added fields to compute the averages in both domains, there is another, much less code-demanding way to do this:
class TempDiff2(Field):
def compute(self, get):
T1 = get("DomainAvg_Temperature-inner")
T2 = get("DomainAvg_Temperature-outer")
return T1-T2
Here, we use the provided get-function to access the fields named as above, and compute the difference. We add an instance of both to the potsprocessor:
pp.add_fields([
TempDiff1(cell_domains, 1, 0, dict(plot=True)),
TempDiff2(dict(plot=True)),
])
Since both these should be the same, we can check this with ErrorNorm
:
pp.add_field(
ErrorNorm("TempDiff1", "TempDiff2", dict(plot=True), name="error"),
)
We ask for the error to be plotted. Since this is a scalar, this will be done using matplotlibs pyplot-module. We also pass the keyword argument name, which overrides the default naming (which would have been ErrorNorm_TempDiff1_TempDiff2) with error.
3.1.2.5. Combining fields¶
Finally, we can also add combination of fields, provided all dependencies have already been added to the postprocessor. For example, we can compute the space average of a time-average of our field Restrict_Temperature the following way:
pp.add_fields([
TimeAverage("Restrict_Temperature"),
DomainAvg("TimeAverage_Restrict_Temperature", params=dict(save=True)),
])
If TimeAverage(“Restrict_Temperature”) is not added first, adding the DomainAvg
-field would
fail with a DependencyException
, since the postprocessor would have no knowledge of the field
TimeAverage_Restrict_Temperature.
3.1.2.6. Saving mesh and parameters¶
We choose to store the mesh, domains and parameters associated with the problem:
pp.store_mesh(mesh, cell_domains=cell_domains)
pp.store_params(params)
These will be stored to mesh.hdf5, params.pickle and params.txt in the case directory.
3.1.3. Solving the problem¶
Solving the problem is done very simply here using simple FEniCS-commands:
solver = KrylovSolver(A, "cg", "hypre_amg")
while t <= params.T+DOLFIN_EPS:
cbc_print("Time: "+str(t))
u0.t = float(t)
assemble(L, tensor=b)
bc.apply(b)
solver.solve(U.vector(), b)
# Update the postprocessor
pp.update_all({"Temperature": lambda: U}, t, timestep)
# Update time
t += float(dt)
timestep += 1
Note the single call to the postprocessor, pp.update_all, which will then execute the logic for the postprocessor. The solution Temperature is passed in a dict as a lambda-function. This lambda-function gives the user flexibility to process the solution in any way before it is used in the postprocessor. This can for example be a scaling to physical units or joining scalar functions to a vector function.
Finally, at the end of the time-loop we finalize the postprocessor through
pp.finalize_all()
This command will finalize and return values for fields such as for example time integrals.