FEniCS Version: 2019.2.0.dev0
Problem Description:
I am solving the Heat Equation in 1D, u_t - u_{xx} = f(x, t), with a manufactured solution, u(x, t) = \cos(t)\sin(\pi x), using a continuous Galerkin method in space-time. This means there are integrals over space and time in the weak form. I use Gaussian quadrature to handle the time integration portion with Lagrangian basis functions. I use FEniCS to handle the spatial integration component.
I solve the Heat equation and store the solution function at each time step in a list. I do this because I need this information later in my broader code to interpolate the solution at times in between time nodes. I have attached a MWE of how this is done in my broader code.
Issue:
I have run into an issue where in a broader code, I run out of MPI Communicators since you can only have 2048 in a single process. I found this thread mentioning this limit when initially investigating the error.
Goal
I want to reduce the number of MPI Communicators necessary to as low as possible while still generating the same solution list unchanged in the MWE.
MWE
from fenics import *
import numpy as np
set_log_level(LogLevel.WARNING)
def Forward_global_solver(t, u_initial, f, a_con, b_con, q, V, bc, Vh, Vh_true_deg):
t_start = t[0]
u_initial.t = t_start
alpha_init = interpolate(u_initial, Vh)
alpha_n = interpolate(u_initial, Vh)
alpha_list = [alpha_init]
for n in range(len(t) - 1):
tL = t[n]
tR = t[n + 1]
alpha_loc = Forward_local_solver(alpha_n, f, a_con, b_con, tL, tR, q, V, bc, Vh_true_deg)
for i in range(q):
if q == 1:
alpha_list.append(alpha_loc)
else:
alpha_list.append(alpha_loc.sub(i, deepcopy=True))
if q == 1:
alpha_n.assign(alpha_loc)
else:
alpha_n.assign(alpha_loc.sub(q - 1, deepcopy=True))
return alpha_list
def Forward_local_solver(alpha_n, f, a_con, b_con, tL, tR, q, V, bc, Vh):
(p, w) = np.polynomial.legendre.leggauss(3 * q)
p = ((tR - tL) / 2.0) * p + (tL + tR) / 2.0
w = ((tR - tL) / 2.0) * w
alpha_tri = TrialFunctions(V)
beta = TestFunctions(V)
F = Constant(0.0)
for j in range(q):
for k in range(len(p)):
TLagr = Constant(TestLagrange(tL, tR, q, j, p[k]))
f_k = interpolate(Expression(f, degree=6, t=p[k], a=a_con, b=b_con), Vh)
p1 = Constant(0.0)
for i in range(q + 1):
Lagr = Constant(Lagrange(tL, tR, q, i, p[k]))
dLagr = Constant(Lagrange_deriv(tL, tR, q, i, p[k]))
if i == 0:
p1 += dLagr * alpha_n * beta[j] + Lagr * alpha_n.dx(0) * beta[j].dx(0)
else:
p1 += dLagr * alpha_tri[i - 1] * beta[j] + Lagr * alpha_tri[i - 1].dx(0) * beta[j].dx(0)
weight = Constant(w[k])
F += weight * TLagr * (p1 - f_k * beta[j])
F = F * dx
a1, L = lhs(F), rhs(F)
alpha_tri = Function(V)
solve(a1 == L, alpha_tri, bc)
return alpha_tri
def Create_vector_element_function_space(msh, basis, q, deg):
if q < 1:
raise Exception('Degree of temporal basis functions must be greater than or equal to 1')
if q == 1:
V = FunctionSpace(msh, basis, deg)
else:
element = Create_vector_element_multiple_single(q, deg)
V = FunctionSpace(msh, element)
return V
def Create_zero_Dirichlet_bc(q, V):
if q < 1:
raise Exception('Degree of temporal basis functions must be greater than or equal to 1')
if q == 1:
tup = (0.0)
else:
tup = ()
tup_single_system = (0.0,)
for i in range(q):
tup = tup + tup_single_system
u_D = Constant(tup)
bc = DirichletBC(V, u_D, boundary_D)
return bc
def boundary_D(x, on_boundary):
return on_boundary
def Create_vector_element_multiple_single(q, deg):
P = FiniteElement('CG', interval, deg)
element_array = []
for i in range(q):
element_array.append(P)
vector_element = MixedElement(element_array)
return vector_element
def Lagrange(initial, final, degree, j, x):
I = np.linspace(initial, final, degree + 1)
z = 1.0
for k in range(degree + 1):
if k != j:
z *= (x - I[k]) / (I[j] - I[k])
return z
def Lagrange_deriv(initial, final, degree, j, x):
I = np.linspace(initial, final, degree + 1)
z = 0.0
for m in range(degree + 1):
v = 1.0
for k in range(degree + 1):
if k != j and k != m:
v *= (x - I[k]) / (I[j] - I[k])
if m != j:
z += (1.0 / (I[j] - I[m])) * v
return z
def TestLagrange(initial, final, degree, j, x):
I = np.linspace(initial, final, degree)
z = 1.0
if degree > 1:
for k in range(degree):
if k != j:
z *= (x - I[k]) / (I[j] - I[k])
return z
# Set up constants for time/space meshes, degrees of basis functions.
t0 = 0.0
T = 1.0
N = 5
Nt = N
Nx = N
q1 = 1
deg1 = 1
deg_true = 6
# Create time and space meshes.
t1 = np.linspace(t0, T, Nt + 1)
mesh = UnitIntervalMesh(Nx)
# Create RHS functions and initial conditions for forward problem.
a = 1.0
b = 1.0
u_init = Expression('cos(a*t)*sin(b*pi*x[0])', degree=6, t=t0, a=a, b=b)
f1 = 'sin(b*pi*x[0])*(b*b*pi*pi*cos(a*t) - a*sin(a*t))'
# Create function spaces for non-vector functions.
Vh1 = FunctionSpace(mesh, 'CG', deg1)
Vh3 = FunctionSpace(mesh, 'CG', deg_true)
# Create vector spaces for solutions.
V1 = Create_vector_element_function_space(mesh, 'CG', q1, deg1)
# Create zero Dirichlet boundary conditions.
bc1 = Create_zero_Dirichlet_bc(q1, V1)
# Forward problem.
alpha = Forward_global_solver(t1, u_init, f1, a, b, q1, V1, bc1, Vh1, Vh3)