Hi everyone!

I’m new in FEniCS and have some troubles in using of the AdaptiveNonlinearVariationalSolver. From some point my solution stop to update with new values. Here is my code

######################################################################################
# Script for calculation of particle shrinkage because of boundary energy minimization
# Differential equation has look
# \dot{\phi} = - M \left( \alpha\Delta\phi - \omega \phi (1 - \phi)  (1 - 2\phi) \right)
######################################################################################

from fenics import *
import numpy as np
import time

## Input parameters for calculation

start = time.time()

L = 100 # Domain size in x and y, microns
alpha = 2 # Gradient energy pre-factor
omega = 1 # Double-well potential heigth
M_phi = 6 # Phase-field mobility
r_0 = 50 # initial particle radius
n_delta = 1 # numbers of elements per width of diffusive zone
d_t = 1 # time increment
t_end = 100 # time of ending of solution

# For mesh building and initiation of phi_0 field it is crucial to calculate boundary width, i.e
delta = 4.0*np.sqrt(alpha/omega)
print('Diffusive zone width:',delta)

# uniform mesh for this task
N_el = np.int64(L * n_delta / delta) # number of elements in calculation
print('Number of elements in 1D:',N_el)
mesh = RectangleMesh(Point(0,0),Point(L,L),N_el,N_el) # Generation of mesh
V = FunctionSpace(mesh,'P',1) #Function space on created mesh

# Boundary conditions
tol = 1E-10 #tolerance value for boundary searching
def boundary_Dir(x, on_boundary):
return (near(x[0],L,tol) and on_boundary) or (near(x[1],L,tol) and on_boundary) # Boundary of the Dirichlet condition
# Boundary on Neumann conditions should be not added because they have no influence of energy
bc = DirichletBC(V,Constant(1), boundary_Dir) #initialization of BC
# Initial conditions
phi_0 = Expression('0.5*(1 + tanh(0.0625*delta*(sqrt(x[0]*x[0]+x[1]*x[1])-r_0)))',
degree=1,r_0=r_0,delta=delta) # particle on the edge of corner in the mesh
phi_0_int = interpolate(phi_0,V) #interpolation of existing function on the mesh

# Defining of functions
phi = Function(V)
var_phi = TestFunction(V)

phi_i = phi_0_int # phi_i - order parameter at previous test
F_dw = M_phi*omega*d_t*phi*(1-phi)*(1-2*phi) # Fraction of double-well potential
F_td = phi - phi_i  #Fraction of time discretization
F = ( F_int + (F_dw + F_td)*var_phi) * dx # Full variational formulation

# Defining a Jacobian of non-linear form
dphi = TrialFunction(V)
J = derivative(F,phi,dphi)

# Defining of goal function and tolerance value for mesh adaption
M = phi * dx
tol_opt = 1e-4 * L**2

t = 0

vtkfile = File('particle-shrinkage_adaptive/solution.pvd') # initialization of file variable
vtkfile << (phi_i,t)

t = d_t

# Calculation itself
while t <= t_end :

#Solution of equation
problem = NonlinearVariationalProblem(F, phi, bc,J)
solver.solve(tol_opt)

print(phi_i.vector() == phi.leaf_node().vector())
if (phi_i.vector() == phi.leaf_node().vector()).all():
phi_i.assign(phi.root_node())
else:
phi_i.assign(phi.leaf_node()) #Update data from current step
print('time=',t)
vtkfile << (phi_i,t) # writing solution in file
t += d_t

end = time.time()
print('Time for calculation = ',(end-start)/60,' min')

Have you found the solution to your problem?
Because with a Nonlinear solver for time dependent, i find an accurate solution but with the adaptive solver, the solution is not updated after some timesteps.

@RossMech, @AgisM, where you’ve got

J = derivative(F,phi,dphi)


I’ve used

F = action (F, phi)
J = derivative(F,phi,dphi)


although I’m using NonlinearVariationalSolver not AdaptiveNonlinearVariationalSolver.

Also, where you’ve used Function phi to define your form, I’m using TrialFunction phit (i.e. dphi in your case). I think I need to use that action statement to link my solution phi to the phit or pdhi in the form.

I’m not certain that’s what you need, but it’s something to try.

I haven’t find any solution of this problem now.

Thank you very much!

I have found the solution myself for my problem.
I did not remember which was my question but for some cases initialization was the problem and in others you had to get the functionspace from the refined variable!