Phase field for sintering || solve converge

[Judy-22]

Hello MOOSE community,

I have tried to create a sintering model for 2 particles, I used all units in SI. I have problem about solve convergence.

I get this error when I launch the input file.
*** ERROR ***
The following error occurred in the object "TimeStepper", of type "ConstantDT".
Solve failed and timestep already at or below dtmin, cannot continue!

I think the problem is related to scaling, when I try to apply automatic scaling, I get this error.
*** ERROR ***
The following error occurred in the object "w", of type "MooseVariable".
VariableAllDoFMap not prepared for w . Check nonlocal coupling requirement for the variable.

Please find below the input file.

Thanks for help.

Judice, PhD student, ICB-COOM Lab.

_______________________________________________________________________________________
#   The phase field model to simulate sintering of two particles in 2D with 
#   constant temperature.

[GlobalParams]        # Automatically define order parameters 
  var_name_base = eta
  op_num = 2
[]

[Mesh]                # Generate mesh (7/3D X 4/3D) 
  type = GeneratedMesh
  dim = 2
  nx = 20             # n°s element in x. 
  ny = 10             # n°s element in y.
  xmin = 0  
  xmax = 70.0e-6      # m
  ymin = 0
  ymax = 40.0e-6      # m
  elem_type = QUAD4
  uniform_refine = 1
[]

[Variables]
  [./c]               # particle concentration (conserved variable)
    order = FIRST
    family = LAGRANGE
  [../]
  [./w]               # functional derivative of the total free energy  
    order = FIRST
    family = LAGRANGE
  [../]
  [./PolycrystalVariables] # Automatically creates order parameter variables
  [../]
[]

[Kernels]
  [./c_res]
  # compute the residual therm of total free energy as function of c.
    type = SplitCHParsed
    variable = c
    f_name = F
    kappa_name = kappa_c
    args = 'eta0 eta1'
    w = w
  [../]
  [./w_res]
  # compute the residual therm of functional free energy as function of w and M.
    type = SplitCHWRes
    variable = w
    mob_name = M
  [../]
  [./time]
  # compute the time derivative therm of w and c.
    type = CoupledTimeDerivative
    variable = w
    v = c
  [../]
  [./motion]
  # kernel to compute the rigid body motion.
    type = MultiGrainRigidBodyMotion
    variable = w
    c = c
    v = 'eta0 eta1'
    grain_force = grain_force
    grain_tracker_object = grain_center
    grain_volumes = grain_volumes
  [../]
  [./eta0_dot]
  # Compute the kernel that implements a simple time derivative of eta.
    type = TimeDerivative
    variable = eta0
  [../]
  [./vadv_eta0]
  # compute advection velocity of particle with index 0.
    type = SingleGrainRigidBodyMotion
    variable = eta0
    c = c
    v = 'eta0 eta1'
    grain_force = grain_force
    grain_tracker_object = grain_center
    grain_volumes = grain_volumes
    op_index = 0
  [../]
  [./acint_eta0]
  # Compute gradient energy Allen-Cahn Kernel
    type = ACInterface
    variable = eta0
    mob_name = L
    #args = 'T'
    kappa_name = kappa_eta
  [../]
  [./acbulk_eta0]
  # Compute Allen-Cahn Kernel that uses a DerivativeMaterial Total Free Energy
    type = AllenCahn
    variable = eta0
    mob_name = L
    f_name = F
    args = 'c eta1'
  [../]
  [./eta1_dot]
  # Compute the kernel that implements a simple time derivative of eta.
    type = TimeDerivative
    variable = eta1
  [../]
  [./vadv_eta1]
  # compute advection velocity of particle with index 1.
    type = SingleGrainRigidBodyMotion
    variable = eta1
    c = c
    v = 'eta0 eta1'
    op_index = 1
    grain_force = grain_force
    grain_tracker_object = grain_center
    grain_volumes = grain_volumes
  [../]
  [./acint_eta1]
    type = ACInterface
    variable = eta1
    mob_name = L
    #args = 'T'
    kappa_name = kappa_eta
  [../]
  [./acbulk_eta1]
    type = AllenCahn
    variable = eta1
    mob_name = L
    f_name = F 
    args = 'c eta0'
  [../]
  
[]

[Materials]  
  [./constants]
    type = GenericConstantMaterial
    prop_names = 'kappa_c   kappa_eta     A      B   
                    Dvol      Dvap       Dsuf    Dgb
                    Qvol      Qvap       Qsuf    Qgb
                    KB         Vm        Mgb     Qm        wGB      T '
    
    prop_values = '8.46e-06     2.385e-06    7.1667e+06   0.3533e+06
                   2.0e-04      0.0          0.4          0.127e-04
                   2.6          0.0          2.21         0.58
                   8.61733e-05  1.21e-29     5.53e-08     0.171  5.0e-06   1200'
                   
    # SI units    J/m          J/m          J/m³         J/m³
    # SI units    m³/s         m³/s         m³/s         m³/s
    # SI units    eV           eV           eV           eV
    # SI units    eV/K         m³/mol       m⁴/J.s       eV    m    K
  [../]
  [./free_energy]            # Total free energy function of c and eta.
    type = DerivativeParsedMaterial
    f_name = F
    args = 'c eta0 eta1'          # Must be changed as op_num changes.
    function = '(A*c^2*(1-c)^2 + B*(c^2+6*(1-c)*(eta0^2+eta1^2)-4*(2-c)*(eta0^3+eta1^3)
                + 3*(eta0^2+eta1^2)^2))'
    material_property_names = 'A  B'
    derivative_order = 2
  [../]
  [./CHmobility]    # Cahn and Hilliard Mobility incluing all diffusin path (m^5/J/s)
    type = DerivativeParsedMaterial
    f_name = M
    args = 'c eta0 eta1'
    function = '(Vm/(KB*T*1.6022e-19))*(Dvol*exp(-Qvol/(KB*T))*c^3*(10 - 15*c + 6*c^2)
                + Dvap*exp(-Qvap/(KB*T))*(1 - (c^3*(10 - 15*c + 6*c^2))) 
                + Dsuf*exp(-Qs/(KB*T))*(1 - c) + Dgb*exp(-Qgb/(KB*T))*(eta0*eta1))'
    material_property_names = 'Dvol Qvol Dvap Qvap Dsuf Qs Dgb Qgb KB Vm d2F:=D[F(c),c,c] T'
    derivative_order = 1
  [../]
  [./L]             # Allen and Cahn Mobility or GBMobility (m^3/J/s)
    type = GBEvolution
    T = 1200
    wGB = 5.0e-06
    GBmob0 = 5.53e-8
    Q = 0.171
    GBenergy = 1.06
    molar_volume = 7.012e-06
    length_scale = 1.0
    time_scale   = 1.0
  [../]
  [./force_density]  # Force between particles caused by material density.
    type = ForceDensityMaterial
    c = c
    etas ='eta0 eta1'
  [../]
[]

[BCs]
  [./bcs]
    #zero flux BC
    type = NeumannBC
    value = 0
    variable = c
    boundary = '0 1 2 3'
  [../]
[]

[AuxVariables]
  [./bnds]
  # Auxivariables that define the location of grain boundaries in a polycrystalline sample
  [../]
  [./df00]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./df01]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./df10]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./df11]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./unique_grains]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./var_indices]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./centroids]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./bnds]
  # Calculate location of grain boundaries in a polycrystalline sample
    type = BndsCalcAux
    variable = bnds
    var_name_base = eta
    op_num = 2
    v = 'eta0 eta1'
  [../]
  [./df01]
    type = MaterialStdVectorRealGradientAux
    variable = df01
    index = 0
    component = 1
    property = force_density
  [../]
  [./df11]
    type = MaterialStdVectorRealGradientAux
    variable = df11
    index = 1
    component = 1
    property = force_density
  [../]
  [./df00]
    type = MaterialStdVectorRealGradientAux
    variable = df00
    index = 0
    component = 0
    property = force_density
  [../]
  [./df10]
    type = MaterialStdVectorRealGradientAux
    variable = df10
    index = 1
    component = 0
    property = force_density
  [../]
  [./unique_grains]
    type = FeatureFloodCountAux
    variable = unique_grains
    flood_counter = grain_center
    field_display = UNIQUE_REGION
    execute_on = timestep_begin
  [../]
  [./var_indices]
    type = FeatureFloodCountAux
    variable = var_indices
    flood_counter = grain_center
    field_display = VARIABLE_COLORING
    execute_on = timestep_begin
  [../]
  [./centroids]
    type = FeatureFloodCountAux
    variable = centroids
    execute_on = timestep_begin
    field_display = CENTROID
    flood_counter = grain_center
  [../]
[]

[ICs]
  [./ic_eta0]
  # initial condition for each order parameter 
    type = SmoothCircleIC
    variable = eta0
    x1 = 20e-6
    y1 = 20e-6
    radius = 15.0e-6
    invalue = 1.0
    outvalue = 0.0
    int_width = 5.0e-6
  [../]
  [./IC_eta1]
    type = SmoothCircleIC
    variable = eta1
    x1 = 50e-6
    y1 = 20e-6
    radius = 15.0e-6
    invalue = 1.0
    outvalue = 0.0
    int_width = 5.0e-6
  [../]
  [./ic_c]
  # initial condition for conserved parameter 
    type = SpecifiedSmoothCircleIC
    invalue = 1.0
    outvalue = 0.0
    int_width = 5.0e-6
    x_positions = '20.0e-6 50.0e-6'
    y_positions = '20.0e-6 20.0e-6'
    z_positions = '0.0 0.0'
    radii       = '15.0e-6 15.0e-6'
    3D_spheres = false
    variable = c
    block = 0
  [../]
[]

[VectorPostprocessors]
  [./forces]
    type = GrainForcesPostprocessor
    grain_force = grain_force
  [../]
  [./grain_volumes]
    type = FeatureVolumeVectorPostprocessor
    flood_counter = grain_center
    execute_on = 'initial timestep_begin'
  [../]
[]

[Postprocessors]
  [./average_c]
    type = ElementAverageValue
    variable = c
    execute_on = 'initial timestep_end'
  [../]
[]

[UserObjects]
  [./grain_center]
    type = GrainTracker
    outputs = none
    compute_var_to_feature_map = true
    execute_on = 'initial timestep_begin'
  [../]
  [./grain_force]
    type = ComputeGrainForceAndTorque
    execute_on = 'linear nonlinear'
    grain_data = grain_center
    force_density = force_density
    c = c
    etas = 'eta0 eta1'
  [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_ksp_type
                         -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm      31                  preonly
                         ilu          2'
  l_tol = 1e-05
  nl_max_its = 30
  l_max_its = 30
  nl_rel_tol = 1e-07
  nl_abs_tol = 1e-09
  start_time = 0.0
  end_time = 0.01
  dt = 0.00125
[]

[Outputs]
  exodus = true
[]

[SudiptaBiswas]

Judice, You can not use dimensions in the order of 1e-6 in the finite element model. It's not numerically tractable. Try changing your length scale to micrormeter. You will have to define your domain and parameterize your material properties accordingly.

[HPL001]

@SudiptaBiswas Excuse me for taking the liberty, but I would like to know how to modify the revision scale you pointed out. I only found time_scale and length_scale in the input file, no energy scale. I am not sure if I can just modify here. Looking forward to your reply.

[Judy-22]

@HPL001,
Compatibility between magnitudes and scales must be taken into account. For your example, you don't need a scale, so it's equal to 1 for both.

[HPL001]

@Judy-22 Thank you for your reply. In fact, I copied your input file above. I have encountered the problem you mentioned above, but I still have no way to solve it now, because I don't understand the problem. In addition, regarding the surface diffusion coefficient here, the coefficient before it is mentioned as c* (1-c) in the literature, but it is all (1-c) in your document, and the derivative is also set to 1. I am not clear about the intention here, please give me some guidance. I look forward to your reply. It will be very meaningful to me. Thank you.

[Judy-22]

@HPL001, I developed some C++ files for the material properties class, which helped me deal with scales, material equations, ...

Maybe I can help, send me an e-mail with more details and to start a good discussion. [email protected]

[HPL001]

@Judy-22 Thank you for your help. I sent an email to you before, but I did not receive a reply. I am not sure whether you have received it or not (maybe I made some mistakes in the process of sending the email). I am looking forward to your reply and guidance, thank you