I am trying to implement the following corrector equation for a fixed \omega.
where I define A(y,\omega) as \kappa below and regard x as y_1, regard y as y_2.
With the inspiration of this link, I wrote the following code.
from dolfin import *
import ufl
import matplotlib.pyplot as plt
import numpy
# Create mesh and define function space
mesh = UnitSquareMesh(4, 6)
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
R = FiniteElement('Real', mesh.ufl_cell(), 0)
X = FunctionSpace(mesh, P1*R)
# Define subdomains
subdomains = MeshFunction("size_t", mesh, 2)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return True if x[1] <= 0.5 else False
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return True if x[1] >= 0.5 else False
# Mark subdomains with numbers 0 and 1
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(subdomains, 0)
subdomain_1.mark(subdomains, 1)
# plot(subdomains, title='subdomains')
# plt.show()
# fill cell value in k
V0 = FunctionSpace(mesh,"DG",0)
k = Function(V0)
k_values = [1.5, 50] # values of k in the two subdomains
print(len(subdomains.array()))
for cell_no in range(len(subdomains.array())):
subdomain_no = subdomains.array()[cell_no]
k.vector()[cell_no] = k_values[subdomain_no]
# print('k degree of freedoms:', k.vector().get_local())
# define corrector equation
u, m = TrialFunctions(X)
v, r = TestFunctions(X)
Fc = u*r*dx + + m*v*ds
# Fl = inner(nabla_grad(u) + Constant((1,0)), k*nabla_grad(v))*dx
Fl = nabla_grad(u) * k * nabla_grad(v) * dx + k * nabla_grad(v) * dx
F = Fc + Fl
# Compute first order corrector
x = Function(X)
solve(F == 0, x)
uh, mh = x.split()
plot(uh)
plot(mh)
plt.show()
But unfortunately, I got this error in Fl = nabla_grad(u) * k * nabla_grad(v) * dx + k * nabla_grad(v) * dx.
Exception has occurred: UFLException
Invalid ranks 1 and 1 in product.
File "/mnt/c/Users/Hanyue/Desktop/periodic_homo/test.py", line 65, in <module>
Fl = nabla_grad(u) * k * nabla_grad(v) * dx + k * nabla_grad(v) * dx
~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~
ufl.log.UFLException: Invalid ranks 1 and 1 in product.
Is anyone know how to fix it? To be honest, I don’t really understand the rank here.