Solution not being updated at solving Picard iteration of coupled problem (Question changed)

variational formulation
1

Hi, I’m quite noob in Fenics 2019.1.0 and I’m having serious problems in the implementation of a Picard iteration method for a nonlinear coupled mixed-primal problem that goes as follow:

Find (\sigma,u,\phi) \in \mathbb{X} \times \mathbf{M} \times H such that

\int_\Omega \frac{1}{\mu(\phi)} \sigma^d : \tau^d + \int_\Omega u : \operatorname{div} \tau = \int_{\Gamma} \tau \nu \cdot u_D \quad \forall \tau \in \mathbb{X}
\int_\Omega v : \operatorname{div} \sigma = - \int_\Omega \phi f \cdot v \quad \forall v \in \mathbf{M}
\int_\Omega \vartheta(|\nabla \phi|) \nabla \phi \cdot \nabla \psi - \int_\Omega \phi u \cdot \nabla \psi = \int_\Omega \gamma(\phi) k \cdot \nabla \psi + \int_\Omega g \psi \quad \forall \psi \in H

where k is a constant vector, f is a vector function, g is a scalar function.

My idea is to iterate between the first two equations and the third one. That is, start with a given \phi_N, solve the first two equations obtaining \sigma_N and u_N. Then, replace u_N and \phi_N in the third equation (\phi_N just inside \gamma) and solve it for a new \phi_N. Then, repeat until some convergence criteria is met.

It seems that the solutions is not updating properly because the plots are empty.

The importants part of my unworking code are

# Load mesh from file
mesh = Mesh() #malla
hdf = HDF5File(mesh.mpi_comm(), "mallas/cuadrado.h5", "r")
hdf.read(mesh, "/mesh", False)

# Referencias de los dominios
domains = MeshFunction("size_t", mesh, dim) #CellFunction("size_t", mesh) #dominio
hdf.read(domains, "/domains")
#Referencias de los lados
ref_sides = MeshFunction("size_t", mesh, dim-1) #FacetFunction("size_t", mesh) #caras
hdf.read(ref_sides, "/bc_markers")

ds = Measure('ds', domain=mesh, subdomain_data=ref_sides)
n = FacetNormal(mesh)

# Problem data
def mu(s):
  return pow(1 - c*s,-2)
def gamma(s):
  return c*s*pow(1-c*s,2)
def theta(s):
  return m1+m2*pow(1+s*s,0.5*m3-1)

#Auxiliar function
def norma(x):
  return sqrt(pow(x[0],2)+pow(x[1],2))
def deviator(A):
  return A - 0.5*tr(A)*Identity(2)

# Define function spaces 
X = VectorElement("RT",mesh.ufl_cell(), degree) 
M = VectorElement("DG",mesh.ufl_cell(), degree)
H = FiniteElement("Lagrange",mesh.ufl_cell(),degree + 1)

#Mixed spaces
WS = FunctionSpace(mesh,MixedElement([X,M]))
WStilde = FunctionSpace(mesh,H)

# Define trial and test functions
(sigma,u) = TrialFunctions(WS)
su = Function(WS)
phi = TrialFunction(WStilde)

(tau,v) = TestFunctions(WS)
psi = TestFunction(WStilde)

#Define iteration (previous step) functions
suN = Function(WS)
sigmaN,uN = split(suN)  

phiN = Function(WStilde)

#Null boundary conditions for phi
bcphi = DirichletBC(WStilde,Constant(0), ref_sides, 0)
bcs = [bcphi]

for i in range(10):
  #S(phi) = (sigma,u)
  lhss1 = inner(deviator(sigma)/mu(phiN),deviator(tau))*dx + inner(u,div(tau))*dx + inner(v,div(sigma))*dx 
  rhss1 = inner(tau*n,u_exact)*ds(0) - inner(f,phiN*v)*dx

  solve(lhss1 == rhss1,suN)
  sigmaN,uN = split(suN)

  #Stilde(phi,u) = phitilde
  lhss2 = inner(theta(norma(grad(phiN)))*grad(phi),grad(psi))*dx - inner(phi*uN,grad(psi))*dx
  rhss2 = inner(gamma(phiN)*kv,grad(psi))*dx + inner(g,psi)*dx

  phiN1 = Function(WStilde)
  solve(lhss2 == rhss2,phiN1,bcs)
  phiN.assign(phiN1)
  print(i)

import matplotlib.pyplot as plt
plot(sigmaN)
plt.show()
plot(uN)
plt.show()
plot(phiN)
plt.show()

What am I doing wrong?

I must notice that I’m not imposing any boundary condition in the first variational problem (although I should impose \int_\Omega \operatorname{tr} \sigma = 0, but I don’t know how to do that)

2

Hi, I do not address all your issues, but maybe it will help to start.

From the first look it seems like you have not defined any TrialFunctions and built your left hand sides with just Functions and TestFunctions. So, the error says that the left hand side is expected to be of the form of an operator (A: dolfin.cpp.la.GenericLinearOperator), instead it detects that you have passed a Vector form (Invoked with: <dolfin.cpp.la.Vector object at 0x7f4addcb67d8>). Refer to the FEniCS demos and book for the correct way to implement.

I believe the below Functions should have been defined as TrialFunctions.

su = Function(WS)
sigma , u = split(su)

phi = Function(WStilde)

Sid

3

Thanks for your quick response. I followed your suggestion and wrote

(sigma,u) = TrialFunctions(WS)

phi = TrialFunctions(WStilde)

(tau,v) = TestFunctions(WS) 
psi = TestFunctions(WStilde)

Nevertheless, now when I construct the variational problem here

lhss2 = inner(theta(norma(grad(phi)))*grad(phi),grad(psi))*dx - inner(phi*uN,grad(psi))*dx

I get the following error

Traceback (most recent call last):
  File "mixed-primal.py", line 184, in <module>
    lhss2 = inner(theta(norma(grad(phi)))*grad(phi),grad(psi))*dx - inner(phi*uN,grad(psi))*dx
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 380, in grad
    f = as_ufl(f)
  File "/usr/lib/python3/dist-packages/ufl/constantvalue.py", line 472, in as_ufl
    " to any UFL type." % str(expression))
ufl.log.UFLValueError: Invalid type conversion: (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 6), FiniteElement('Lagrange', triangle, 2)), 1, None),) can not be converted to any UFL type.
4

There seem to be errors in the implementation. I would suggest going through the tutorials and approaching the problem step by step.
For example in the iteration loop: In the 1st problem you need to compute for a solution which would be a Function object and not a TrialFunction. So you should pass uN instead of su when asking to solve(because you have changed su to a TrialFunction). After solving the 1st problem you have to split the solution to form the new sub-functions (since it is MixedFunction) and then you can proceed with updating/assigning the solution to be used as a ‘known’ in the 2nd problem, and so on.

It would be difficult to read your question this way. People who answer regularly prefer a Minimum working example to replicate your issue (see this) sure someone with more expertise will pass by to answer in detail.

Sid

5

Ok, I will add some mathematical background explaining my desire and will reduce the code.

6

I have just modified the code giving some math background. Please let me know if it is more understandable now.

7

Now it is compiling but no solution is shown in the plots. I reformulated the question.