Hi there,
I'm looking at Gryphon project, but find that none of the given examples are mixed function space problems. I'm wondering if Gryphon works with mixed function space. Actually, I did try to solve a dynamic version of the mixed poisson problem with Gryphon, but haven't succeeded yet. It will be great if you could give me some idea where the problem is.
Provided are two minimal test codes: one with DOLFIN, working but fixed time step; the other with Gryphon, not working.
DOLFIN version:
'''
Status: working
This is to solve the dynamic Heat conduction problem
with uniform Neumann boundary conditions:
sigma - grad(u) = 0
dudt - div(sigma) = f
where
f = 10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)
Initial condition:
u = 0
sigma = (0, 0)
Boundary condition:
n*sigma = 0 on all boundary
Time interval:
T = [0, 1]
'''
from dolfin import *
t = 0
dt = 1.0e-2
T = 1.0
# Create mesh
mesh = UnitSquareMesh(32, 32)
# Define function spaces and mixed (product) space
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
W = BDM * DG
# Define trial and test functions and initial solution w0
(sigma, u) = TrialFunctions(W)
w0 = Function(W)
(sigma0, u0) = split(w0)
(tau, v) = TestFunctions(W)
# Define source function
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
# Define variational form
a = (dot(sigma, tau) + div(tau)*u + u*v/dt - div(sigma)*v)*dx
L = u0*v/dt*dx + f*v*dx
# Define function G such that G \cdot n = g
class BoundarySource(Expression):
def __init__(self, mesh):
self.mesh = mesh
def eval_cell(self, values, x, ufc_cell):
cell = Cell(self.mesh, ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
g = 0.0
values[0] = g*n[0]
values[1] = g*n[1]
def value_shape(self):
return (2,)
G = BoundarySource(mesh)
# Define essential boundary
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(W.sub(0), G, boundary)
# Compute solution
w = Function(W)
while t < T:
t += dt
solve(a == L, w, bc)
w0.vector()[:] = w.vector()
# Plot sigma and u
(sigma, u) = w.split()
plot(sigma)
plot(u)
interactive()
Gryphon version:
'''
Status: not working
This is to solve the dynamic Heat conduction problem
with uniform Neumann boundary conditions:
sigma - grad(u) = 0
dudt - div(sigma) = f
where
f = 10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)
Initial condition:
u = 0
sigma = (0, 0)
Boundary condition:
n*sigma = 0 on all boundary
Time interval:
T = [0, 1]
'''
from gryphon import ESDIRK
from dolfin import *
# initial value classes
class InitialConditions(Expression):
def eval(self, values, x):
values[0] = 0.0
values[1] = 0.0
values[2] = 0.0
def value_shape(self):
return (3,)
# Create mesh
mesh = UnitSquareMesh(32, 32)
# Define function spaces and mixed (product) space
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
W = BDM * DG
# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)
# Define source function
source = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
f = div(sigma)*v*dx + source*v*dx
g = -dot(sigma, tau)*dx - div(tau)*u*dx
# Define function G such that G \cdot n = g
class BoundarySource(Expression):
def __init__(self, mesh):
self.mesh = mesh
def eval_cell(self, values, x, ufc_cell):
cell = Cell(self.mesh, ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
g = 0.0
values[0] = g*n[0]
values[1] = g*n[1]
def value_shape(self):
return (2,)
G = BoundarySource(mesh)
# Define essential boundary
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(W.sub(0), G, boundary)
w = Function(W)
w.interpolate(InitialConditions())
T = [0,1.0]
# Create the ESDIRK object.
obj = ESDIRK(T,w,f,g=[g],bcs=bc)
obj.parameters['output']['reuseoutputfolder'] = True
# obj.parameters['timestepping']['dt'] = 1.0e-6
# Turn on extra terminal output
obj.parameters["verbose"] = True
# Turn on runtime plot of current time step
obj.parameters["drawplot"] = True
# Save runtime statistics.
obj.parameters["output"]["statistics"] = False
# Save plot of each time step in VTK format.
obj.parameters["output"]["plot"] = True
# Set that the plot of selected step sizes should be saved in jpg.
# Available choices are jpg, png and eps.
obj.parameters["output"]["imgformat"] = "jpg"
# Call the solver which will do the actual calculation.
obj.solve()
