Hello,
Following your advices I have tried to apply the cahn-hilliard demo to my problem with four constitutive variables(the solid displacement, the gas pressure and the temperature fields). I have corrected all the "formulation" errors, but now I have a problem when I launch the program, here is an extract of my code, I have removed the parts where I define my boundary conditions and the constants/parameters in order to reduce it as much as I could:
import random
from dolfin import *
import numpy
cell = triangle
mesh = Mesh("bonmesh.xml")
V1 = VectorFunctionSpace(mesh, 'CG', 2)
V2 = FunctionSpace(mesh, 'CG', 1)
V3 = FunctionSpace(mesh, 'CG', 1)
V4 = FunctionSpace(mesh, 'CG', 1)
W = MixedFunctionSpace([V1, V2, V3, V4])
[...] define boundary conditions
# Class representing the intial conditions
class InitialConditions(Expression):
def __init__(self):
random.seed(2 + MPI.process_number())
def eval(self, values, x):
values[0] = 0.0
values[1] = 0.0
values[2] = 1.0
values[3] = 40.0
values[4] = 40.0
def value_shape(self):
return (5,)
# Class for interfacing with the Newton solver
class ECRTCEquation(NonlinearProblem):
def __init__(self, a, L):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.reset_sparsity = True
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A, reset_sparsity=self.reset_sparsity)
self.reset_sparsity = False
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"
# Define trial and test functions
du = TrialFunction(W)
v, q, tS, tG = TestFunctions(W)
# Define functions
u = Function(W) # current solution
u0 = Function(W) # solution from previous converged step
# Split mixed functions
duS, dpGR, dthetaS, dthetaG = split(du)
uS, pGR, thetaS, thetaG = split(u)
uS0, pGR0, thetaS0, thetaG0 = split(u0)
# Create intial conditions and interpolate
u_init = InitialConditions()
u.interpolate(u_init)
u0.interpolate(u_init)
[...] define constants/parameters
# Weak statement of the equations
L0 = rhoSR(thetaS)*nS(uS)*inner(b, v)*dx+rhoGR(thetaG, pGR)*nG(uS)*inner(b, v)*dx-\
inner(sigmaS_E(uS, thetaS), grad(v))*dx+inner(nS(uS)*pGR*(thetaS/thetaG)*I, grad(v))*dx+\
inner(nG(uS)*pGR*I, grad(v))*dx-inner(TG_E(uS), grad(v))*dx
L1 = (1/muGR)*rhoGR(thetaG, pGR)*rhoGR(thetaG, pGR)*inner(dot(KS, b), grad(q))*dx-\
nG(uS)*(1/(RGR*thetaG))*DpGR(uS, pGR)*q*dx+\
nG(uS)*(1/(RGR*thetaG*thetaG))*pGR*DthetaG(uS, thetaG)*q*dx-\
(1/muGR)*rhoGR(thetaG, pGR)*inner(dot(KS, grad(pGR)), grad(q))*dx-\
(1/muGR*nG(uS))*rhoGR(thetaG, pGR)*pGR*(1-thetaS/thetaG)*inner(dot(KS, grad(nS(uS))), grad(q))*dx-\
nG(uS)*rhoGR(thetaG, pGR)*div((uS-uS0)/dt)*q*dx+\
DnS_S(uS, thetaS)*rhoGR(thetaG, pGR)*q*dx
L2=[...]
L3=[...]
L = L0+L1+L2+L3
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
# Create nonlinear problem and Newton solver
problem = ECRTCEquation(a, L)
solver = NewtonSolver("lu")
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
# Output file
file = File("output.pvd", "compressed")
# Step in time
t = 0.0
T = 80*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u.split()[0], t)
plot(u.split()[0])
L2 and L3 are built the same way as L0 and L1 (they are even longer and more complex).
If I launch the program like this, I get:
WARNING: The number of integration points for each cell will be: 1156
Consider using the option 'quadrature_degree' to reduce the number of points
How could I have 1156 points for each cell? Is there a problem with MixedFunctionSpace? Because if we add the integration points for each functionspace we are only at 15?(6+3+3+3?)
Therefore I have add: parameters["form_compiler"]["quadrature_degree"] = 4, and with that I finally obtain after a long time:
Traceback (most recent call last):
Form argument must be restricted.
Traceback (most recent call last):
File "mesh.py", line 294, in
solver.solve(problem, u.vector())
File "mesh.py", line 131, in F
assemble(self.L, tensor=b)
[...]
Memory error
Where do you think it comes from in my code? Is there an obvious error I have done?
Previously, I am not very entrusted in my initial conditions, because values[0] and values[1] represent uS0=(0,0), but if I put "values[0] = Constant((0.0, 0.0))" and thus "return (4,)" I have:
Traceback (most recent call last):
File "mesh.py", line 157, in
u.interpolate(u_init)
Error: Unable to interpolate function into function space.
Reason: Dimension 0 of function (4) does not match dimension 0 of function space (5).
Maybe should I try another solver?
So I go a bit round in circles, if someone wants more precisions and details about my problem I would be glad to answer, and in any case thank you in advance for any time and effort you will take in looking at my question.
Remi
