Automate system of PDEs without using MixedElement

1

Hello, I am trying to use FEniCS to solve a system of N equations without using MixedElement, as the equations are not necessary to be solved all at once. I am attempting to use Python list to achieve the automation:

from dolfin import *
class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(MPI.comm_world))
        super().__init__(**kwargs)
    def eval(self, values, x):
        values = random.uniform(-0.001,0.001)

parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

N = 99 # some large number of equations

mesh = UnitSquareMesh(300, 300)

du = [0]*N
v  = [0]*N
u  = [0]*N
u_n= [0]*N
u_init = InitialConditions()
V  = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, V)

for i in range(N):
    # Define trial and test functions
    du[i]    = TrialFunction(ME)
    v[i]  = TestFunctions(ME)
    # Define functions
    u[i]   = Function(ME)  # current solution
    u_n[i]  = Function(ME)  # solution from previous converged step
    # Create intial conditions and interpolate
    u[i].interpolate(u_init)
    u_n[i].interpolate(u_init)

L = [0]*N
ap = [0] *N  ## ap == L , used ap instead of a
problem = [0] *N
for i in range(N):
    dsum = 0
    for a in range(N):
        for b in range(a+1,N):
            dsum += u_n[a]*u_n[b]**2

    dFdu  =  dsum *v[i] + dot(grad(u_n[i]), grad(v[i]))
    L[i] = ( u[i]- u_n[i] )*v[i]*dx  +  dt*( dsum *v[i] + dot(grad(u_n[i]), grad(v[i])) )*dx
    # Compute directional derivative about u in the direction of du (Jacobian)
    ap[i] = derivative(L[i], u[i], du[i])
    problem[i] = NSol(ap[i], L[i])  # Nsol is a class to interface with newton solver

I haven’t gotten down to the time stepping part yet, I am currently stuck at grad(v[i]) where it give me an error saying that cannot be converted to any UFL type.

Beside, there is an warning from my class of initial condition, mentioning that: WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.

Does anyone have some inside on how to fix these issues and automate this process? Thanks!

I am also including my tentative codes for the time stepping part with Newton solver:

# Class for interfacing with the Newton solver
class NSol(NonlinearProblem): ## this will be moved to the beginning part of the code
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)


# neglected some codes setting up PETSc 

t = 0.0
dt = 1e-2
T = 10*dt
while (t < T):
    t += dt
    for i in range(N):
        u_n[i].vector()[:] = u[i].vector()
    for i in range(N):
        solver.solve(problem[i], u[i].vector())
        # neglected some code exporting the solutions
2

For the first problem, there is an insidious typo, where you wrote TestFunctions instead of TestFunction, and the former returns a tuple, which cannot be converted to a UFL type.

For the second problem, you can add the method

    def value_shape(self):
        return ()

to get rid of the warning. You probably also want to change eval to something like

    def eval(self, values, x):
        values[0] = random.uniform(-0.001,0.001)

where an element of the NumPy array values is written to.