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)
# Defining of variational task
phi_i = phi_0_int # phi_i - order parameter at previous test
F_int = M_phi*alpha*d_t*dot(grad(phi),grad(var_phi)) # Fraction of boundary energy
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 = AdaptiveNonlinearVariationalSolver(problem, M)
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')