Petsc DOF mismatch

[mhelel]

I have gotten the error Split 'upc' has size 3585 but the parent split expected size 3244. Make sure that you have non-overlapping complete sets for variables and blocks as well as consistency in sides/unsides, contacts/uncontacts, etc. while running the following code:

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmin = 0
  xmax = 3.0
  ymin = 0
  ymax = 1.0
  nx = 30
  ny = 10
  elem_type = QUAD9
[]

[Variables]
  [./velocity]
    order = SECOND
    family = LAGRANGE_VEC
  [../]

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

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

[Kernels]
    [./mass]
      type = INSMassVector
      variable = pressure
      velocity = velocity
    [../]

    [./momentum]
      type = INSMomentumVector
      variable = velocity
      pressure = pressure
    [../]
    [./marangoni]
      type = MarangoniForce
      variable = velocity
      concentration = concentration
    [../]

    [./diffusion]
      type = Diffusion
      variable = concentration
    [../]
[]

[BCs]
    [./NoSlip]
        type = VectorDirichletBC
        variable = velocity
        boundary = 'top bottom'
        values = '0 0 0'
    [../]

    [./inlet_p]
        type = DirichletBC
        variable = pressure
        boundary = left
        value = 1.0
    [../]
    [./outlet_p]
        type = DirichletBC
        variable = pressure
        boundary = right
        value = 0.0
    [../]
    [./inlet_c]
        type = DirichletBC
        variable = concentration
        boundary = left
        value = 1.0
    [../]
    [./outlet_c]
        type = DirichletBC
        variable = concentration
        boundary = right
        value = 0.0
    [../]
[]

[Preconditioning]
  [./FSP]
    type = FSP
    topsplit = 'upc'
    [./upc]
      splitting = 'up c'
      splitting_type = additive
    [../]
    [./up]
      splitting = 'u p' # 'u' and 'p' are the names of subsolvers
      splitting_type  = schur
      petsc_options_iname = '-pc_fieldsplit_schur_fact_type  -pc_fieldsplit_schur_precondition'
      petsc_options_value = 'full selfp'
    [../]
    [./u]
      vars = 'velocity'
      petsc_options_iname = '-pc_type -ksp_type -ksp_rtol'
      petsc_options_value = '     hypre gmres 1e-4'
    [../]
    [./p]
      vars = 'pressure'
      petsc_options_iname = '-pc_type -ksp_type -ksp_rtol'
      petsc_options_value = '   jacobi    gmres     1e-4'
    [../]
    [./c]
      vars = 'concentration'
      petsc_options_iname = '-pc_type -ksp_type -ksp_rtol'
      petsc_options_value = '   jacobi    gmres     1e-4'
    [../]
  [../]
[]

[Executioner]
  type = Steady
  solve_type = PJFNK
  nl_rel_tol = 1e-12
  nl_max_its = 6
  l_tol = 1e-6
  l_max_its = 300
[]

[Outputs]
  exodus = true
[]

The problem must come from adding the concentration variable since I have Code without this variable that works perfectly fine:

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmin = 0
  xmax = 3.0
  ymin = 0
  ymax = 1.0
  nx = 30
  ny = 10
  elem_type = QUAD9
[]

[Variables]
  [./velocity]
    order = SECOND
    family = LAGRANGE_VEC
  [../]

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

[Kernels]
    [./mass]
        type = INSMassVector
        variable = pressure
        velocity = velocity
    [../]

    [./momentum]
        type = INSMomentumVector
        variable = velocity
        pressure = pressure
    [../]
[]

[BCs]
    [./NoSlip]
        type = VectorDirichletBC
        variable = velocity
        boundary = 'top bottom'
        values = '0 0 0'
    [../]

    [./inlet_p]
        type = DirichletBC
        variable = pressure
        boundary = left
        value = 1.0
    [../]
    [./outlet_p]
        type = DirichletBC
        variable = pressure
        boundary = right
        value = 0.0
    [../]
[]

[Preconditioning]
    [./FSP]
        type = FSP
        topsplit = up

        [./up]
            splitting = 'u p'
            splitting_type = schur
            petsc_options_iname = '-pc_fieldsplit_schur_fact_type  -pc_fieldsplit_schur_precondition'
            petsc_options_value = 'full selfp'
        [../]
        [./u]
            vars = 'velocity'
            petsc_options_iname = '-pc_type -ksp_type -ksp_rtol'
            petsc_options_value = '     hypre gmres 1e-4'
        [../]
        [./p]
            vars = 'pressure'
            petsc_options_iname = '-pc_type -ksp_type -ksp_rtol'
            petsc_options_value = '   jacobi    gmres     1e-4'
        [../]
    [../]
[]

[Executioner]
  type = Steady
  solve_type = PJFNK
  nl_rel_tol = 1e-12
  nl_max_its = 6
  l_tol = 1e-6
  l_max_its = 300
[]

[Outputs]
  exodus = true
  [dof]
    type = DOFMap
    execute_on = 'initial'
  []
[]

I narrowed down a few things:

1. With DOFMap I found out that there are 660 DOFs for pressure, which makes sense, because there are 31 * 11 = 341 nodes with the Lagrange Basis of first oder, which should have 2 DOFs per node? With the DirichletBC on the left and right side I remove 22 DOFs. So this exactly matches.

2. From the error message and the DOFMap I know that adding the concentration variable adds 3585-2903 = 682 DOFs, which is exactly what I would expect that adding a variable does.

3. I also thought that the Kernel for Concentration might not work as intended or Boundary Conditions are missing, but I tested it all and none of this seems to be the case.

So what does this missmatch mean?

What is the parent split in this case? What does it mean, that the split has has more DOF's then it is supposed to have, and isn't the parent split just miscalculating? And how to resolve this issue.

Notes: I don't think it should be necessary, but this is the self written MarangoniForce Kernel

#include "MarangoniForce.h"

registerMooseObject("otterApp", MarangoniForce);

InputParameters MarangoniForce::validParams() {
    InputParameters params = VectorKernel::validParams();
    params.addClassDescription("Marangoni surface force term");

    params.addRequiredCoupledVar("concentration", "The Gradient of this is the driving force for surface flows");
    
    params.addParam<Real>("marangoniCoefficient", 1.0, "Absolute Value of Derivative of surface tension with respect to concentration");

    return params;
}

MarangoniForce::MarangoniForce(const InputParameters & parameters) :
    VectorKernel(parameters),
    _marangoniCoefficient(getParam<Real>("marangoniCoefficient")),

    _grad_concentration(coupledGradient("concentration")),
    _concentration_id(coupled("concentration")),
    _concentration_var(*getVar("concentration", 0)),
    _standard_grad_phi(_assembly.gradPhi(_concentration_var))

{}

Real MarangoniForce::computeQpResidual() {
    return -_marangoniCoefficient * _grad_concentration[_qp] * _test[_i][_qp];
}


Real MarangoniForce::computeQpJacobian() {
    return 0.0;
}

Real MarangoniForce::computeQpOffDiagJacobian(unsigned jvar) {
    // might need to re-add mask, look into git! I have saved a version where it still has a mask

    if (jvar == _concentration_id) {
        return - _marangoniCoefficient * _standard_grad_phi[_j][_qp] * _test[_i][_qp];
    } else return 0.0;
}

[lindsayad]

Try adding the vars line in your up block

[Preconditioning]
  [./FSP]
    type = FSP
    topsplit = 'upc'
    [./upc]
      splitting = 'up c'
      splitting_type = additive
    [../]
    [./up]
      vars = 'velocity pressure'
      splitting = 'u p' # 'u' and 'p' are the names of subsolvers
      splitting_type  = schur
      petsc_options_iname = '-pc_fieldsplit_schur_fact_type  -pc_fieldsplit_schur_precondition'
      petsc_options_value = 'full selfp'
    [../]
    [./u]
      vars = 'velocity'
      petsc_options_iname = '-pc_type -ksp_type -ksp_rtol'
      petsc_options_value = '     hypre gmres 1e-4'
    [../]
    [./p]
      vars = 'pressure'
      petsc_options_iname = '-pc_type -ksp_type -ksp_rtol'
      petsc_options_value = '   jacobi    gmres     1e-4'
    [../]
    [./c]
      vars = 'concentration'
      petsc_options_iname = '-pc_type -ksp_type -ksp_rtol'
      petsc_options_value = '   jacobi    gmres     1e-4'
    [../]
  [../]
[]