CoupledVarThresholdElementSubdomainModifier - moving boundary error

[czadeh]

Hi,

I made sure MOOSE and conda were fully updated before testing and posting this and I could not find this issue posted elsewhere.

I am using CoupledVarThresholdElementSubdomainModifier in my input file to modify my active domain. I start with a thin layer and add more layers in the Z direction. The function used to modify the boundary is called 'layerZ' and its value is displayed in exodus.

Input file is here:

[Problem]
  solve = false
  kernel_coverage_check = false
[]
[Mesh]
    construct_side_list_from_node_list = true
    [gen]
      type = GeneratedMeshGenerator
      dim = 3
      nx = 20
      ny = 20
      nz = 20
      xmax = 2
      ymax = 1
      zmax = 1
    []
    [powder_domain] #the initial domain, the first layer
      type = SubdomainBoundingBoxGenerator
      input = 'gen'
      block_id = 1
      bottom_left = '0 0 0'
      top_right = '2 1 0.1'
    []
    [gas_domain] #the initially inactive domain
      type = SubdomainBoundingBoxGenerator
      input = 'powder_domain'
      block_id = 2
      bottom_left = '0 0 0.1'
      top_right = '2 1 1'
    []
[]
[Variables]
  [u]
    initial_condition = 1
    family = LAGRANGE
    block = '1 2'
  []
[]
[Kernels]
  [dummy]
    type = NullKernel
    variable = u
  []
[]
[UserObjects]
  [addLayer] #subdomain modifier
    type = CoupledVarThresholdElementSubdomainModifier
    coupled_var = 'layerZ' #the z height of the active layers
    criterion_type = BELOW
    threshold = 0
    subdomain_id = 1
    # complement_subdomain_id = 2
    moving_boundary_name = topLayer #the top boundary of the top layer
    execute_on = 'TIMESTEP_BEGIN'
  []
[]
[Functions]
  [LayerZFunc] #function to define layerZ
    type = ParsedFunction
    value = 'z-0.1*floor(t)-0.1'
    # value = 'int(10*(z-0.1*floor(t)-0.1))/10'
  []
[]
[AuxVariables]
  [layerZ]
    [InitialCondition]
      type = ConstantIC
      value = 0.1
      variable = layerZ
    []
  []
[]
[AuxKernels]
  [layerZ]
    type = FunctionAux
    variable = layerZ
    function = LayerZFunc
    execute_on = 'TIMESTEP_BEGIN'
  []
[]
[Executioner]
  type = Transient
  solve_type = JFNK
  num_steps = 5
  dt = 1
[]
[Outputs]
  exodus = true
[]

The bulk of the domain seems to behave correctly, but some unusual behavior appears on the boundary.

time = 2

According to the legend it seems like a rounding error. Not a big deal, but the moving boundary labeled 'topLayer' shows some inconsistent behavior.

time = 2

At first it matches block 1 almost exactly.

But at some point it diverges from block 1 and somehow the moving boundary is comprised of separate z-planes.

time = 4

time = 4

I need to use the moving boundary as part of my problem but I do not think it will work like this. I am struggling to figure out why this error occurs.

I appreciate your help in figuring this out.

Thanks,
Cameron

[GiudGiud]

Hello

Since you have already identified the issue as being a rounding error, just change the value field in LayerZfunc to be offset by some small epsilon like 1e-8 and it will resolve this.

Rounding errors can be inconsistent between different computations at different timesteps, there s very little that can done apart from using small offsets

Guillaume

[czadeh]

Hi Guillaume,

I tried adding an epsilon starting with e-8 up to e-2. I tried it both within the argument of the floor function and outside of it, but the error persists the same way.

Would changing the value of
execute_on
in UserObjects or AuxKernels affect this?

[GiudGiud]

It can affect this if there is an issue with dependency resolution. There should not be.

you can try to shift the order of execution of the User object if you want using:
force_preaux = true
or
force_postaux = true
in the user object parameter

By default, user objects are executed after aux kernels, so what you did should work since the layer criterion is being computed by an aux kernel.

[hugary1995]

I can test to see if I can reproduce this. I've already recorded several bugs, if this is also a bug in element subdomain modifier, I'll also add this to my to-do list. But right now I don't have a clear idea of when I can get to them.

[hugary1995]

Pinging @dewenyushu here because she might have some resources to work on this. I'm not sure.

[dewenyushu]

Yes, we can potentially use part of a ldrd project to fix the issues associated with the subdomain modifier, if you may have the availability.

[czadeh]

I tried to circumvent this issue by using SideSetsAroundSubdomainGenerator when defining my mesh and blocks initially. I used it to define the boundaries around block 1. But now the problem is that

1. The top boundary between blocks 1 and 2 works correctly until the subdomain is modified by CoupledVarThresholdElementSubdomainModifier. At that point, the boundary between the subdomains disappears.

2. The sides of block 1 (everything except z = 0 and the boundary between blocks 1 and 2) stay at their initially defined geometry and do not update with increasing z.

Is there a way to make the mesh generator update the sidesets when CoupledVarThresholdElementSubdomainModifier makes a change?

Thanks,
Cameron

[GiudGiud]

Iirc we use the initial condition objects to initialize variables on elements that just changed domains.

We could create a new IC that uses the neighbor value to initialize.

There is an issue opened by @hugary1995 that summarizes both the current options for initialization and what would be desired

[keassyguang]

If I want to develop my own IC or modify the current IC, where should I start? What source code might need to be modified？and also for the "growing" side sets

[GiudGiud]

You should likely create a new IC object for this. It's not obvious and most of the existing objects are simple

I would first try to write the algorithm on paper before writing a new IC

[keassyguang]

You should likely create a new IC object for this. It's not obvious and most of the existing objects are simple

Thank you very much.

There is an issue opened by @hugary1995 that summarizes both the current options for initialization and what would be desired

I can't find the issue you mentioned, could you kindly provide the links?

[GiudGiud]

Could not find it either. @hugary1995 if you remember if you actually created this issue or if I imagined it