# Unable to access vector of degrees of freedom

**URL:** <https://fenicsproject.discourse.group/t/unable-to-access-vector-of-degrees-of-freedom/15601>\
**Category:** General\
**Created:** [August 21, 2024, 7:25pm UTC](https://fenicsproject.discourse.group/t/unable-to-access-vector-of-degrees-of-freedom/15601 "2024-08-21T19:25:34Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![sinhon2030](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sinhon2030/32/7442_2.png) [@sinhon2030](https://fenicsproject.discourse.group/u/sinhon2030)\
**Post date:** [August 21, 2024, 7:25pm UTC](https://fenicsproject.discourse.group/t/unable-to-access-vector-of-degrees-of-freedom/15601/1 "2024-08-21T19:25:34Z")

</div>

I am trying to model a phase field model (allen cahn) being advected by a viscoelastic oldroyd-b model. Here is the entire code. Could you let me know the mistake i am making?

```auto
from fenics import *
import numpy as np

#mesh and function space
mesh = UnitSquareMesh(64, 64)
V = VectorFunctionSpace(mesh, 'P', 2)#for velocity u
Q = FunctionSpace(mesh, 'P', 1) #for pressure p
S = TensorFunctionSpace(mesh, 'P', 1) #for stress tensor tau
P = FunctionSpace(mesh, 'P', 1) #for phase field phi

#Define current solution functions
W = FunctionSpace(mesh, MixedElement([V.ufl_element(), Q.ufl_element(), S.ufl_element()]))
w = Function(W) # Current solution (u, p, tau)
u, p, tau = split(w)
phi = Function(P) # Current solution (phi)

# Define previous solution functions
u_n = Function(V)
p_n = Function(Q)
tau_n = Function(S)
phi_n = Function(P)

# Define test functions
v, q, chi = TestFunctions(W)
psi = TestFunction(P)

#define material parameters
gamma = 0.1 # Adjusted gamma to control the phase field evolution
def eta_s(phi):
    return 1.0 + 0.5 * (1.0 + phi) # Tumor has higher eta_s, phi=1
def eta_p(phi):
    return 1.0 + 0.5 * (1.0 - phi) # Non-tumor has higher eta_p, phi=-1
def lambda_(phi):
    return 1.0 + 0.5 * (1.0 - phi) # Non-Tumor has higher lambda_, phi=1
def rho(phi):
    return 1.0 + 0.5 * (1.0 - phi) # Non-Tumor has higher rho, phi=1
def alpha(phi):
    return 0.5 * (1 + phi) # Tumor has higher alpha, phi=1

# Body force (proportional to velocity)
f_body = -alpha(phi)*u # Corrected the sign of the body force

# Define the stress tensor
D = sym(grad(u))
sigma = -p * Identity(len(u)) + 2 * eta_s(phi) * D + tau

# Variational form for polymeric stress evolution
F_tau = (inner((tau - tau_n) / dt, chi) * dx 
         + inner(dot(u, nabla_grad(tau)), chi) * dx 
         - inner(dot(tau, grad(u)), chi) * dx 
         - inner(dot(grad(u).T, tau), chi) * dx 
         - (2 * eta_p(phi) / lambda_(phi)) * inner(D, chi) * dx 
         + (1 / lambda_(phi)) * inner(tau, chi) * dx)

# Variational form for momentum balance (mixed formulation)
F_u = (rho(phi) * inner((u - u_n) / dt, v) * dx 
       + rho(phi) * inner(dot(u, nabla_grad(u)), v) * dx 
       + inner(sigma, grad(v)) * dx # Stress tensor applied to gradient of test function
       - dot(f_body, v) * dx) # Body force term

# variational form for the Allen-Cahn phase field
epsilon = 0.01 # Interface thickness parameter
F_phi = (inner((phi - phi_n) / dt, psi) * dx 
         + inner(dot(u, grad(phi)), psi) * dx 
         + gamma * epsilon**2 * inner(grad(phi), grad(psi)) * dx 
         - gamma * inner(phi**3 - phi, psi) * dx)

# Define the Jacobians for the nonlinear solvers
J_tau = derivative(F_tau, w.sub(2))
J_u = derivative(F_u, w)
J_phi = derivative(F_phi, phi)

# Define solvers with adjusted parameters
newton_solver_parameters = {
    "nonlinear_solver": "newton",
    "newton_solver": {
        "relative_tolerance": 1e-6,
        "absolute_tolerance": 1e-10,
        "maximum_iterations": 100,
        "relaxation_parameter": 1.0,
    }
}
# Define solvers
tau_problem = NonlinearVariationalProblem(F_tau, w.sub(2), J=J_tau)
tau_solver = NonlinearVariationalSolver(tau_problem)
tau_solver.parameters.update(newton_solver_parameters)

u_problem = NonlinearVariationalProblem(F_u, w, J=J_u)
u_solver = NonlinearVariationalSolver(u_problem)
u_solver.parameters.update(newton_solver_parameters)

phi_problem = NonlinearVariationalProblem(F_phi, phi, J=J_phi)
phi_solver = NonlinearVariationalSolver(phi_problem)
phi_solver.parameters.update(newton_solver_parameters)

# Initialize the previous time step solutions
u_n.assign(Constant((0.1, 0.1))) # Small initial velocity
p_n.assign(Constant(0.0)) # Initial pressure
tau_n.assign(Constant(((0.01, 0.01), (0.01, 0.01)))) # Small initial stress

# Initialize the phase field with the generalized tanh function
class InitialCondition(UserExpression):
    def eval(self, values, x):
        R0 = 0.2 # Critical radius
        epsilon_ = 0.01 # Controls the width of the transition
        r = np.sqrt((x[0] - 0.5)**2 + (x[1] - 0.5)**2)
        values[0] = np.tanh((R0 - r) / epsilon_)
        
    def value_shape(self):
        return ()

phi.interpolate(InitialCondition(degree=1))
phi_n.assign(phi)

# Parameters
dt = 0.005 # Smaller time step to improve convergence
num_steps = 200 # Increased number of steps due to reduced dt
t = 0

phi_values = []
for n in range(num_steps):
    
    t += dt
    print(n,t)

    tau_solver.solve()

    
    u_solver.solve()

   
    phi_solver.solve()

   
    phi_values.append(phi.copy(deepcopy=True))

    u_sol, p_sol, tau_sol = w.split()

  
    u_n.assign(u_sol.copy(deepcopy=True))   
    p_n.assign(p_sol.copy(deepcopy=True))   
    tau_n.assign(tau_sol.copy(deepcopy=True))  
    phi_n.assign(phi.copy(deepcopy=True))

```

The error i got is this: ---------------------------------------------------------------------------  
RuntimeError Traceback (most recent call last)  
 in \<cell line: 2\>()  
5 print(n,t)  
6  
----\> 7 tau\_solver.solve()  
8  
9

RuntimeError:

\*\*\* -------------------------------------------------------------------------  
\*\*\* DOLFIN encountered an error. If you are not able to resolve this issue  
\*\*\* using the information listed below, you can ask for help at

* * *

\*\*\* [https://fenicsproject.discourse.group/](https://fenicsproject.discourse.group/)

* * *

\*\*\* Remember to include the error message listed below and, if possible,  
\*\*\* include a _minimal_ running example to reproduce the error.

* * *

\*\*\* -------------------------------------------------------------------------  
\*\*\* Error: Unable to access vector of degrees of freedom.  
\*\*\* Reason: Cannot access a non-const vector from a subfunction.  
\*\*\* Where: This error was encountered inside Function.cpp.  
\*\*\* Process: 0

* * *

\*\*\* DOLFIN version: 2019.2.0.dev0  
\*\*\* Git changeset: 1c52e837eb54cc34627f90bde254be4aa8a2ae17  
\*\*\* -------------------------------------------------------------------------

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [August 22, 2024, 5:08am UTC](https://fenicsproject.discourse.group/t/unable-to-access-vector-of-degrees-of-freedom/15601/2 "2024-08-22T05:08:28Z")

</div>

The issue is that `w.sub(2)` creates a view into the function in the space `W` which has three components, while you’d want to create a function in the space `S`.

I made the following changes and are able to get the code running

```auto
...

# Variational form for polymeric stress evolution
tau_guess = Function(S)
chi_TEST_IN_S_RATHER_THAN_W = TestFunction(S)
F_tau = inner((tau_guess - tau_n) / dt, chi_TEST_IN_S_RATHER_THAN_W) * dx
# of course, choose a better name than chi_TEST_IN_S_RATHER_THAN_W in your actual code ;)

...

# Define the Jacobians for the nonlinear solvers
J_tau = derivative(F_tau, tau_guess)

...

# Define solvers
tau_problem = NonlinearVariationalProblem(F_tau, tau_guess, J=J_tau)
tau_solver = NonlinearVariationalSolver(tau_problem)
tau_solver.parameters.update(newton_solver_parameters)

...

for n in range(num_steps):

    t += dt
    print(n,t)

    tau_solver.solve()
    assign(w.sub(2), tau_guess)

    ...

```

If you have further issues, please make sure to post the updated code upon integrating these changes.

---

<div class="post-metadata">

**Author:** ![sinhon2030](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sinhon2030/32/7442_2.png) [@sinhon2030](https://fenicsproject.discourse.group/u/sinhon2030)\
**Post date:** [August 22, 2024, 2:37pm UTC](https://fenicsproject.discourse.group/t/unable-to-access-vector-of-degrees-of-freedom/15601/3 "2024-08-22T14:37:48Z")

</div>

Thanks a lot for your help. Now, I could go past the tau\_solver, but now there is an error with the u\_solver. I am pasting both the revised code and the error. Please let me know if I need to modify the u\_solver. In general, I think I am confused due to the mixed formulation or it might just be due to convergence. Please let know what you think.

```auto
from fenics import *
import numpy as np

#sim parameters
dt = 0.005 # Smaller time step to improve convergence
num_steps = 200 # Increased number of steps due to reduced dt
t = 0

#material_parameters
gamma = 0.1 # Adjusted gamma to control the phase field evolution
def eta_s(phi):
    return 1.0 + 0.5 * (1.0 + phi) # Tumor has higher eta_s, phi=1
def eta_p(phi):
    return 1.0 + 0.5 * (1.0 - phi) # Non-tumor has higher eta_p, phi=-1
def lambda_(phi):
    return 1.0 + 0.5 * (1.0 - phi) # Non-Tumor has higher lambda_, phi=1
def rho(phi):
    return 1.0 + 0.5 * (1.0 - phi) # Non-Tumor has higher rho, phi=1
def alpha(phi):
    return 0.5 * (1 + phi) # Tumor has higher alpha, phi=1
# Create mesh and define function spaces
mesh = UnitSquareMesh(64, 64)
V = VectorFunctionSpace(mesh, 'P', 2)#for velocity u
Q = FunctionSpace(mesh, 'P', 1) #for pressure p
S = TensorFunctionSpace(mesh, 'P', 1) #for stress tensor tau
P = FunctionSpace(mesh, 'P', 1) #for phase field phi

# Define current solution functions
W = FunctionSpace(mesh, MixedElement([V.ufl_element(), Q.ufl_element(), S.ufl_element()]))
w = Function(W) # Current solution (u, p, tau)
u, p, tau = split(w)
phi = Function(P) # Current solution (phi)

# Define previous solution functions
u_n = Function(V)
p_n = Function(Q)
tau_n = Function(S)
phi_n = Function(P)

# Define test functions
v, q, chi = TestFunctions(W)
chi_test=TestFunction(S)
psi = TestFunction(P)

# Body force (proportional to velocity)
f_body = -alpha(phi)*u # Corrected the sign of the body force

# Define the stress tensor
D = sym(grad(u))
sigma = -p * Identity(len(u)) + 2 * eta_s(phi) * D + tau

# Variational form for polymeric stress evolution
tau_guess = Function(S)
F_tau = (inner((tau_guess - tau_n) / dt, chi_test) * dx 
         + inner(dot(u, nabla_grad(tau_guess)), chi_test) * dx 
         - inner(dot(tau_guess, grad(u)), chi_test) * dx 
         - inner(dot(grad(u).T, tau_guess), chi_test) * dx 
         - (2 * eta_p(phi) / lambda_(phi)) * inner(D, chi_test) * dx 
         + (1 / lambda_(phi)) * inner(tau_guess, chi_test) * dx)

# Variational form for momentum balance (mixed formulation)
F_u = (rho(phi) * inner((u - u_n) / dt, v) * dx 
       + rho(phi) * inner(dot(u, nabla_grad(u)), v) * dx 
       + inner(sigma, grad(v)) * dx # Stress tensor applied to gradient of test function
       - dot(f_body, v) * dx) # Body force term

# Corrected variational form for the Allen-Cahn phase field
epsilon = 0.01 # Interface thickness parameter
F_phi = (inner((phi - phi_n) / dt, psi) * dx 
         + inner(dot(u, grad(phi)), psi) * dx 
         + gamma * epsilon**2 * inner(grad(phi), grad(psi)) * dx 
         - gamma * inner(phi**3 - phi, psi) * dx)

# Define the Jacobians for the nonlinear solvers
J_tau = derivative(F_tau, tau_guess)
J_u = derivative(F_u, w)
J_phi = derivative(F_phi, phi)

# Define solvers with adjusted parameters
newton_solver_parameters = {
    "nonlinear_solver": "newton",
    "newton_solver": {
        "relative_tolerance": 1e-6,
        "absolute_tolerance": 1e-10,
        "maximum_iterations": 100,
        "relaxation_parameter": 1.0,
    }
}
# Define solvers
tau_problem = NonlinearVariationalProblem(F_tau, tau_guess, J=J_tau)
tau_solver = NonlinearVariationalSolver(tau_problem)
tau_solver.parameters.update(newton_solver_parameters)

u_problem = NonlinearVariationalProblem(F_u, w, J=J_u)
u_solver = NonlinearVariationalSolver(u_problem)
u_solver.parameters.update(newton_solver_parameters)

phi_problem = NonlinearVariationalProblem(F_phi, phi, J=J_phi)
phi_solver = NonlinearVariationalSolver(phi_problem)
phi_solver.parameters.update(newton_solver_parameters)

# Initialize the previous time step solutions
u_n.assign(Constant((0.1, 0.1))) # Small initial velocity
p_n.assign(Constant(0.0)) # Initial pressure
tau_n.assign(Constant(((0.01, 0.01), (0.01, 0.01)))) # Small initial stress

# Initialize the phase field with the generalized tanh function
class InitialCondition(UserExpression):
    def eval(self, values, x):
        R0 = 0.2 # Critical radius
        epsilon_ = 0.01 # Controls the width of the transition
        r = np.sqrt((x[0] - 0.5)**2 + (x[1] - 0.5)**2)
        values[0] = np.tanh((R0 - r) / epsilon_)
        
    def value_shape(self):
        return ()

phi.interpolate(InitialCondition(degree=1))
phi_n.assign(phi)

phi_values = []
for n in range(num_steps):
    
    t += dt
    print(n,t)

    tau_solver.solve()
    assign(w.sub(2), tau_guess)
    print("tau solved")

    
    u_solver.solve()
    print("u solved")

   
    phi_solver.solve()
    print("phi solved")

   
    phi_values.append(phi.copy(deepcopy=True))

    u_sol, p_sol, tau_sol = w.split()

  
    u_n.assign(u_sol.copy(deepcopy=True))   
    p_n.assign(p_sol.copy(deepcopy=True))   
    tau_n.assign(tau_sol.copy(deepcopy=True))  
    phi_n.assign(phi.copy(deepcopy=True))

```

This is the error I got: ---------------------------------------------------------------------------  
RuntimeError Traceback (most recent call last)  
 in \<cell line: 2\>()  
10  
11  
—\> 12 u\_solver.solve()  
13 print(“u solved”)  
14

RuntimeError:

\*\*\* Error: Unable to solve nonlinear system with NewtonSolver.  
\*\*\* Reason: Newton solver did not converge because maximum number of iterations reached.  
\*\*\* Where: This error was encountered inside NewtonSolver.cpp.  
\*\*\* Process: 0

* * *

\*\*\* DOLFIN version: 2019.2.0.dev0  
\*\*\* Git changeset: 1c52e837eb54cc34627f90bde254be4aa8a2ae17  
\*\*\* -------------------------------------------------------------------------

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [August 22, 2024, 4:36pm UTC](https://fenicsproject.discourse.group/t/unable-to-access-vector-of-degrees-of-freedom/15601/4 "2024-08-22T16:36:57Z")

</div>

That typically means that the initial guess is not good enough. With your domain specific knowledge you should come up with a good initial guess.

Some general suggestions are also discussed in [Default absolute tolerance and relative tolerance - #4 by nate](https://fenicsproject.discourse.group/t/default-absolute-tolerance-and-relative-tolerance/3829/4) by @nate.

---

<div class="post-metadata">

**Author:** ![sinhon2030](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sinhon2030/32/7442_2.png) [@sinhon2030](https://fenicsproject.discourse.group/u/sinhon2030)\
**Post date:** [August 23, 2024, 1:37am UTC](https://fenicsproject.discourse.group/t/unable-to-access-vector-of-degrees-of-freedom/15601/5 "2024-08-23T01:37:21Z")

</div>

Thanks a lot for your help.