Notebook

ChEn-5310: Computational Continuum Transport Phenomena Spring 2021 UMass Lowell; Prof. V. F. de Almeida 13Apr21

12. Peclet Bulk-Coupled 1D with Dirichlet Boundary Conditions¶

$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}} \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}} \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}} \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}} \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}} \newcommand{\Smtrx}{\boldsymbol{\mathsf{S}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \newcommand{\avec}{\boldsymbol{\mathsf{a}}} \newcommand{\bvec}{\boldsymbol{\mathsf{b}}} \newcommand{\cvec}{\boldsymbol{\mathsf{c}}} \newcommand{\rvec}{\boldsymbol{\mathsf{r}}} \newcommand{\fvec}{\boldsymbol{\mathsf{f}}} \newcommand{\mvec}{\boldsymbol{\mathsf{m}}} \newcommand{\gvec}{\boldsymbol{\mathsf{g}}} \newcommand{\zerovec}{\boldsymbol{\mathsf{0}}} \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert} \newcommand{\transpose}[1]{{#1}^\top} \DeclareMathOperator{\rank}{rank} \newcommand{\Reals}{\mathbb{R}} \newcommand{\thetavec}{\boldsymbol{\theta}} $

Objectives ¶

Introduce bulk-coupling in the context of a Peclet problem studied in the notebook series 09/10/11; all past notebooks must be thoroughly reviewed.
Present the Galerkin weak form of the Peclet 1D problem with two unknown variables as described below (OneNote notes here).
Some initial code is provided in the course repository but no full source code is given out. A significant effort in programing is often necessary to learn the subject well. However the material in this course is helpful with this task. Hands-on work during lectures will try to fill in existing gap. The steps in this notebook are necessary for a basic understanding of the subject.
The reader is supposed to consult the MOOSE source documentation to fill in gaps in reproducing the steps below.

Plotting Functions ¶

This is an auxiliary section for holding plotting functions used later.

In [1]:

'''Plot function for FEM Solution'''

def plot_solution(df, 
                  title='No Title', 
                  u1_legend='no u1 legend',
                  u2_legend='no u2 legend',
                  u1_flux_legend='no u1 flux legend',
                  u2_flux_legend='no u2 flux legend'):
    
    import matplotlib.pyplot as plt
    %matplotlib inline
    plt.style.use('dark_background')

    (fig, ax1) = plt.subplots(1, figsize=(15, 6))

    ax1.plot(df['x'], df['u'],'r*-',label=u1_legend)
    
    if 'u2' in df.columns:
        ax1.plot(df['x'], df['u2'],'r*--',label=u2_legend)

    ax1.set_xlabel(r'$x$ [cm]', fontsize=18)
    ax1.set_ylabel(r'$u_h(x)$ [g/cc]', fontsize=18, color='red')
    ax1.tick_params(axis='y', labelcolor='red', labelsize=14)
    ax1.tick_params(axis='x', labelsize=14)
    ax1.legend(loc='center left', fontsize=12)
    #ax1.set_ylim(0,1)
    ax1.grid(True)

    if 'diffFluxU_x' in df.columns:
        # create a twin x axis to be shared
        ax2 = ax1.twinx()

        ax2.plot(df['x'], df['diffFluxU_x'],'*-', color='yellow', label=u1_flux_legend)
        
        if 'diffFluxU2_x' in df.columns:
            ax2.plot(df['x'], df['diffFluxU2_x'],'*--', color='yellow', label=u2_flux_legend)

        ax2.set_ylabel(r"$q_h(x)$ [g/cm$^2$-s]", fontsize=16, color='yellow')
        ax2.tick_params(axis='y', labelcolor='yellow', labelsize=14)
        ax2.legend(loc='center right', fontsize=12)
        #ax2.set_ylim(0,2)
        #ax2.grid(True)

    plt.title(title, fontsize=20)
    fig.tight_layout()  # otherwise the right y-label is slightly clipped
    plt.show()
    print('')

Problem Statement ¶

There exists two kinds of coupling of unknown variables (or fields), one kind is when the fields share the same domain $\Omega$ and are coupled tightly over the entire $\Omega$. This is the kind that will be studied here. The other kind is when the fields do not share the same domain $\Omega$ and the coupling takes place on a portion of the boundary of $\Omega$ shared by the fields. This latter case will be described in a future notebook.

The following sections describe an extension of the Peclet problem described in previous notebooks for the case when two fields are coupled through a source term.

Strong Form of Problem Statement ¶

Solve the Peclet model problem. Find $u_1:[a,b]\subset\Reals\rightarrow\Reals$ and $u_2:[a,b]\subset\Reals\rightarrow\Reals$ for $D_1 > 0$ and $D_2 > 0$ such that:

\begin{align*} v\, u_1' &= -\bigl(-D_1\, u_1'\bigr)'(x) + S(u_1, u_2) \quad \forall \quad x\in [a,b], \\ u_1(a) &= A_1, \\ u_1(b) &= B_1. \end{align*}

and

\begin{align*} v\, u_2' &= -\bigl(-D_2\, u_2'\bigr)'(x) - S(u_1, u_2) \quad \forall \quad x\in [a,b], \\ u_2(a) &= A_2, \\ u_2(b) &= B_2. \end{align*}

The diffusion flux associated to the quantity $u_i, \, \ i=1,2$ is denoted $q_i := -D_i\,u_i'$, and it is often of interest as a derived quantity. Here a point-wise convective sink (or sweep) is given by $v\,u_i'$. There exists two Peclet numbers:

Peclet number: $\frac{v\,L}{D_i}$.

whose effects has been described in earlier notebooks.

Likewise in the single-field Peclet 1-D problem (Notebook 09), the values of the dependent variables are given on the two end points of the domain (essential boundary conditions or Dirichlet boundary conditions).

Source Coupling ¶

Consider the following source coupling:

\begin{align*} S(u_1, u_2) = S_1 - h_1\,\bigl(u_1-u_1^*\bigr) - \Bigl( S_2 - h_2\, \bigl(u_2 - u_2^*\bigr) \Bigr) \end{align*}

where $S_i$ is a fixed source (sink) for each field, $h_i$ is a transfer coefficient, and $u_i^*$ is a saturation value.

Weak Form of the Problem Statement ¶

The Galerkin weak formulation is as follows. Find $u_1 \in H^1\!\bigl([a,b]\bigr)$ and $u_2 \in H^1\!\bigl([a,b]\bigr)$ so that

\begin{align*} \int\limits_a^b v\, u_1'(x)\, w(x)\,dx + \int\limits_a^b D\, u_1'(x)\,w'(x)\,dx - \int\limits_a^b S(u_1, u_2)\,w(x)\,dx &= 0 \quad \forall \quad w \in H^1_0\!\bigl([a,b]\bigr), \text{and} \\ \int\limits_a^b v\, u_2'(x)\, w(x)\,dx + \int\limits_a^b D\, u_2'(x)\,w'(x)\,dx + \int\limits_a^b S(u_1, u_2)\,w(x)\,dx &= 0 \quad \forall \quad w \in H^1_0\!\bigl([a,b]\bigr), \end{align*}

where $H^1\!\bigl([a,b]\bigr) := \bigl\{ u:[a,b]\subset\Reals\rightarrow \Reals \mid \int_a^b u'^2\,dx < \infty\bigr\}$ and $H^1_0\!\bigl([a,b]\bigr) := \bigl\{ w \mid w \in H^1(a,b), w(a) = 0, w(b) =0 \bigr\}$. Both function sets as just defined are Hilbert spaces. The function $w$ is called a test function. Because $w$, $u_1$, $u_2$ are sought in very similar sets of functions, this weak form is called Galerkin's weak form.

The new form of the source term is the key term to be computed here. Since the MOOSE framework performs the integration, and provides the implementation of the test function, we need to provide the integrand of the integral, that is, the kernel. Therefore the kernel needed is an expansion of what has been covered so far in the course:

$\pm S(u_1, u_2)\,w(x)$.

The kernels are to be evaluated at quadrature points provided by the MOOSE framework.

Problem Solution ¶

We will leverage the Peclet 1D development (Notebook 09) to modify the source term where the coupling takes place.

Code a Source Term Kernel ¶

Starting from the Peclet 1D development, the source kernel file SourceTerm.h can be modified as follows:

/// The variables associated to the source const Real _sourceS; const Real _transferCoeff; const Real _saturation;

const VariableValue & _uCoupled; const Real _sourceSCoupled; const Real _transferCoeffCoupled; const Real _saturationCoupled;

1. Likewise the `SourceTerm.C` class implementation can be also changed to:

Residual:
  ```c++
Real
SourceTerm::computeQpResidual()
{
 return  ( (_sourceS - _transferCoeff * (_u[_qp] - _saturation))              \
       -(_sourceSCoupled -_transferCoeffCoupled * (_uCoupled[_qp] - _saturationCoupled)) )  \
      * _test[_i][_qp];
}

Jacobian diagonal:

Real
SourceTerm::computeQpJacobian()
{
 return (- _transferCoeff + _transferCoeffCoupled) * _phi[_j][_qp] * _test[_i][_qp];
}

Compile and Link Application ¶

Make sure you are in the problem project:

cd ../..
pwd

You should see: ..../engy5310p1
Compile and link the application

make

If all is sucessfull you should see among other things in the screen output:

Linking Library .../engy5310p1/lib/libengy5310p1-opt.la...
Linking Executable .../engy5310p1/engy5310p1-opt...

Input File ¶

The newly revised source kernel can now be used for both fields:

  [source-term-1]
    type = SourceTerm
    variable = u1     # add to produced quantity
    sourceS = ${replace source_s_1}
    transferCoeff = ${replace source_transfer_coeff_1}
    saturation = ${replace source_saturation_1}
    coupledVariable = u2
    sourceSCoupled = ${replace source_s_2}
    transferCoeffCoupled = ${replace source_transfer_coeff_2}
    saturationCoupled = ${replace source_saturation_2}
  []
  [source-term-2]
    type = SourceTerm
    variable = u2     # add to produced quantity
    sourceS = ${replace source_s_2}
    transferCoeff = ${replace source_transfer_coeff_2}
    saturation = ${replace source_saturation_2}
    coupledVariable = u1
    sourceSCoupled = ${replace source_s_1}
    transferCoeffCoupled = ${replace source_transfer_coeff_1}
    saturationCoupled = ${replace source_saturation_1}
  []

On the working examples below a cleaner version with defaults is used to avoid writing entries with zero values. Additional blocks in the input file can be used to create the weak form for both variables as described below.

Run Application ¶

In the engy5310p1/ directory run the application with the Linux shell command:

./engy5310p1-opt -i input.hit

Compare your results with this notebook results below.
Return here to follow instructions on how to implement the calculation of the total energy.

Linear Lagrange FEM ¶

Solve problem with parameter values:

$a = 0$ cm

$b = 25$ cm

Pe$_\text{ave} = 10$

$u_1$ Parameter	Value	$u_2$ Parameter	Value
$A_1$	3 g/cc	$A_2$	0 g/cc
$B_1$	0 g/cc	$B_2$	0 g/cc
$D_1$	0.1 cm^2/s	$D_2$	0.5 cm^2/s
$S_1$	$1\times 10^{-3}$ g/cc-s	$S_2$	0 g/cc-s
$h_1$	$5\times 10^{-3}$ cm/s	$h_2$	0 cm/s
$u_1^*$	1.5 g/cc	$u_2^*$	0 g/cc

FEM parameters:

Basis Functions: First Order Lagrangian

num. of finite elements: 20

Comment: This data tests the case when there is only transfer from $u_1$ to $u_2$. Other cases will be presented in another notebook.

In [2]:

'''Domain'''

x_a = 0
x_b = 25

x_length = x_b - x_a

In [3]:

'''Parameters and data'''

Pe_ave = 10 # mildly convective dominated

diff_coeff_1 = 0.1
source_s_1 = 1e-3
source_transfer_coeff_1 = 5e-3
source_saturation_1 = 1.0

diff_coeff_2 = 0.5

velocity = (Pe_ave * (diff_coeff_1+diff_coeff_2)/2/x_length, 0, 0)  # length scale is the x length

u_a = 3
u_b = 0

u2_a = 0
u2_b = 0

In [4]:

'''FEM Solution'''

n_felem = 20

order = 'first'

n_plot_pts = n_felem + 1

try:    
    from engy_5310.toolkit import write_engy5310_p1_1d_input_file  
except ModuleNotFoundError:
    assert False, 'You need to provide your own code here. Bailing out.'

write_engy5310_p1_1d_input_file(x_left=x_a, x_right=x_b, 
                                u_left=u_a, u_right=u_b, 
                                diff_coeff=diff_coeff_1, source_s=source_s_1,
                                source_transfer_coeff=source_transfer_coeff_1, 
                                source_saturation=source_saturation_1,
                                u2_left=u2_a, u2_right=u2_b, 
                                diff_coeff_2=diff_coeff_2,
                                velocity=velocity, 
                                n_felem=n_felem, order=order, 
                                n_plot_pts=n_plot_pts,
                                compute_diffusion_flux=True,
                                solver='fdp-newt-full')

In [5]:

'''Display MOOSE input file created'''

!cat engy5310p1/input.hit

# Engy-5310 Problem 1: Poisson 1D FEM
# UMass Lowell Nuclear Chemical Engineering
# Prof. Valmor F. de Almeida
# 02May21 18:24:20

# Parameters
xmin = 0.00000e+00
xmax = 2.50000e+01
diff_coeff = 1.00000e-01
source_s = 1.00000e-03
source_transfer_coeff = 5.00000e-03
source_saturation = 1.00000e+00
u_left = 3.00000e+00
u_right = 0.00000e+00
u2_left = 0.00000e+00
u2_right = 0.00000e+00
diff_coeff_2 = 5.00000e-01
velocity = '1.20000e-01 0.00000e+00 0.00000e+00'

[Problem]
  type = FEProblem
  coord_type = XYZ
[]

[Mesh]
  [omega-1d]
    type = GeneratedMeshGenerator
    dim = 1
    xmin = ${replace xmin}
    xmax = ${replace xmax}
    nx = 20
  []
[]

[Variables]
  [u]
    order = first
    family = lagrange
    initial_condition = ${fparse (u_right+u_left)/2}
  []
  [u2]
    order = first
    family = lagrange
    initial_condition = ${fparse (u2_right+u2_left)/2}
  []
[]

[AuxVariables]
  [diffFluxU]
    order = CONSTANT
    family = MONOMIAL_VEC
  []
  [diffFluxU2]
    order = CONSTANT
    family = MONOMIAL_VEC
  []
  [diffFluxU_x]
    order = CONSTANT
    family = MONOMIAL
  []
  [diffFluxU2_x]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [diffusion-term]
    type = DiffusionTerm
    variable = u     # produced quantity
    diffCoeff = ${replace diff_coeff}
  []
  [source-term]
    type = SourceTerm
    variable = u     # add to produced quantity
    sourceS = ${replace source_s}
    transferCoeff = ${replace source_transfer_coeff}
    saturation = ${replace source_saturation}
    coupledVariable = u2
  []
  [convection-term]
    type = ConvectionTerm
    variable = u     # produced quantity
    velocity = ${replace velocity}
  []
  [diffusion-term-2]
    type = DiffusionTerm
    variable = u2     # produced quantity
    diffCoeff = ${replace diff_coeff_2}
  []
  [source-term-2]
    type = SourceTerm
    variable = u2     # add to produced quantity
    coupledVariable = u
    sourceSCoupled = ${replace source_s}
    transferCoeffCoupled = ${replace source_transfer_coeff}
    saturationCoupled = ${replace source_saturation}
  []
  [convection-term-2]
    type = ConvectionTerm
    variable = u2     # produced quantity
    velocity = ${replace velocity}
  []
[]

[AuxKernels]
  [diffusion-flux]
    execute_on = timestep_end
    type = DiffusionFlux
    field = u
    diffCoeff = ${replace diff_coeff}
    variable = diffFluxU     # produced quantity
  []
  [diffusion-flux-x]
    execute_on = timestep_end
    type = VectorVariableComponentAux
    variable = diffFluxU_x    # produced quantity
    component = x
    vector_variable = diffFluxU   
  []
  [diffusion-flux-2]
    execute_on = timestep_end
    type = DiffusionFlux
    field = u2
    diffCoeff = ${replace diff_coeff_2}
    variable = diffFluxU2     # produced quantity
  []
  [diffusion-flux-x-2]
    execute_on = timestep_end
    type = VectorVariableComponentAux
    variable = diffFluxU2_x    # produced quantity
    component = x
    vector_variable = diffFluxU2   
  []
[]

[BCs]
  [entry-u]
    type = DirichletBC
    variable = u
    boundary = left
    value = ${replace u_left}
  []
  [entry-u2]
    type = DirichletBC
    variable = u2
    boundary = left
    value = ${replace u2_left}
  []
  [exit-u]
    type = DirichletBC
    variable = u
    boundary = right
    value = ${replace u_right}
  []
  [exit-u2]
    type = DirichletBC
    variable = u2
    boundary = right
    value = ${replace u2_right}
  []
[]

[Preconditioning]
  active = 'fdp-newt-full'
  [fdp-newt-full]
    type = FDP
    full = true
    solve_type = 'NEWTON'
    petsc_options_iname = '-pc_type -mat_fd_coloring_err -mat_fd_type'
    petsc_options_value = 'lu   1.000e-06               ds'
  []
[]

[Executioner]
  type = Steady
[]

[VectorPostprocessors]
  [omega-data]
    type = LineValueSampler
    execute_on = 'timestep_end final'
    variable = 'u u2 diffFluxU_x diffFluxU2_x'    # output data
    start_point = '${replace xmin} 0 0'
    end_point = '${replace xmax} 0 0'
    num_points = 21
    sort_by = id
  []
[]

[Outputs]
  console = true
  [csv]
    type = CSV
    file_base = 'output'
    execute_on = 'final'
  []
[]

In [6]:

'''Run Engy5310P1 MOOSE App'''

!engy5310p1/engy5310p1-opt -i engy5310p1/input.hit


Framework Information:
MOOSE Version:           git commit a7f499ed31 on 2021-04-30
LibMesh Version:         27141d18f3137f77e33cdb3d565fd38ebfbfc46f
PETSc Version:           3.15.0
SLEPc Version:           3.14.2
Current Time:            Sun May  2 18:24:20 2021
Executable Timestamp:    Sun May  2 16:17:43 2021

Parallelism:
  Num Processors:          1
  Num Threads:             1

Mesh: 
  Parallel Type:           replicated
  Mesh Dimension:          1
  Spatial Dimension:       1
  Nodes:                   
    Total:                 21
    Local:                 21
  Elems:                   
    Total:                 20
    Local:                 20
  Num Subdomains:          1
  Num Partitions:          1

Nonlinear System:
  Num DOFs:                42
  Num Local DOFs:          42
  Variables:               { "u" "u2" } 
  Finite Element Types:    "LAGRANGE" 
  Approximation Orders:    "FIRST" 

Auxiliary System:
  Num DOFs:                80
  Num Local DOFs:          80
  Variables:               { "diffFluxU" "diffFluxU2" } { "diffFluxU_x" "diffFluxU2_x" } 
  Finite Element Types:    "MONOMIAL_VEC" "MONOMIAL" 
  Approximation Orders:    "CONSTANT" "CONSTANT" 

Execution Information:
  Executioner:             Steady
  Solver Mode:             NEWTON
  PETSc Preconditioner:    lu 
  MOOSE Preconditioner:    FDP

 0 Nonlinear |R| = 2.090894e-01
      0 Linear |R| = 2.090894e-01
      1 Linear |R| = 3.244198e-16
 1 Nonlinear |R| = 1.387873e-08
      0 Linear |R| = 1.387873e-08
      1 Linear |R| = 3.268868e-23
 2 Nonlinear |R| = 1.982895e-16
 Solve Converged!
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There is one unused database option. It is:
Option left: name:-i value: engy5310p1/input.hit

In [7]:

'''Show FEM Solution'''

import pandas as pd
df = pd.read_csv('output_omega-data_0002.csv')
    
plot_solution(df, title='Peclet Coupled Variables w/ Dirichlet BC FEM Solution', 
              u1_legend=r'$u_1$ Linear Lagrange', u2_legend=r'$u_2$ Linear Lagrange',
              u1_flux_legend=r'$u_1$ Diff. Flux Constant Monomial',
              u2_flux_legend=r'$u_2$ Diff. Flux Constant Monomial')

Comments:

The concentration of $u_1$ decreases as $u_2$ increases since there is only transfer from 1 to 2.

FEM Error ¶

In [8]:

'''Error Compared to Exact Dimensionless Solution'''

'''coming...'''

Out[8]:

'coming...'

Comments:

Quadratic Lagrange FEM ¶

Solve problem with the same parameter values as above.

FEM parameters:

Basis Functions: Second Order Lagrangian

num. of finite elements: 20

In [9]:

'''FEM Solution'''

n_felem = 20

order = 'second'

n_plot_pts = 2*n_felem + 1

try:    
    from engy_5310.toolkit import write_engy5310_p1_1d_input_file  
except ModuleNotFoundError:
    assert False, 'You need to provide your own code here. Bailing out.'

write_engy5310_p1_1d_input_file(x_left=x_a, x_right=x_b, 
                                u_left=u_a, u_right=u_b, 
                                diff_coeff=diff_coeff_1, source_s=source_s_1,
                                source_transfer_coeff=source_transfer_coeff_1, 
                                source_saturation=source_saturation_1,
                                u2_left=u2_a, u2_right=u2_b, 
                                diff_coeff_2=diff_coeff_2,
                                velocity=velocity, 
                                n_felem=n_felem, order=order, 
                                n_plot_pts=n_plot_pts,
                                compute_diffusion_flux=True,
                                solver='fdp-newt-full')

In [10]:

'''Display MOOSE input file created'''

!cat engy5310p1/input.hit

# Engy-5310 Problem 1: Poisson 1D FEM
# UMass Lowell Nuclear Chemical Engineering
# Prof. Valmor F. de Almeida
# 02May21 18:24:21

# Parameters
xmin = 0.00000e+00
xmax = 2.50000e+01
diff_coeff = 1.00000e-01
source_s = 1.00000e-03
source_transfer_coeff = 5.00000e-03
source_saturation = 1.00000e+00
u_left = 3.00000e+00
u_right = 0.00000e+00
u2_left = 0.00000e+00
u2_right = 0.00000e+00
diff_coeff_2 = 5.00000e-01
velocity = '1.20000e-01 0.00000e+00 0.00000e+00'

[Problem]
  type = FEProblem
  coord_type = XYZ
[]

[Mesh]
  [omega-1d]
    type = GeneratedMeshGenerator
    dim = 1
    xmin = ${replace xmin}
    xmax = ${replace xmax}
    nx = 20
    elem_type = edge3
  []
[]

[Variables]
  [u]
    order = second
    family = lagrange
    initial_condition = ${fparse (u_right+u_left)/2}
  []
  [u2]
    order = second
    family = lagrange
    initial_condition = ${fparse (u2_right+u2_left)/2}
  []
[]

[AuxVariables]
  [diffFluxU]
    order = FIRST
    family = MONOMIAL_VEC
  []
  [diffFluxU2]
    order = FIRST
    family = MONOMIAL_VEC
  []
  [diffFluxU_x]
    order = FIRST
    family = MONOMIAL
  []
  [diffFluxU2_x]
    order = FIRST
    family = MONOMIAL
  []
[]

[Kernels]
  [diffusion-term]
    type = DiffusionTerm
    variable = u     # produced quantity
    diffCoeff = ${replace diff_coeff}
  []
  [source-term]
    type = SourceTerm
    variable = u     # add to produced quantity
    sourceS = ${replace source_s}
    transferCoeff = ${replace source_transfer_coeff}
    saturation = ${replace source_saturation}
    coupledVariable = u2
  []
  [convection-term]
    type = ConvectionTerm
    variable = u     # produced quantity
    velocity = ${replace velocity}
  []
  [diffusion-term-2]
    type = DiffusionTerm
    variable = u2     # produced quantity
    diffCoeff = ${replace diff_coeff_2}
  []
  [source-term-2]
    type = SourceTerm
    variable = u2     # add to produced quantity
    coupledVariable = u
    sourceSCoupled = ${replace source_s}
    transferCoeffCoupled = ${replace source_transfer_coeff}
    saturationCoupled = ${replace source_saturation}
  []
  [convection-term-2]
    type = ConvectionTerm
    variable = u2     # produced quantity
    velocity = ${replace velocity}
  []
[]

[AuxKernels]
  [diffusion-flux]
    execute_on = timestep_end
    type = DiffusionFlux
    field = u
    diffCoeff = ${replace diff_coeff}
    variable = diffFluxU     # produced quantity
  []
  [diffusion-flux-x]
    execute_on = timestep_end
    type = VectorVariableComponentAux
    variable = diffFluxU_x    # produced quantity
    component = x
    vector_variable = diffFluxU   
  []
  [diffusion-flux-2]
    execute_on = timestep_end
    type = DiffusionFlux
    field = u2
    diffCoeff = ${replace diff_coeff_2}
    variable = diffFluxU2     # produced quantity
  []
  [diffusion-flux-x-2]
    execute_on = timestep_end
    type = VectorVariableComponentAux
    variable = diffFluxU2_x    # produced quantity
    component = x
    vector_variable = diffFluxU2   
  []
[]

[BCs]
  [entry-u]
    type = DirichletBC
    variable = u
    boundary = left
    value = ${replace u_left}
  []
  [entry-u2]
    type = DirichletBC
    variable = u2
    boundary = left
    value = ${replace u2_left}
  []
  [exit-u]
    type = DirichletBC
    variable = u
    boundary = right
    value = ${replace u_right}
  []
  [exit-u2]
    type = DirichletBC
    variable = u2
    boundary = right
    value = ${replace u2_right}
  []
[]

[Preconditioning]
  active = 'fdp-newt-full'
  [fdp-newt-full]
    type = FDP
    full = true
    solve_type = 'NEWTON'
    petsc_options_iname = '-pc_type -mat_fd_coloring_err -mat_fd_type'
    petsc_options_value = 'lu   1.000e-06               ds'
  []
[]

[Executioner]
  type = Steady
[]

[VectorPostprocessors]
  [omega-data]
    type = LineValueSampler
    execute_on = 'timestep_end final'
    variable = 'u u2 diffFluxU_x diffFluxU2_x'    # output data
    start_point = '${replace xmin} 0 0'
    end_point = '${replace xmax} 0 0'
    num_points = 41
    sort_by = id
  []
[]

[Outputs]
  console = true
  [csv]
    type = CSV
    file_base = 'output'
    execute_on = 'final'
  []
[]

In [11]:

'''Run Engy5310P1 MOOSE App'''

!engy5310p1/engy5310p1-opt -i engy5310p1/input.hit


Framework Information:
MOOSE Version:           git commit a7f499ed31 on 2021-04-30
LibMesh Version:         27141d18f3137f77e33cdb3d565fd38ebfbfc46f
PETSc Version:           3.15.0
SLEPc Version:           3.14.2
Current Time:            Sun May  2 18:24:22 2021
Executable Timestamp:    Sun May  2 16:17:43 2021

Parallelism:
  Num Processors:          1
  Num Threads:             1

Mesh: 
  Parallel Type:           replicated
  Mesh Dimension:          1
  Spatial Dimension:       1
  Nodes:                   
    Total:                 41
    Local:                 41
  Elems:                   
    Total:                 20
    Local:                 20
  Num Subdomains:          1
  Num Partitions:          1

Nonlinear System:
  Num DOFs:                82
  Num Local DOFs:          82
  Variables:               { "u" "u2" } 
  Finite Element Types:    "LAGRANGE" 
  Approximation Orders:    "SECOND" 

Auxiliary System:
  Num DOFs:                160
  Num Local DOFs:          160
  Variables:               { "diffFluxU" "diffFluxU2" } { "diffFluxU_x" "diffFluxU2_x" } 
  Finite Element Types:    "MONOMIAL_VEC" "MONOMIAL" 
  Approximation Orders:    "FIRST" "FIRST" 

Execution Information:
  Executioner:             Steady
  Solver Mode:             NEWTON
  PETSc Preconditioner:    lu 
  MOOSE Preconditioner:    FDP

 0 Nonlinear |R| = 4.871389e-01
      0 Linear |R| = 4.871389e-01
      1 Linear |R| = 9.612425e-16
 1 Nonlinear |R| = 3.801380e-07
      0 Linear |R| = 3.801380e-07
      1 Linear |R| = 8.456680e-21
 2 Nonlinear |R| = 6.874920e-15
 Solve Converged!
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There is one unused database option. It is:
Option left: name:-i value: engy5310p1/input.hit

In [12]:

'''Show FEM Solution'''

import pandas as pd
df = pd.read_csv('output_omega-data_0002.csv')
    
plot_solution(df, title='Peclet Coupled Variables w/ Dirichlet BC FEM Solution', 
              u1_legend=r'$u_1$ Quadratic Lagrange', u2_legend=r'$u_2$ Quadratic Lagrange',
              u1_flux_legend=r'$u_1$ Diff. Flux Linear Monomial',
              u2_flux_legend=r'$u_2$ Diff. Flux Linear Monomial')

FEM Error ¶

In [13]:

'''Error Compared to Exact Dimensionless Solution'''

'''comming...'''

Out[13]:

'comming...'

Comments:

High Peclet Number ¶

Highly convective problems may lead to numerical difficulties. Below, enough finite elements are chosen so to avoid oscillations.

Quadratic Lagrange FEM ¶

Solve problem with the same parameter values above.

FEM parameters:

Basis Functions: Second Order Lagrangian

num. of finite elements: 10

In [14]:

'''Parameters and data'''

Pe_ave = 50 # convective dominated

velocity = (Pe_ave * (diff_coeff_1+diff_coeff_2)/2/x_length, 0, 0)  # length scale is the x length

In [15]:

'''FEM Solution'''

n_felem = 10

order = 'second'

n_plot_pts = 2*n_felem + 1

try:    
    from engy_5310.toolkit import write_engy5310_p1_1d_input_file  
except ModuleNotFoundError:
    assert False, 'You need to provide your own code here. Bailing out.'

write_engy5310_p1_1d_input_file(x_left=x_a, x_right=x_b, 
                                u_left=u_a, u_right=u_b, 
                                diff_coeff=diff_coeff_1, source_s=source_s_1,
                                source_transfer_coeff=source_transfer_coeff_1, 
                                source_saturation=source_saturation_1,
                                u2_left=u2_a, u2_right=u2_b, 
                                diff_coeff_2=diff_coeff_2,
                                velocity=velocity, 
                                n_felem=n_felem, order=order, 
                                n_plot_pts=n_plot_pts,
                                compute_diffusion_flux=True,
                                solver='fdp-newt-full')

In [16]:

'''Display MOOSE input file created'''

!cat engy5310p1/input.hit

# Engy-5310 Problem 1: Poisson 1D FEM
# UMass Lowell Nuclear Chemical Engineering
# Prof. Valmor F. de Almeida
# 02May21 18:24:22

# Parameters
xmin = 0.00000e+00
xmax = 2.50000e+01
diff_coeff = 1.00000e-01
source_s = 1.00000e-03
source_transfer_coeff = 5.00000e-03
source_saturation = 1.00000e+00
u_left = 3.00000e+00
u_right = 0.00000e+00
u2_left = 0.00000e+00
u2_right = 0.00000e+00
diff_coeff_2 = 5.00000e-01
velocity = '6.00000e-01 0.00000e+00 0.00000e+00'

[Problem]
  type = FEProblem
  coord_type = XYZ
[]

[Mesh]
  [omega-1d]
    type = GeneratedMeshGenerator
    dim = 1
    xmin = ${replace xmin}
    xmax = ${replace xmax}
    nx = 10
    elem_type = edge3
  []
[]

[Variables]
  [u]
    order = second
    family = lagrange
    initial_condition = ${fparse (u_right+u_left)/2}
  []
  [u2]
    order = second
    family = lagrange
    initial_condition = ${fparse (u2_right+u2_left)/2}
  []
[]

[AuxVariables]
  [diffFluxU]
    order = FIRST
    family = MONOMIAL_VEC
  []
  [diffFluxU2]
    order = FIRST
    family = MONOMIAL_VEC
  []
  [diffFluxU_x]
    order = FIRST
    family = MONOMIAL
  []
  [diffFluxU2_x]
    order = FIRST
    family = MONOMIAL
  []
[]

[Kernels]
  [diffusion-term]
    type = DiffusionTerm
    variable = u     # produced quantity
    diffCoeff = ${replace diff_coeff}
  []
  [source-term]
    type = SourceTerm
    variable = u     # add to produced quantity
    sourceS = ${replace source_s}
    transferCoeff = ${replace source_transfer_coeff}
    saturation = ${replace source_saturation}
    coupledVariable = u2
  []
  [convection-term]
    type = ConvectionTerm
    variable = u     # produced quantity
    velocity = ${replace velocity}
  []
  [diffusion-term-2]
    type = DiffusionTerm
    variable = u2     # produced quantity
    diffCoeff = ${replace diff_coeff_2}
  []
  [source-term-2]
    type = SourceTerm
    variable = u2     # add to produced quantity
    coupledVariable = u
    sourceSCoupled = ${replace source_s}
    transferCoeffCoupled = ${replace source_transfer_coeff}
    saturationCoupled = ${replace source_saturation}
  []
  [convection-term-2]
    type = ConvectionTerm
    variable = u2     # produced quantity
    velocity = ${replace velocity}
  []
[]

[AuxKernels]
  [diffusion-flux]
    execute_on = timestep_end
    type = DiffusionFlux
    field = u
    diffCoeff = ${replace diff_coeff}
    variable = diffFluxU     # produced quantity
  []
  [diffusion-flux-x]
    execute_on = timestep_end
    type = VectorVariableComponentAux
    variable = diffFluxU_x    # produced quantity
    component = x
    vector_variable = diffFluxU   
  []
  [diffusion-flux-2]
    execute_on = timestep_end
    type = DiffusionFlux
    field = u2
    diffCoeff = ${replace diff_coeff_2}
    variable = diffFluxU2     # produced quantity
  []
  [diffusion-flux-x-2]
    execute_on = timestep_end
    type = VectorVariableComponentAux
    variable = diffFluxU2_x    # produced quantity
    component = x
    vector_variable = diffFluxU2   
  []
[]

[BCs]
  [entry-u]
    type = DirichletBC
    variable = u
    boundary = left
    value = ${replace u_left}
  []
  [entry-u2]
    type = DirichletBC
    variable = u2
    boundary = left
    value = ${replace u2_left}
  []
  [exit-u]
    type = DirichletBC
    variable = u
    boundary = right
    value = ${replace u_right}
  []
  [exit-u2]
    type = DirichletBC
    variable = u2
    boundary = right
    value = ${replace u2_right}
  []
[]

[Preconditioning]
  active = 'fdp-newt-full'
  [fdp-newt-full]
    type = FDP
    full = true
    solve_type = 'NEWTON'
    petsc_options_iname = '-pc_type -mat_fd_coloring_err -mat_fd_type'
    petsc_options_value = 'lu   1.000e-06               ds'
  []
[]

[Executioner]
  type = Steady
[]

[VectorPostprocessors]
  [omega-data]
    type = LineValueSampler
    execute_on = 'timestep_end final'
    variable = 'u u2 diffFluxU_x diffFluxU2_x'    # output data
    start_point = '${replace xmin} 0 0'
    end_point = '${replace xmax} 0 0'
    num_points = 21
    sort_by = id
  []
[]

[Outputs]
  console = true
  [csv]
    type = CSV
    file_base = 'output'
    execute_on = 'final'
  []
[]

In [17]:

'''Run Engy5310P1 MOOSE App'''

!engy5310p1/engy5310p1-opt -i engy5310p1/input.hit


Framework Information:
MOOSE Version:           git commit a7f499ed31 on 2021-04-30
LibMesh Version:         27141d18f3137f77e33cdb3d565fd38ebfbfc46f
PETSc Version:           3.15.0
SLEPc Version:           3.14.2
Current Time:            Sun May  2 18:24:23 2021
Executable Timestamp:    Sun May  2 16:17:43 2021

Parallelism:
  Num Processors:          1
  Num Threads:             1

Mesh: 
  Parallel Type:           replicated
  Mesh Dimension:          1
  Spatial Dimension:       1
  Nodes:                   
    Total:                 21
    Local:                 21
  Elems:                   
    Total:                 10
    Local:                 10
  Num Subdomains:          1
  Num Partitions:          1

Nonlinear System:
  Num DOFs:                42
  Num Local DOFs:          42
  Variables:               { "u" "u2" } 
  Finite Element Types:    "LAGRANGE" 
  Approximation Orders:    "SECOND" 

Auxiliary System:
  Num DOFs:                80
  Num Local DOFs:          80
  Variables:               { "diffFluxU" "diffFluxU2" } { "diffFluxU_x" "diffFluxU2_x" } 
  Finite Element Types:    "MONOMIAL_VEC" "MONOMIAL" 
  Approximation Orders:    "FIRST" "FIRST" 

Execution Information:
  Executioner:             Steady
  Solver Mode:             NEWTON
  PETSc Preconditioner:    lu 
  MOOSE Preconditioner:    FDP

 0 Nonlinear |R| = 9.005892e-01
      0 Linear |R| = 9.005892e-01
      1 Linear |R| = 1.085705e-15
 1 Nonlinear |R| = 1.660617e-07
      0 Linear |R| = 1.660617e-07
      1 Linear |R| = 1.928338e-22
 2 Nonlinear |R| = 4.792070e-16
 Solve Converged!
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There is one unused database option. It is:
Option left: name:-i value: engy5310p1/input.hit

In [18]:

'''Show FEM Solution'''

import pandas as pd
df = pd.read_csv('output_omega-data_0002.csv')
    
plot_solution(df, title='Peclet Coupled Variables w/ Dirichlet BC FEM Solution', 
              u1_legend=r'$u_1$ Quadratic Lagrange', u2_legend=r'$u_2$ Quadratic Lagrange',
              u1_flux_legend=r'$u_1$ Diff. Flux Linear Monomial',
              u2_flux_legend=r'$u_2$ Diff. Flux Linear Monomial')

Comments:

As convection dominates, the boundary layer sharpens and numerical instability appears in the form of oscillations. Either an increase of number of finite element basis functions or mesh adaptivity could resolve the boundary layer as demonstrated below.

Compute Error ¶

In [19]:

'''Flux Error Compared to Exact Dimensionless Solution'''

'''coming...'''

Out[19]:

'coming...'

Quadratic Lagrange FEM ¶

Solve problem with the same parameter values above.

FEM parameters:

Basis Functions: Second Order Lagrangian

num. of finite elements: 20

In [20]:

'''FEM Solution'''

n_felem = 20

x_bias = .5  # mesh adaptivity

order = 'second'

n_plot_pts = 2*n_felem + 1
n_plot_pts *= 2

try:    
    from engy_5310.toolkit import write_engy5310_p1_1d_input_file  
except ModuleNotFoundError:
    assert False, 'You need to provide your own code here. Bailing out.'

write_engy5310_p1_1d_input_file(x_left=x_a, x_right=x_b, 
                                u_left=u_a, u_right=u_b, 
                                diff_coeff=diff_coeff_1, source_s=source_s_1,
                                source_transfer_coeff=source_transfer_coeff_1, 
                                source_saturation=source_saturation_1,
                                u2_left=u2_a, u2_right=u2_b, 
                                diff_coeff_2=diff_coeff_2,
                                velocity=velocity, 
                                n_felem=n_felem, order=order, 
                                x_bias=x_bias,
                                n_plot_pts=n_plot_pts,
                                compute_diffusion_flux=True,
                                solver='fdp-newt-full')

In [21]:

'''Display MOOSE input file created'''

!cat engy5310p1/input.hit

# Engy-5310 Problem 1: Poisson 1D FEM
# UMass Lowell Nuclear Chemical Engineering
# Prof. Valmor F. de Almeida
# 02May21 18:24:23

# Parameters
xmin = 0.00000e+00
xmax = 2.50000e+01
diff_coeff = 1.00000e-01
source_s = 1.00000e-03
source_transfer_coeff = 5.00000e-03
source_saturation = 1.00000e+00
u_left = 3.00000e+00
u_right = 0.00000e+00
u2_left = 0.00000e+00
u2_right = 0.00000e+00
diff_coeff_2 = 5.00000e-01
velocity = '6.00000e-01 0.00000e+00 0.00000e+00'

[Problem]
  type = FEProblem
  coord_type = XYZ
[]

[Mesh]
  [omega-1d]
    type = GeneratedMeshGenerator
    dim = 1
    xmin = ${replace xmin}
    xmax = ${replace xmax}
    nx = 20
    elem_type = edge3
    bias_x = 5.000e-01
  []
[]

[Variables]
  [u]
    order = second
    family = lagrange
    initial_condition = ${fparse (u_right+u_left)/2}
  []
  [u2]
    order = second
    family = lagrange
    initial_condition = ${fparse (u2_right+u2_left)/2}
  []
[]

[AuxVariables]
  [diffFluxU]
    order = FIRST
    family = MONOMIAL_VEC
  []
  [diffFluxU2]
    order = FIRST
    family = MONOMIAL_VEC
  []
  [diffFluxU_x]
    order = FIRST
    family = MONOMIAL
  []
  [diffFluxU2_x]
    order = FIRST
    family = MONOMIAL
  []
[]

[Kernels]
  [diffusion-term]
    type = DiffusionTerm
    variable = u     # produced quantity
    diffCoeff = ${replace diff_coeff}
  []
  [source-term]
    type = SourceTerm
    variable = u     # add to produced quantity
    sourceS = ${replace source_s}
    transferCoeff = ${replace source_transfer_coeff}
    saturation = ${replace source_saturation}
    coupledVariable = u2
  []
  [convection-term]
    type = ConvectionTerm
    variable = u     # produced quantity
    velocity = ${replace velocity}
  []
  [diffusion-term-2]
    type = DiffusionTerm
    variable = u2     # produced quantity
    diffCoeff = ${replace diff_coeff_2}
  []
  [source-term-2]
    type = SourceTerm
    variable = u2     # add to produced quantity
    coupledVariable = u
    sourceSCoupled = ${replace source_s}
    transferCoeffCoupled = ${replace source_transfer_coeff}
    saturationCoupled = ${replace source_saturation}
  []
  [convection-term-2]
    type = ConvectionTerm
    variable = u2     # produced quantity
    velocity = ${replace velocity}
  []
[]

[AuxKernels]
  [diffusion-flux]
    execute_on = timestep_end
    type = DiffusionFlux
    field = u
    diffCoeff = ${replace diff_coeff}
    variable = diffFluxU     # produced quantity
  []
  [diffusion-flux-x]
    execute_on = timestep_end
    type = VectorVariableComponentAux
    variable = diffFluxU_x    # produced quantity
    component = x
    vector_variable = diffFluxU   
  []
  [diffusion-flux-2]
    execute_on = timestep_end
    type = DiffusionFlux
    field = u2
    diffCoeff = ${replace diff_coeff_2}
    variable = diffFluxU2     # produced quantity
  []
  [diffusion-flux-x-2]
    execute_on = timestep_end
    type = VectorVariableComponentAux
    variable = diffFluxU2_x    # produced quantity
    component = x
    vector_variable = diffFluxU2   
  []
[]

[BCs]
  [entry-u]
    type = DirichletBC
    variable = u
    boundary = left
    value = ${replace u_left}
  []
  [entry-u2]
    type = DirichletBC
    variable = u2
    boundary = left
    value = ${replace u2_left}
  []
  [exit-u]
    type = DirichletBC
    variable = u
    boundary = right
    value = ${replace u_right}
  []
  [exit-u2]
    type = DirichletBC
    variable = u2
    boundary = right
    value = ${replace u2_right}
  []
[]

[Preconditioning]
  active = 'fdp-newt-full'
  [fdp-newt-full]
    type = FDP
    full = true
    solve_type = 'NEWTON'
    petsc_options_iname = '-pc_type -mat_fd_coloring_err -mat_fd_type'
    petsc_options_value = 'lu   1.000e-06               ds'
  []
[]

[Executioner]
  type = Steady
[]

[VectorPostprocessors]
  [omega-data]
    type = LineValueSampler
    execute_on = 'timestep_end final'
    variable = 'u u2 diffFluxU_x diffFluxU2_x'    # output data
    start_point = '${replace xmin} 0 0'
    end_point = '${replace xmax} 0 0'
    num_points = 82
    sort_by = id
  []
[]

[Outputs]
  console = true
  [csv]
    type = CSV
    file_base = 'output'
    execute_on = 'final'
  []
[]

In [22]:

'''Run Engy5310P1 MOOSE App'''

!engy5310p1/engy5310p1-opt -i engy5310p1/input.hit


Framework Information:
MOOSE Version:           git commit a7f499ed31 on 2021-04-30
LibMesh Version:         27141d18f3137f77e33cdb3d565fd38ebfbfc46f
PETSc Version:           3.15.0
SLEPc Version:           3.14.2
Current Time:            Sun May  2 18:24:23 2021
Executable Timestamp:    Sun May  2 16:17:43 2021

Parallelism:
  Num Processors:          1
  Num Threads:             1

Mesh: 
  Parallel Type:           replicated
  Mesh Dimension:          1
  Spatial Dimension:       1
  Nodes:                   
    Total:                 41
    Local:                 41
  Elems:                   
    Total:                 20
    Local:                 20
  Num Subdomains:          1
  Num Partitions:          1

Nonlinear System:
  Num DOFs:                82
  Num Local DOFs:          82
  Variables:               { "u" "u2" } 
  Finite Element Types:    "LAGRANGE" 
  Approximation Orders:    "SECOND" 

Auxiliary System:
  Num DOFs:                160
  Num Local DOFs:          160
  Variables:               { "diffFluxU" "diffFluxU2" } { "diffFluxU_x" "diffFluxU2_x" } 
  Finite Element Types:    "MONOMIAL_VEC" "MONOMIAL" 
  Approximation Orders:    "FIRST" "FIRST" 

Execution Information:
  Executioner:             Steady
  Solver Mode:             NEWTON
  PETSc Preconditioner:    lu 
  MOOSE Preconditioner:    FDP

 0 Nonlinear |R| = 1.690715e+04
      0 Linear |R| = 1.690715e+04
      1 Linear |R| = 4.774424e-12
 1 Nonlinear |R| = 3.952108e-06
 Solve Converged!
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There is one unused database option. It is:
Option left: name:-i value: engy5310p1/input.hit

In [23]:

'''Show FEM Solution'''

import pandas as pd
df = pd.read_csv('output_omega-data_0002.csv')
    
plot_solution(df, title='Peclet Coupled Variables w/ Dirichlet BC FEM Solution', 
              u1_legend=r'$u_1$ Quadratic Lagrange', u2_legend=r'$u_2$ Quadratic Lagrange',
              u1_flux_legend=r'$u_1$ Diff. Flux Linear Monomial',
              u2_flux_legend=r'$u_2$ Diff. Flux Linear Monomial')

Comments:

The mesh adaptivity with refinement on the right side of the domain captures the exit boundary layer. An increase of the number of finite element would also correct the oscillation issue but this strategy is often undesirable.

Compute Error ¶

In [24]:

'''Error Compared to Exact Dimensionless Solution'''

'''coming'''

Out[24]:

'coming'

Application Tree ¶

This tree printout helps the understanding of various pieces of the MOOSE application repository created after all the above steps including future implementations in the notebooks following the present one that covers various boundary conditions.

In [25]:

!tree engy5310p1

engy5310p1
├── LICENSE
├── Makefile
├── README.md
├── __pycache__
│   └── chigger.cpython-38.pyc
├── build
│   ├── header_symlinks
│   │   ├── BoundaryEnergy.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/postprocessors/BoundaryEnergy.h
│   │   ├── BulkEnergy.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/postprocessors/BulkEnergy.h
│   │   ├── ConvectionTerm.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/kernels/ConvectionTerm.h
│   │   ├── DiffusionFlux.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/auxkernels/DiffusionFlux.h
│   │   ├── DiffusionFluxComponent.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/auxkernels/DiffusionFluxComponent.h
│   │   ├── DiffusionTerm.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/kernels/DiffusionTerm.h
│   │   ├── Engy5310P1App.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/base/Engy5310P1App.h
│   │   ├── InterfaceDiffusion.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/interfkernels/InterfaceDiffusion.h
│   │   ├── InterfaceNormalFluxContinuity.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/interfkernels/InterfaceNormalFluxContinuity.h
│   │   ├── InterfacePartition.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/interfkernels/InterfacePartition.h
│   │   ├── NormalFluxBC.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/bcs/NormalFluxBC.h
│   │   └── SourceTerm.h -> /home/dealmeida/OneDrive/uml-courses/engy-5310/2021-01-05-spring/jupynb-repo/notebooks/engy5310p1/include/kernels/SourceTerm.h
│   └── unity_src
│       ├── auxkernels_Unity.C
│       ├── auxkernels_Unity.x86_64-pc-linux-gnu.opt.lo
│       ├── auxkernels_Unity.x86_64-pc-linux-gnu.opt.lo.d
│       ├── bcs_Unity.C
│       ├── bcs_Unity.x86_64-pc-linux-gnu.opt.lo
│       ├── bcs_Unity.x86_64-pc-linux-gnu.opt.lo.d
│       ├── interfkernels_Unity.C
│       ├── interfkernels_Unity.x86_64-pc-linux-gnu.opt.lo
│       ├── interfkernels_Unity.x86_64-pc-linux-gnu.opt.lo.d
│       ├── kernels_Unity.C
│       ├── kernels_Unity.x86_64-pc-linux-gnu.opt.lo
│       ├── kernels_Unity.x86_64-pc-linux-gnu.opt.lo.d
│       ├── postprocessors_Unity.C
│       ├── postprocessors_Unity.x86_64-pc-linux-gnu.opt.lo
│       └── postprocessors_Unity.x86_64-pc-linux-gnu.opt.lo.d
├── engy5310p1-opt
├── include
│   ├── auxkernels
│   │   ├── DiffusionFlux.h
│   │   └── DiffusionFluxComponent.h
│   ├── base
│   │   └── Engy5310P1App.h
│   ├── bcs
│   │   └── NormalFluxBC.h
│   ├── interfkernels
│   │   ├── InterfaceDiffusion.h
│   │   ├── InterfaceNormalFluxContinuity.h
│   │   ├── InterfaceNormalFluxContinuity.h~
│   │   ├── InterfacePartition.h
│   │   └── InterfaceReaction
│   ├── kernels
│   │   ├── ConvectionTerm.h
│   │   ├── ConvectionTerm.h~
│   │   ├── DiffusionTerm.h
│   │   ├── SourceTerm.h
│   │   └── SourceTerm.h~
│   └── postprocessors
│       ├── BoundaryEnergy.h
│       └── BulkEnergy.h
├── input-test.hit
├── input.hit
├── lib
│   ├── libengy5310p1-opt.la
│   ├── libengy5310p1-opt.so -> libengy5310p1-opt.so.0.0.0
│   ├── libengy5310p1-opt.so.0 -> libengy5310p1-opt.so.0.0.0
│   └── libengy5310p1-opt.so.0.0.0
├── mesh_test.i
├── partition-test.i
├── src
│   ├── auxkernels
│   │   ├── DiffusionFlux.C
│   │   └── DiffusionFluxComponent.C
│   ├── base
│   │   ├── Engy5310P1App.C
│   │   ├── Engy5310P1App.x86_64-pc-linux-gnu.opt.lo
│   │   └── Engy5310P1App.x86_64-pc-linux-gnu.opt.lo.d
│   ├── bcs
│   │   ├── NormalFluxBC.C
│   │   └── NormalFluxBC.C~
│   ├── interfkernels
│   │   ├── InterfaceDiffusion.C
│   │   ├── InterfaceNormalFluxContinuity.C
│   │   ├── InterfacePartition.C
│   │   └── InterfaceReaction
│   ├── kernels
│   │   ├── ConvectionTerm.C
│   │   ├── ConvectionTerm.C~
│   │   ├── DiffusionTerm.C
│   │   ├── DiffusionTerm.C~
│   │   ├── SourceTerm.C
│   │   └── SourceTerm.C~
│   ├── main.C
│   ├── main.C~
│   ├── main.x86_64-pc-linux-gnu.opt.lo
│   ├── main.x86_64-pc-linux-gnu.opt.lo.d
│   └── postprocessors
│       ├── BoundaryEnergy.C
│       └── BulkEnergy.C
├── subdom_test2.i
├── test
│   └── lib
│       ├── libengy5310p1_test-opt.la
│       ├── libengy5310p1_test-opt.so -> libengy5310p1_test-opt.so.0.0.0
│       ├── libengy5310p1_test-opt.so.0 -> libengy5310p1_test-opt.so.0.0.0
│       └── libengy5310p1_test-opt.so.0.0.0
└── vtkviz.py

21 directories, 85 files

In [ ]:

12. Peclet Bulk-Coupled 1D with Dirichlet Boundary Conditions¶

Table of Contents¶

Table of Contents ¶