"Solve did not converge"

[abc-hy]

I have some questions about "Solve did not converge".
I want to set the mobility function as c dependent variable, which equals to M=abs(1-c^2), which is the absolute value of (1-c^2). However, my simulation always cannot converge, even though the time steps go downs to 10e(-14). When I do not set mobility as a function of c and just set M=1 or other constant, using the same parameters, it can work and can get results. I wonder what's the problem? Why can't it run when using c dependent mobility? Could you give me some advice about this?
Thank you very much.
The following is my input file.

[Mesh]
  type = GeneratedMesh
  dim = 2
  elem_type = QUAD4
  nx = 60
  ny = 24
  nz = 0
  xmin = 0
  xmax = 6.0
  ymin = 0
  ymax = 2.4
  zmin = 0
  zmax = 0
  uniform_refine = 3
  nemesis = true
  skip_partitioning = true
[]

[Variables]

  [./c]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = IsolatedBoundingBoxIC
      smaller_coordinate_corners = '0.0 0.3 0   0.0 1.5 0  4.0 0.3 0   4.0 1.5 0  2.0 0.0 0'
      larger_coordinate_corners =  '2.0 0.9 0   2.0 2.1 0  6.0 0.9 0   6.0 2.1 0  4.0 2.4 0'
      inside = '1.0 1.0 1.0 1.0 0.0'
      outside = -1.0
      int_width = 0
    [../]
  [../]

  [./w]
    order = FIRST
    family = LAGRANGE
  [../]

[]


[BCs]
  [./Periodic]
    [./cx]
      variable = c
      auto_direction = 'x,y'
    [../]
    [./wx]
      variable = w
      auto_direction = 'x,y'
    [../]
  [../]

[]


[Kernels]
  [./c_res]
    type = SplitCHParsed
    variable = c
    f_name = F
    kappa_name = kappa_c
    w = w
  [../]
  [./w_res]
    type = SplitCHWRes
    variable = w
    mob_name = M
  [../]
  [./time]
    type = CoupledTimeDerivative
    variable = w
    v = c
  [../]
[]


[Materials]
  [./mobility]
    type = DerivativeParsedMaterial
    f_name = M
    args = 'c'
    function = abs(1-c^2)
    derivative_order = 1
  [../]

  [./kappa_c]
    type = GenericConstantMaterial
    prop_names  = 'kappa_c'
    prop_values = '8.0e-4'
  [../]

  [./free_energy]
    # equivalent to `MathFreeEnergy`
    type = DerivativeParsedMaterial
    f_name = F
    args = 'c'
    function = (27/8)*c^2*(2-4*c^2+2*c^4)
  [../]
[]

[Postprocessors]
  [./step_size]             # Size of the time step
    type = TimestepSize
  [../]
  [./iterations]            # Number of iterations needed to converge timestep
    type = NumNonlinearIterations
  [../]
  [./nodes]                 # Number of nodes in mesh
    type = NumNodes
  [../]
  [./evaluations]           # Cumulative residual calculations for simulation
    type = NumResidualEvaluations
  [../]
  [./active_time]           # Time computer spent on simulation
    type = PerfGraphData
    section_name = "Root"
    data_type = total
  [../]
  [./c]            # output concentration file
    type = ElementAverageValue
    variable = 'c'
  [../]

[]

[VectorPostprocessors]
  [./x_direc]
   type =  LineValueSampler
    start_point = '0 1.2 0'
    end_point = '6 1.2 0'
    variable = 'c'
    num_points = 600
    sort_by =  id
  [../]
[]



[Preconditioning]
  # active = ' '
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  scheme = 'BDF2'
  #petsc_options = '-snes_mf'

  #Preconditioned JFNK (default)
  solve_type = 'NEWTON'

  petsc_options_iname = '-pc_type -ksp_grmres_restart -sub_ksp_type -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm      31                  preonly       lu           1'
  #petsc_options_iname = '-pc_type'
  #petsc_options_value = 'lu'

  l_max_its = 30 # maximum linear iterations
  l_tol = 1.0e-6 # 0.001 Linear Tolerance

  nl_max_its = 50 # maximum number of nonlinear iterations. exceed will cut dt.
  nl_rel_tol = 1.0e-9 # -8 nonlinear relative tolerance
  #num_steps = 2

  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-6
    cutback_factor = 0.67
    growth_factor = 1.5   # 1.6
    optimal_iterations = 10
    iteration_window = 2
  [../]

  [./Adaptivity]
    coarsen_fraction = 0.1
    refine_fraction = 0.7
    max_h_level = 3
    interval = 1
  [../]
[]

[Debug]
  show_var_residual_norms = true
[]

[Outputs]
  exodus = true
  console = true
  csv = true
  interval = 2
  checkpoint = true
  [./console]
    type = Console
    max_rows = 10
  [../]
[]

[jiangwen84]

Add line_search = 'none' in [Executioner] block.

[abc-hy]

The following is the last step's concentration data shown in the left column. I found that some of the concentration is a little bit less than -1. Is that the reason why the simulation cannot continue? So, all the concentration should be exactly between -1 and 1? In this case, how can I smooth the function?

NMC0_inconstant_deltaf0_R2_out_x_direc_0020.csv

[abc-hy]

Do I need to smooth the interface in the initial condition? Do you know how can I smooth the interface in the initial condition? Currently, in my input file, it is not diffuse interface, it is sharp interface...Is this the reason why I cannot run the simulation with that mobility?

[GiudGiud]

You can use the Bounds system to force it to stay between -1 and 1 if that is what is physical for this variable
https://mooseframework.inl.gov/source/auxkernels/ConstantBoundsAux.html

[abc-hy]

I have another question, which is about the interface width setting in the initial condition block.
Currently, in my input file, the interface width I set is 0,does this mean the model is sharp interface model? But I actually want to investigate the diffuse interface model, is it ok to set 0 interface width at the beginning? I think in the simulation, there will automatically generate an interface even if the interface width is set as 0, right? Do I still have to set an interface width in the initial condition block to make the model become diffuse interface model? If so, how large should I set for the interface width and will be maintain as constant during the simulation?

[jiangwen84]

@abc-hy

You should not set interface_width to be 0 in the ICs. It will cause issues in computing the phase-field value. The interface_width should be consistent with the kappa_c value. The actual width is typically determined by the interfacial energy.