Maintaining Variable Dependencies, Using ADMatDiffusion

[swonner]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

I am currently writing an input file to solve to following system of equations. The input file for this system is attached below.

Many of the variables are connected, such as the material property terms L which are functions of XV and XB, and the auxiliary variable terms $\mu$ terms which are functions of X as well. I created a diagram that shows how my code currently handles this system of equations.

I am mostly concerned about the implementation of my ADMatDiffusion Kernels. My concern is that the material properties and auxillary variables being used are functions of my main variables (XA, XB, XV and XI), and I am unsure if right now that dependency is taken into account. How might I ensure that this dependency is properly taken into account? Is it okay that I split the system up in this way?

Input File

[Mesh] 
    [gen]
        type = GeneratedMeshGenerator
        dim = 2
        nx = 200
        ny = 50
        xmax = 50 # um
        ymax = 10 #um
    []
    [left_box]
        type = SubdomainBoundingBoxGenerator
        input = gen
        bottom_left = '0 0 0' 
        top_right = '25 10 0' # (xmax/2, ymax, 0)
        block_id = 1
    []
    [grainboundary]
        type = SideSetsBetweenSubdomainsGenerator
        input = left_box
        paired_block = 1
        primary_block = 0
        new_boundary = grain_boundary
    [] 
[]




[Variables] 
    # XA, XB, XV, XI - the concentration of atoms (A,B) and defects (V,I)
    [XA]
        order = FIRST
        family = LAGRANGE
    []

    [XB]
        order = FIRST
        family = LAGRANGE
    []

    [XV]
        order = FIRST
        family = LAGRANGE
    []

    [XI]
        order = FIRST
        family = LAGRANGE
    []  
[]

[AuxVariables] 
    # Chemical Potential of A, B atoms and V, I defects
    [muA]
        order = FIRST 
        family = MONOMIAL
    []

    [muB]
        order = FIRST 
        family = MONOMIAL
    []

    [muV]
        order = FIRST 
        family = MONOMIAL
    []

    [muI]
        order = FIRST 
        family = MONOMIAL
    []


[]


[ICs]
    [XA_ic]
        type = ConstantIC
        variable = XA
        value = 0.95
    []

    [XB_ic]
        type = ConstantIC
        variable = XB
        value = 0.05 
    []  

    [XV_ic]
        type = ConstantIC
        variable = XV
        value = 9.73e-15  
    []

    [XI_ic]
        type = ConstantIC
        variable = XI
        value = 595e-33
    []
    [muA_ic]
        type = ConstantIC
        variable = muA
        value = -2.21e-03   
    []
    [muB_ic]
        type = ConstantIC
        variable = muB
        value = -1.29e-01 
    []
    [muI_ic]
        type = ConstantIC
        variable = muI
        value = 0 
    []
    [muV_ic]
        type = ConstantIC
        variable = muV
        value = 0 
    []


[]

[Kernels]
    #################### Main variable: XA ####################
    ## 
    # XA + XB - XV + XI - 1 = 0
    [conservation_XA] # XA
        type =ADMatReaction
        variable = XA
        reaction_rate = 1
    []

    [conservation_XB] # XB
        type = ADCoupledForce
        variable = XA
        v = XB
    []

    [conservation_XV] # -XV
        type = ADCoupledForce
        variable = XA
        v = XV
        coef = -1
    []  

    [conservation_XI] # XI
        type = ADCoupledForce
        variable = XA
        v = XI
    []

    [conservation_one] # -1
        type = ADBodyForce
        variable = XA
        value = -1
    []

    
    #################### Main Variable: XB ####################
    ## (𝜕𝑋𝐵)/𝜕𝑡 + Vat*∇∙𝐽𝐵𝑥 = 0 
    [timeXB] # (𝜕𝑋𝐵)/𝜕𝑡
        type = ADTimeDerivative
        variable = XB
    []

    # Vat*∇∙𝐽𝐵 = Vat*∇∙1/(𝑘𝐵*𝑇)*[𝐿𝐵𝐵V*(∇𝜇𝐵−∇𝜇𝑉) + 𝐿𝐴𝐵𝑉*(∇𝜇𝐴−∇𝜇𝑉)]
    [gradfluxB_1]
    # Vat*∇∙1/(𝑘𝐵*𝑇)*[𝐿𝐵𝐵V*(∇𝜇𝐵)]
        type = ADMatDiffusion
        variable = XB
        v = muA  # coupled variable
        diffusivity = LBBV_kbT #Vat*LBBV/kb*T
    []   

    [gradfluxB_2]
    # Vat*∇∙1/(𝑘𝐵*𝑇)*[𝐿𝐵𝐵V*(-∇𝜇V)]
        type = ADMatDiffusion
        variable = XB
        v = muV # coupled variable
        diffusivity = LBBV_kbT_neg #-Vat*LBBV/kb*T
    []  

    [gradfluxB_3]
    # Vat*∇∙1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐵𝑉*(∇𝜇𝐴)]
        type = ADMatDiffusion
        variable = XB
        v = muA # coupled variable
        diffusivity = LABV_kbT #Vat*LABV/kb*T
    []  

    [gradfluxB_4]
    # Vat*∇∙1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐵𝑉*(-∇𝜇V)]
        type = ADMatDiffusion
        variable = XB
        v = muV # coupled variable
        diffusivity = LABV_kbT_neg #-Vat*LABV/kb*T
    []  



    #################### Main Variable: XV ####################
    ## (𝜕𝑋𝑉)/𝜕𝑡 + Vat*∇∙(𝐽𝑉x) - 𝐾0 + 𝐾𝐼𝑉*𝑋𝑉*𝑋𝐼 = 0
    [timeXV] #(𝜕𝑋𝑉)/𝜕𝑡
        type = ADTimeDerivative
        variable = XV
    []

   #Vat* ∇∙𝐽𝑉 = Vat*∇∙(−1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐴V*(∇𝜇𝐴−∇𝜇𝑉) + 𝐿𝐴𝐵𝑉*(∇𝜇𝐵−∇𝜇𝑉) + 𝐿𝐴𝐵𝑉*(∇𝜇𝐴−∇𝜇𝑉) + 𝐿𝐵𝐵V*(∇𝜇𝐵−∇𝜇𝑉))
    [gradfluxV_1]
    # Vat*∇∙(−1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐴V*(∇𝜇𝐴)]
        type = ADMatDiffusion
        variable = XV
        v = muA # coupled variable
        diffusivity = LAAV_kbT_neg #-Vat*LAAV/kb*T
    []  

    [gradfluxV_2]
    # Vat*∇∙(−1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐴V*(-∇𝜇V)]
        type = ADMatDiffusion
        variable = XV
        v = muV # coupled variable
        diffusivity = LAAV_kbT #Vat*LAAV/kb*T
    []  

    [gradfluxV_3]
    # Vat*∇∙(−1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐵𝑉*(∇𝜇𝐵)]
        type = ADMatDiffusion
        variable = XV
        v = muB # coupled variable
        diffusivity = LABV_kbT_neg #-Vat*LABV/kb*T
    []  

    [gradfluxV_4]
    # Vat*∇∙(−1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐵𝑉*(−∇𝜇𝑉)]
        type = ADMatDiffusion
        variable = XV
        v = muV # coupled variable
        diffusivity = LABV_kbT #Vat*LABV/kb*T
    []  

    [gradfluxV_5]
    #Vat*∇∙(−1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐵𝑉*(∇𝜇𝐴)]
        type = ADMatDiffusion
        variable = XV
        v = muA # coupled variable
        diffusivity = LABV_kbT_neg #-Vat*LABV/kb*T
    [] 

    [gradfluxV_6]
    # Vat*∇∙(−1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐵𝑉*(−∇𝜇𝑉)]
        type = ADMatDiffusion
        variable = XV
        v = muV # coupled variable
        diffusivity = LABV_kbT #Vat*LABV/kb*T
    [] 
 
    [gradfluxV_7]
    # Vat*∇∙(−1/(𝑘𝐵*𝑇)*[𝐿𝐵𝐵V*(∇𝜇𝐵)]
        type = ADMatDiffusion
        variable = XV
        v = muB # coupled variable
        diffusivity = LBBV_kbT_neg #-Vat*LBBV/kb*T
    []
    
    [gradfluxV_8]
    # Vat*∇∙(−1/(𝑘𝐵*𝑇)*[𝐿𝐵𝐵V*(−∇𝜇𝑉)]
        type = ADMatDiffusion
        variable = XV
        v = muV # coupled variable
        diffusivity = LBBV_kbT #Vat*LBBV/kb*T
    []

    [defectgenV] #- 𝐾0
        type = ADMatBodyForce
        variable = XV
        material_property = K0 
    []

    [defectgrecombinationV] #𝐾𝐼𝑉*𝑋𝑉*𝑋𝐼*Vat
        type = ADMatReaction
        variable = XV
        reaction_rate = RR_I #𝐾𝐼𝑉*𝑋𝐼
    []

    #################### Main Variable: XI ####################
    ## (𝜕𝑋𝐼)/𝜕𝑡 + Vat*∇∙(𝐽𝐼x) - 𝐾0 + 𝐾𝐼𝑉*𝑋𝑉*𝑋𝐼 = 0
    [timeXI] #(𝜕𝑋𝐼)/𝜕𝑡
        type = ADTimeDerivative
        variable = XI
    []

    #Vat*∇∙𝐽𝐼 = Vat*∇∙(1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐴𝐼*(∇𝜇𝐴+∇𝜇𝐼)])
    [gradfluxI_1]
    # Vat*∇∙(1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐴𝐼*(∇𝜇𝐴)]
        type = ADMatDiffusion
        variable = XI
        v = muA # coupled variable
        diffusivity = LAAI_kbT #Vat*LAAI/kb*T
    []
    [gradfluxI_2]
    # Vat*∇∙(1/(𝑘𝐵*𝑇)*[𝐿𝐴𝐴𝐼*(∇𝜇I)]
        type = ADMatDiffusion
        variable = XI
        v = muI # coupled variable
        diffusivity = LAAI_kbT #Vat*LAAI/kb*T
    []

    [defectgenI] #- 𝐾0
        type = ADMatBodyForce
        variable = XV
        material_property = K0 # K0 
    []
    [defectgrecombinationI] #𝐾𝐼𝑉*𝑋𝑉*𝑋𝐼  
        type = ADMatReaction
        variable = XI
        reaction_rate = RR_V
    []
[]



[AuxKernels]
    # Guidance using older (outdated) ParsedAux, https://mooseframework.inl.gov/source/auxkernels/ParsedAux.html
    # refer to modules/phase_field/examples/anisotropic_interfaces/echebarria_iso.i for old syntax
    # new term == old term
    # coupled_variables = args 
    # material_properties == I do not think the older version included a way to import material properties
        # instead you should use constant_names and constant_expressions 
    # expression == function 

    
    [chempotA]
        ## 𝜇A = 𝐺𝐴 + 𝑘B*𝑇*𝑙𝑛(𝑋𝐴) ##
        type = ParsedAux 
        constant_names = 'kb T G_A'
        constant_expressions = '8.61733326e-5 500 0'
        args = 'XA' # this is correct
        variable = muA
        function = 'G_A + kb*T*log(XA)' 

        # written for the newer ParsedAux
        # type = parsedAux
        # material_properties = 'kb T G_A'
        # coupled_variables = 'XA'
        # variable = muA
        # expression = 'G_A + kb*T*log(XA)'
    []

    [chempotB] 
        ## 𝜇B = 𝐺B + 𝑘B*𝑇*𝑙𝑛(𝑋B) ##
        type = ParsedAux 
        constant_names = 'kb T G_B'
        constant_expressions = '8.61733326e-5 500 0'
        args = 'XB'
        variable = muB
        function = 'G_B + kb*T*log(XB)' 

        # written for the newer ParsedAux
        # type = parsedAux
        # material_properties = 'kb T G_B'
        # coupled_variables = 'XB'
        # variable = muB
        # expression = 'G_B + kb*T*log(XB)' 
    []
    
            
    [chempotV] 
        ## 𝜇V =  𝑘B*𝑇*𝑙𝑛(𝑋V/𝑋Veq) ##
        type = ParsedAux 
        constant_names = 'kb T XV_eq'
        constant_expressions = '8.61733326e-5 500 1'
        args = 'XV'
        variable = muV
        function = 'kb*T*log(XV/XV_eq)'

        # written for the newer ParsedAux
        # type = parsedAux
        # material_properties = 'kb T XV_eq'
        # coupled_variables = 'XV'
        # variable = muV
        # expression = 'kb*T*log(XV/XV_eq)' 
        
    []

    [chempotI] 
        ## 𝜇I =  𝑘B*𝑇*𝑙𝑛(𝑋I/𝑋Ieq) ##
        type = ParsedAux 
        constant_names = 'kb T XI_eq'
        constant_expressions = '8.61733326e-5 500 1'
        args = 'XI'
        variable = muI
        function = 'kb*T*log(XI/XI_eq)'

        # written for the newer ParsedAux
        # type = parsedAux
        # material_properties = 'kb T XI_eq'
        # coupled_variables = 'XI'
        # variable = muI
        # expression = 'kb*T*log(XI/XI_eq)' 
    []
    
[]

[Materials]
    # https://mooseframeworkinlgov/source/materials/GenericConstantMaterialhtml
    [SettingGibbs]
        type = ADGenericConstantMaterial 
        prop_names = 'G_A G_B'
        prop_values = '0 0' 
    []

    [Settingboltsmannandtemp]
        type = ADGenericConstantMaterial
        # kb has units of ev/K
        # T has units of K
        prop_names = 'kb T'
        prop_values ='8.61733326E-5 500' 
    []
    
    [SettingKineCluEConstants]
        type = ADGenericConstantMaterial
        prop_names = 'LVV_v LVV_vb LVB_vb LBB_vb'
        prop_values = '1 1 1 1' #ADD VALUES LATER
    []
    [Settingparitionfunctions]
        type = ADGenericConstantMaterial
        prop_names = 'z_v z_vb'
        prop_values = '1 1' #ADD VALUES LATER
    []
    [SettingEQCon]
        type = ADGenericConstantMaterial
        prop_names = 'XI_eq XV_eq'
        prop_values = '1 1' #ADD VALUES LATER
    []

    [SettingDI]
        type = ADGenericConstantMaterial
        prop_names = 'DI'
        prop_values = '1' #ADD VALUES LATER
    []  
    [SettingVat]
        type = ADGenericConstantMaterial
        prop_names = 'Vat'
        prop_values = '1' #ADD VALUES LATER
    []  
    [SettingIRRates]
        type = ADGenericConstantMaterial
        prop_names = 'K0 KIV'
        prop_values = '1 1' #ADD VALUES LATER
    []  

    # Guidance using older (outdated) ADDerivativeParsedMaterial
    # https://mooseframework.inl.gov/source/materials/DerivativeParsedMaterial.html
    # new term == old term
    # coupled_variables = args see example modules/phase_field/test/tests/GrandPotentialPFM/GrandPotentialMultiphase_AD.i
    # property_name == f_name 
    # expression == function 

    [clustercon_B] 
        # 0 = Ω^3*zvb*CB^2 + (Ω^2*zv + XV*Ω^2*zvb - XB*Ω^2*zvb)*CB - Ω*zv*XB
        # a = Ω^3*zvb = Vat*Vat*Vat*z_vb
        # b = (Ω^2*zv + XV*Ω^2*zvb - XB*Ω^2*zvb) = (Vat*Vat*z_v + XV*Vat*Vat*z_vb - XB*Vat*Vat*z_vb)
        # c = - Ω*zv = -1*Vat*z_v
        # CB = (-b + (b^(2) - 4*a*c)^(1/2)) / (2*a)
        # CB = '(-(Vat*Vat*z_v + XV*Vat*Vat*z_vb - XB*Vat*Vat*z_vb) + ((Vat*Vat*z_v + XV*Vat*Vat*z_vb - XB*Vat*Vat*z_vb)^(2) - 4*Vat*Vat*Vat*z_vb*-1*Vat*z_v)^(1/2)) / (2*Vat*Vat*Vat*z_vb)' 
        type = ADDerivativeParsedMaterial
        f_name = CB
        args = 'XV XB'
        material_property_names = 'z_v z_vb Vat'
        function = '(-(Vat*Vat*z_v + XV*Vat*Vat*z_vb - XB*Vat*Vat*z_vb) + ((Vat*Vat*z_v + XV*Vat*Vat*z_vb - XB*Vat*Vat*z_vb)^(2) - 4*Vat*Vat*Vat*z_vb*-1*Vat*z_v)^(1/2)) / (2*Vat*Vat*Vat*z_vb)' 
        
        # written for the new ADDerivativeParsedMaterial
        # property_name = CB
        # coupled_variables = 'XV XB'
        # material_property_names = 'z_v z_vb Vat'
        # expression = '(-(Vat*Vat*z_v + XV*Vat*Vat*z_vb - XB*Vat*Vat*z_vb) + ((Vat*Vat*z_v + XV*Vat*Vat*z_vb - XB*Vat*Vat*z_vb)^(2) - 4*Vat*Vat*Vat*z_vb*-1*Vat*z_v)^(1/2)) / (2*Vat*Vat*Vat*z_vb)' 
    []


    [clustercon_V] # CV = XV / (Ω*zv + Ω^2*CB*zvb)
    type = ADDerivativeParsedMaterial 
        f_name = CV
        material_property_names = 'z_v z_vb Vat CB(XV,XB)'
        args = 'XV'
        function = 'XV / (Vat*z_v + Vat*Vat*CB*z_vb)' 

        # written for the new ADDerivativeParsedMaterial
        # property_name = CV
        # coupled_variables = 'XV'
        # material_property_names =  'z_v z_vb Vat CB'
        # expression = 'XV / (z_v + Vat*CB*z_vb)' 
    []


    [SolvingLVV] ## 𝐿𝑉𝑉 = 𝐶𝑉*𝑧v*𝐿𝑉𝑉𝑒𝑞(𝑣) + 𝐶𝑉*𝐶𝐵*𝑧𝑣𝑏*𝐿VVeq(𝑣𝐵)*Ω ##
        type = ADDerivativeParsedMaterial
        f_name = LVV
        material_property_names = 'z_v LVV_v LVV_vb CV(XV,XB) CB(XV,XB) Vat' 
        function = 'CV*z_v*LVV_v + CV*CB*LVV_vb*Vat'

        
        # written for the new ADDerivativeParsedMaterial
        # property_name = LVV
        # material_property_names =  'z_v LVV_v LVV_vb CV CB'
        # expression = 'CV*z_v*LVV_v + CV*CB*LVV_vb'
    []


    [SolvingLBBV] ##  𝐿𝐵𝐵𝑉 = 𝐶𝑉*𝐶𝐵*𝑧𝑣𝑏*𝐿𝐵𝐵(𝑣,𝑒𝑞)(𝑣𝐵)*Ω ##
        type =  ADDerivativeParsedMaterial
        f_name = LBBV
        material_property_names = 'z_vb LBB_vb CV(XV,XB) CB(XV,XB) Vat' 
        function =  'CV*CB*z_vb*LBB_vb*Vat'

                
        # written for the new ADDerivativeParsedMaterial
        # property_name = LBBV
        # material_property_names = 'z_vb LBB_vb CV CB'
        # expression =  'CV*CB*z_vb*LBB_vb'
    []        

    [SolvingLBV] ## 𝐿𝐵𝑉 = 𝐶𝑉*𝐶𝐵*𝑧𝑣𝑏*𝐿𝐵𝑉𝑒𝑞(vB)*Ω ##
        type =  ADDerivativeParsedMaterial
        f_name = LBV
        material_property_names = 'z_vb LVB_vb CV(XV,XB) CB(XV,XB) Vat'  #  CV(XV, XB) CB(XV, XB)?
        function =  'CV*CB*z_vb*LVB_vb*Vat'

        # written for the new ADDerivativeParsedMaterial
        # property_name = LBV
        # material_property_names = 'z_vb LVB_vb CV CB'
        # expression =  'CV*CB*z_vb*LVB_vb'

    []

    [SolvingLABV] ## 𝐿𝐴𝐵𝑉 = −(𝐿B𝑉 + 𝐿𝐵𝐵𝑉) ##
        type =  ADDerivativeParsedMaterial
        f_name = LABV
        material_property_names = 'LBV(XV,XB) LBBV(XV,XB)'  
        function = '-1*(LBV + LBBV)'

        # written for the new ADDerivativeParsedMaterial
        # property_name = LABV
        # material_property_names = 'LBV LBBV'  
        # expression = '-1*(LBV + LBBV)'
    []
    
    [SolvingLAAV] ## 𝐿𝐴𝐴𝑉 = 𝐿𝑉𝑉 − 𝐿𝐵𝐵𝑉 + 2*(𝐿𝑉𝐵 + 𝐿𝐵𝐵𝑉)
        type = ADDerivativeParsedMaterial
        f_name = LAAV
        material_property_names = 'LVV(XV,XB) LBBV(XV,XB)  LBV(XV,XB)' 
        function = 'LVV - LBBV + 2*(LBV + LBBV)'

        # written for the new ADDerivativeParsedMaterial
        # property_name = LAAV
        # material_property_names = 'LVV LBBV  LBV'   
        # expression = 'LVV - LBBV + 2*(LBV + LBBV)'

    []
    
    [SolvingLAAI] ##𝐿𝐴𝐴𝐼 = X𝐼*𝐷𝐼 
        type = ADDerivativeParsedMaterial  
        f_name = LAAI
        material_property_names = 'DI' 
        args = 'XI'
        function = 'XI*DI'

        # written for the new ADDerivativeParsedMaterial
        # property_name = LAAI
        # material_property_names = 'DI'  
        # expression = 'XI*DI'
    []



    #### Calculating Material Properties that are directly input into the kernels block
    [laav_dividekbt]  
        type = ADDerivativeParsedMaterial  
        f_name = LAAV_kbT
        material_property_names = 'LAAV(XV,XB) kb T Vat' 
        function = 'Vat*LAAV/(kb*T)'
        # written for the new ADDerivativeParsedMaterial
        # property_name = LAAV_kbT
        # material_property_names =  'LAAV kb T Vat'  
        # expression = 'Vat*LAAV/(kb*T)'
    []
    [negativelaav_dividekbt]  
        type = ADDerivativeParsedMaterial  
        f_name = LAAV_kbT_neg
        material_property_names = 'LAAV(XV,XB) kb T Vat' 
        function = '-1*Vat*LAAV/(kb*T)'
        # written for the new ADDerivativeParsedMaterial
        # property_name = LAAV_kbT_neg
        # material_property_names =  'LAAV kb T Vat'  
        # expression = '-1*Vat*LAAV/(kb*T)'
    []

    [labv_dividekbt]  
        type = ADDerivativeParsedMaterial  
        f_name = LABV_kbT
        material_property_names = 'LABV(XV,XB) kb T Vat' 
        function = 'Vat*LABV/(kb*T)'
        # written for the new ADDerivativeParsedMaterial
        # property_name = LAAV_kbT
        # material_property_names =  'LABV kb T Vat'  
        # expression = 'Vat*LABV/(kb*T)'
    []
    [negativelabv_dividekbt]  
        type = ADDerivativeParsedMaterial  
        f_name = LABV_kbT_neg
        material_property_names = 'LABV(XV,XB) kb T Vat' 
        function = '-1*Vat*LABV/(kb*T)'
        # written for the new ADDerivativeParsedMaterial
        # property_name = LABV_kbT_neg
        # material_property_names =  'LABV kb T Vat'  
        # expression = '-1*Vat*LABV/(kb*T)'
    []

    [lbbv_dividekbt]  
        type = ADDerivativeParsedMaterial  
        f_name = LBBV_kbT
        material_property_names = 'LBBV(XV,XB) kb T Vat' 
        function = 'Vat*LBBV/(kb*T)'
        # written for the new ADDerivativeParsedMaterial
        # property_name = LBAV_kbT
        # material_property_names =  'LBBV kb T Vat'  
        # expression = 'Vat*LBBV/(kb*T)'
    []
    [negativelbbv_dividekbt]  
        type = ADDerivativeParsedMaterial  
        f_name = LBBV_kbT_neg
        material_property_names = 'LBBV(XV,XB) kb T Vat' 
        function = '-1*Vat*LBBV/(kb*T)'
        # written for the new ADDerivativeParsedMaterial
        # property_name = LBBV_kbT_neg
        # material_property_names =  'LBBV kb T'  
        # expression = '-1*LBBV/(kb*T)'
    []

    [laai_dividekbt]  
        type = ADDerivativeParsedMaterial  
        f_name = LAAI_kbT
        material_property_names = 'LAAI kb T Vat' 
        function = 'Vat*LAAI/(kb*T)'
        # written for the new ADDerivativeParsedMaterial
        # property_name = LAAI_kbT
        # material_property_names =  'LAAI kb T'  
        # expression = 'LAAI/(kb*T)'
    []

    
    [Recombination_I] # XI*KIV
        type = ADDerivativeParsedMaterial 
        f_name = RR_I
        args = 'XI'
        material_property_names = 'KIV Vat'
        function = 'XI*KIV*Vat'

        
        # written for the new ADDerivativeParsedMaterial
        # property_name = RR_I
        # coupled_variables = 'XI'
        # material_property_names = 'KIV Vat'  
        # expression = 'XI*KIV*Vat'
    []

    [Recombination_V] # XV*KIV
        type = ADDerivativeParsedMaterial 
        f_name = RR_V
        args = 'XV'
        material_property_names = 'KIV Vat'
        function = 'XV*KIV*Vat'

        # written for the new ADDerivativeParsedMaterial
        # property_name = RR_V
        # coupled_variables = 'XV'
        # material_property_names = 'KIV Vat'  
        # expression = 'XV*KIV*Vat'
    []    
     
[]


[BCs]
    [GB_XV]
        type = ADDirichletBC
        variable = XV
        boundary = grain_boundary
        value = 0
    []
    [GB_XI]
        type = ADDirichletBC
        variable = XI
        boundary = grain_boundary
        value = 0
    []
    
    [left_xa]
        type =  ADNeumannBC
        variable = XA
        boundary = left
        value = 0
    []
    [left_xb]
        type = ADNeumannBC
        variable = XB
        boundary = left
        value = 0
    []
    [left_xi]
        type = ADNeumannBC
        variable = XI
        boundary = left
        value = 0
    []
    [left_xv]
        type = ADNeumannBC
        variable = XV
        boundary = left
        value = 0
    []
    [right_xa]
        type = ADNeumannBC
        variable = XA
        boundary = right
        value = 0
    []
    [right_xb]
        type = ADNeumannBC
        variable = XB
        boundary = right
        value = 0
    []
    [right_xi]
        type = ADNeumannBC
        variable = XI
        boundary = right
        value = 0
    []
    [right_xv]
        type = ADNeumannBC
        variable = XV
        boundary = right
        value = 0
    []

    [Periodic]
        [XA_y]
            variable = XA
            auto_direction = 'y'
        []
    []
    [Periodic]
        [XB_y]
            variable = XB
            auto_direction = 'y'
        []
    []
    [Periodic]
        [XI_y]
            variable = XI
            auto_direction = 'y'
        []
    []
    [Periodic]
        [XV_y]
            variable = XV
            auto_direction = 'y'
        []
    []     
[]



[Executioner]

  solve_type = NEWTON
  l_max_its = 20 
  l_tol = 1e-6 
  nl_max_its = 15 
  nl_rel_tol = 1e-8 
  nl_abs_tol = 1e-8 
  end_time = 50 

  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = ' hypre    boomeramg'

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-9  
    cutback_factor = 0.8 
    growth_factor = 1.2 
    optimal_iterations = 7
  []
[]

[Outputs]
    file_base = ristest_11252024
    [./exodus]
        type = Exodus
    [../]
    [./outfile]
        type = CSV
        time_column = true
        show = 'XA XB XI XV'
    [../]
[]

[Debug]
    show_var_residual_norms = true
[]

Also, as a note about why I am using old syntax. I am unable to install an up-to-date copy of MOOSE on the cluster that I have access to. Therefore, I am using an older (late 2021) install. Sorry for this inconvenience.

[GiudGiud]

Hello

you’ll want to use ADMatDiffusion for the kernel to take into account any variable dependency of the diffusion coefficient automatically.

[swonner]

Hi GiudGiud, thank you for your response!

To follow up on this, is it okay for ADMatDiffusion to be used in the following way:

$\nabla$ $\cdot$ D(u) $\nabla$ v(u) , where

• D is a function of the main variable being solved in the residual (u) and D is a material property defined with ADDerivativeParsedMaterial
• v is a coupled variable, a function of the main variable, and is defined as an auxiliary variable

Will this kernel be able to still take into account the variable dependencies on u in this form?

[GiudGiud]

no this wont work.
v must be a solver variable to be able to use ADMatDiffusion

You ll need a create a new kernel. v should be written down in the kernel as a function of u, otherwise derivatives (to form the jacobian) will not be propagated

[swonner]

I see, okay. Thank you! So the material property is okay? Is this because variable dependency preserved with ADDerivativeParsedMaterial?

Regarding the variable I want to couple, its form is v = $\mu$ = ln(X).When I browsed the existing source code, I could not find a kernel that uses a natural log. Could you please advise if there might be any existing AD kernel implementations that use natural log that I might refer to?

[GiudGiud]

So the material property is okay? Is this because variable dependency preserved with ADDerivativeParsedMaterial?

because AD works well with material properties, keeps track of derivatives as they are computed
It does not through auxiliary variables

You can just use std::log
https://cplusplus.com/reference/cmath/log/
We've implemented AD support for it