SolutionIC with specific subdomains

[rnpratoori]

I have a solution with 2 subdomains that I want to read into another simulation. The new simulation has 3 subdomains - 2 remaining exactly the same as the previous solution. To initialize this, I use the below Mesh, Variable and IC block -

[GlobalParams]
    block = '0 2'
[]

[Mesh]
    add_subdomain_ids = '1'
    [2d]
        # generate a 2D mesh
        type = GeneratedMeshGenerator
        dim = 2
        nx = ${nx}
        ny = ${ny}
        xmax = ${dx}
        ymax = ${dy}
    []
    # Subdomain for ramp
    [c3_domain]
        type = ParsedSubdomainMeshGenerator
        block_id = 2
        combinatorial_geometry = 'y > 1'
        input = 2d
    []
[]

[MeshModifiers]
    [void]
        type = CoupledVarThresholdElementSubdomainModifier
        coupled_var = eta
        criterion_type = BELOW
        subdomain_id = 0
        complement_subdomain_id = 1
        threshold = 1
        execute_on = 'INITIAL'
    []
[]

[Variables]
    # polymer volume fraction
    [c1]
        order = FIRST
        family = LAGRANGE
    []
    # Chemical potential (nJ/mol)
    [w1]
        order = FIRST
        family = LAGRANGE
    []
    # polymer volume fraction
    [c2]
        order = FIRST
        family = LAGRANGE
    []
    # Chemical potential (nJ/mol)
    [w2]
        order = FIRST
        family = LAGRANGE
    []
    # AC variable
    [eta]
        order = FIRST
        family = LAGRANGE
        block = '0  1   2'
    []
[]

[ICs]
    [c1]
        type = SolutionIC
        from_variable = 'c'
        solution_uo = 2phase_polymer
        variable = c1
        block = 0
        from_subdomains = 0
    []
    [c2]
        type = SolutionIC
        from_variable = 'c2'
        solution_uo = 2phase_polymer
        variable = c2
        block = 0
        from_subdomains = 0
    []
    [etaIC]
        type = SolutionIC
        from_variable = 'eta'
        solution_uo = 2phase_void
        variable = eta
        block = '0  1'
        from_subdomains = '0    1'
    []
    [top_c1]
        type = ConstantIC
        value = ${delta}
        variable = c1
        block = 2
    []
    [top_c2]
        type = ConstantIC
        value = ${delta}
        variable = c2
        block = 2
    []
    [top_eta]
        type = ConstantIC
        value = ${delta}
        variable = eta
        block = 2
    []
[]

[UserObjects]
    [2phase_polymer]
        type = SolutionUserObject
        mesh = 'output/2p_void.e'
        system_variables = 'c   c2'
        timestep = LATEST
        # block = '0'
    []
    [2phase_void]
        type = SolutionUserObject
        mesh = 'output/2p_void.e'
        system_variables = 'eta'
        timestep = LATEST
        # block = '0  1'
    []
[]

This should read the c variable only from subdomain 0 of the solution to the subdomain 0 of the new simulation. But it is accessing points in subdomain 1 and throwing an error that variable c does not exist. Could it be because the subdomain 1 in the new simulation is not yet initialized? Or is it something else?

[GiudGiud]

hello

it's likely coming from the generic projector.
Use this branch for now please
#30826
It should have the fix you need

It's not merged because it's not as general as it needs to be. But for your variable type it works

Guillaume

[GiudGiud]

Use an IC for eta and an auxkernel executed on INITIAL for c1, and force_preic = false for the SubdomainModifier

This should obey the sequence that you desire

[rnpratoori]

I did this -

[GlobalParams]
    block = '0 2'
[]

[Mesh]
    add_subdomain_ids = '1'
    [2d]
        # generate a 2D mesh
        type = GeneratedMeshGenerator
        dim = 2
        nx = ${nx}
        ny = ${ny}
        xmax = ${dx}
        ymax = ${dy}
        show_info = true
    []
    # Subdomain for ramp
    [c3_domain]
        type = ParsedSubdomainMeshGenerator
        block_id = 2
        combinatorial_geometry = 'y > 1'
        input = 2d
        show_info = true
    []
[]

[MeshModifiers]
    [void]
        type = CoupledVarThresholdElementSubdomainModifier
        coupled_var = eta
        criterion_type = BELOW
        subdomain_id = 0
        complement_subdomain_id = 1
        threshold = 1
        execute_on = 'INITIAL'
        force_preic = false
    []
[]

[Variables]
    # polymer volume fraction
    [c1]
        order = FIRST
        family = LAGRANGE
    []
    # Chemical potential (nJ/mol)
    [w1]
        order = FIRST
        family = LAGRANGE
    []
    # polymer volume fraction
    [c2]
        order = FIRST
        family = LAGRANGE
    []
    # Chemical potential (nJ/mol)
    [w2]
        order = FIRST
        family = LAGRANGE
    []
    # AC variable
    [eta]
        order = FIRST
        family = LAGRANGE
        block = '0  1   2'
    []
[]

[ICs]
    [eta]
        type = SolutionIC
        from_variable = 'eta'
        solution_uo = 2phase_void
        variable = eta
        block = '0'
        from_subdomains = '0    1'
    []
    [top_c1]
        type = ConstantIC
        value = ${delta}
        variable = c1
        block = 2
    []
    [top_c2]
        type = ConstantIC
        value = ${delta}
        variable = c2
        block = 2
    []
    [top_eta]
        type = ConstantIC
        value = ${delta}
        variable = eta
        block = 2
    []
[]

[UserObjects]
    [2phase_polymer]
        type = SolutionUserObject
        mesh = 'output/2p_void.e'
        system_variables = 'c   c2'
        timestep = LATEST
        # block = '0'
    []
    [2phase_void]
        type = SolutionUserObject
        mesh = 'output/2p_void.e'
        system_variables = 'eta'
        timestep = LATEST
        # block = '0  1'
    []
[]

[AuxVariables]
    [f_density]
        order = CONSTANT
        family = MONOMIAL
    []
    [f_int_density]
        order = CONSTANT
        family = MONOMIAL
    []
[]

[Kernels]
    [w1_dot]
        type = CoupledTimeDerivative
        variable = w1
        v = c1
    []
    [coupled_res1]
        type = SplitCHWRes
        variable = w1
        mob_name = M1
    []
    [coupled_parsed1]
        type = SplitCHParsed
        variable = c1
        coupled_variables = 'c2'
        f_name = f_mix
        kappa_name = kappa
        w = w1
    []
    [w2_dot]
        type = CoupledTimeDerivative
        variable = w2
        v = c2
    []
    [coupled_res2]
        type = SplitCHWRes
        variable = w2
        mob_name = M2
    []
    [coupled_parsed2]
        type = SplitCHParsed
        variable = c2
        coupled_variables = 'c1'
        f_name = f_mix
        kappa_name = kappa
        w = w2
    []
    [null]
        type = NullKernel
        variable = eta
        block = '0  1 2'
    []
[]

This gives the error - Variable 'c1' does not exist in this system. Is this because the AuxKernel is looking for c1 in AuxVariables? So, should I change c1 and c2 to AuxVariables?

Also, using force_preic (be it true or false), doesn't seem to do anything. void is still being executed before eta.

[DBG][ACT] TASK (         add_user_object) TYPE (             AddUserObjectAction) NAME (  2phase_polymer) Memory usage 190MB
[DBG][ACT] TASK (         add_user_object) TYPE (             AddUserObjectAction) NAME (     2phase_void) Memory usage 190MB
[DBG][ACT] TASK (       add_mesh_modifier) TYPE (          AddMeshModifiersAction) NAME (            void) Memory usage 190MB
[DBG][ACT] TASK (     load_surrogate_data) TYPE (         LoadSurrogateDataAction) NAME (                ) Memory usage 190MB
[DBG][ACT] TASK (    load_covariance_data) TYPE (        LoadCovarianceDataAction) NAME (                ) Memory usage 190MB
[DBG][ACT] TASK (                  add_ic) TYPE (       AddInitialConditionAction) NAME (             eta) Memory usage 190MB
[DBG][ACT] TASK (                  add_ic) TYPE (       AddInitialConditionAction) NAME (          top_c1) Memory usage 190MB
[DBG][ACT] TASK (                  add_ic) TYPE (       AddInitialConditionAction) NAME (          top_c2) Memory usage 190MB
[DBG][ACT] TASK (                  add_ic) TYPE (       AddInitialConditionAction) NAME (         top_eta) Memory usage 190MB

[GiudGiud]

oh right yes. if c1 is a nonlinear variable you cannot use auxkernels.
You could write a custom 'Predictor' object

[rnpratoori]

I am trying to do this, and I was wondering if there is a way to separate variables being worked on with the Predictor object. The base class only implements complete transfer of the NumericVector, but that may cause problems as I am only interested some of the variables and only parts of the mesh.

Also, is there a way I can merge the blocks of a solution in a postprocessor or some way. Although it is more of a duct-taping way, I feel, it would solve the major problem of transferring the solution to the new mesh. Then dividing it into different domains could be easier.

[GiudGiud]

I am trying to do this, and I was wondering if there is a way to separate variables being worked on with the Predictor object. The base class only implements complete transfer of the NumericVector, but that may cause problems as I am only interested some of the variables and only parts of the mesh.

Yes your need is more complicated than basic operations on the solution vector. But it's doable, you just use the relevant APIs on the system to implement variable-restriction in the operations on the solution vector.
There's an example here doing block and variable-restriction while working on the residual vector.
https://github.com/idaholab/moose/pull/31184/files#diff-68a2a7da0b0bb59f1141555b01113186418e23eec75c5935f68b884fdd2ae8fa

Also, is there a way I can merge the blocks of a solution in a postprocessor or some way. Although it is more of a duct-taping way, I feel, it would solve the major problem of transferring the solution to the new mesh. Then dividing it into different domains could be easier.

This does not make sense in moose-speak. What do you mean by "merging" blocks of a solution?
A postprocessor can only output a single scalar. Do you mean a predictor?

[rnpratoori]

I was trying to do this and it got quickly complicated for me to handle it. But I got a new idea to probably circumvent this.

A summary of what I have -
Old simulation: (1x1 domain)

• c1 and c2 are restricted on block 0
• eta (LatticeSmoothCircle) is used to determine block 1

Since I am facing issues in creating block 1 in the new simulation (based on eta), I added 2 additional aux variables (c1_total and c2_total) which are present across the domain. This should remove the necessity to create block 1 to be able to read the variables correctly.

Now, ideally I should be able to read the variables correctly, then divide the subdomains based on eta. But, now I get an error that c1 is trying to access elements in domain 2. But this miraculously works when I comment out the creation of block 1 (CoupledVarThresholdElementSubdomainModifier), albeit the wrong physics. Is there anything I am missing in my approach?

Here is the code I have -

[Mesh]
    add_subdomain_ids = '1'
    [2d]
        # generate a 2D mesh
        type = GeneratedMeshGenerator
        dim = 2
        nx = ${nx}
        ny = ${ny}
        xmax = ${dx}
        ymax = ${dy}
        # show_info = true
    []
    # Subdomain for ramp
    [c3_domain]
        type = ParsedSubdomainMeshGenerator
        block_id = 2
        combinatorial_geometry = 'y > 1'
        input = 2d
        # show_info = true
    []
[]

[MeshModifiers]
    [void]
        type = CoupledVarThresholdElementSubdomainModifier
        coupled_var = eta
        criterion_type = BELOW
        subdomain_id = 0
        complement_subdomain_id = 1
        threshold = 1
        execute_on = 'INITIAL'
    []
[]

[Variables]
    # polymer volume fraction
    [c1]
        order = FIRST
        family = LAGRANGE
    []
    # Chemical potential (nJ/mol)
    [w1]
        order = FIRST
        family = LAGRANGE
    []
    # polymer volume fraction
    [c2]
        order = FIRST
        family = LAGRANGE
    []
    # Chemical potential (nJ/mol)
    [w2]
        order = FIRST
        family = LAGRANGE
    []
    # AC variable
    [eta]
        order = FIRST
        family = LAGRANGE
        block = '0  1   2'
    []
[]

[ICs]
    [c1]
        type = SolutionIC
        from_variable = 'c1_total'
        solution_uo = 2phase
        variable = c1
        from_subdomains = '0 1'
        block = '0'
    []
    [c2]
        type = SolutionIC
        from_variable = 'c2_total'
        solution_uo = 2phase
        variable = c2
        from_subdomains = '0 1'
        block = '0'
    []
    [eta]
        type = SolutionIC
        from_variable = 'eta'
        solution_uo = 2phase
        variable = eta
        from_subdomains = '0 1'
        block = '0'
    []
    [top_c1]
        type = ConstantIC
        value = ${delta}
        variable = c1
        block = 2
    []
    [top_c2]
        type = ConstantIC
        value = ${delta}
        variable = c2
        block = 2
    []
    [top_eta]
        type = ConstantIC
        value = ${delta}
        variable = eta
        block = 2
    []
[]

[UserObjects]
    [2phase]
        type = SolutionUserObject
        mesh = 'output/2p_void.e'
        system_variables = 'c1_total c2_total eta'
        timestep = LATEST
    []
[]

[rnpratoori]

The issue was in CoupledVarThresholdElementSubdomainModifier. I was using the wrong subdomain_id