MOOSE 系统介绍

[Kernels]
[diffusion]
type = ADDiffusion # 使用自动微分的拉普拉斯算子
variable = pressure # 对上面的 "pressure" 变量进行操作
[]
[]

[Mesh]
[gmg]
type = GeneratedMeshGenerator # 可以生成简单的线、矩形和长方体
dim = 2 # 网格维度
nx = 100 # x 方向的单元数
ny = 10 # y 方向的单元数
xmax = 0.304 # 测试室长度
ymax = 0.0257 # 测试室半径
[]
coord_type = RZ # 轴对称 RZ
rz_coord_axis = X # 对称轴
[]

[Variables/pressure]

默认添加一个线性拉格朗日变量
[]

[Kernels/diffusion]
type = ADDiffusion # 使用自动微分的拉普拉斯算子
variable = pressure # 对上面的 "pressure" 变量进行操作
[]

[BCs]
[inlet]
type = DirichletBC # 简单的 u=value 边界条件
variable = pressure # 要设置的变量
boundary = left # 网格中的一个边集名称
value = 4000 # (Pa) 来自论文图 2。1mm 球体的第一个数据点。
[]
[outlet]
type = DirichletBC
variable = pressure
boundary = right
value = 0 # (Pa) 根据图 2 给出 1mm 球体的正确压降
[]
[]

[Problem]
type = FEProblem # 这是 MOOSE 中“正常”类型的有限元问题
[]

[Executioner]
type = Steady # 稳态问题
solve_type = NEWTON # 执行牛顿求解，使用 AD 计算雅可比项
petsc_options_iname = '-pc_type -pc_hypre_type' # PETSc 选项对，值在下面
petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
exodus = true # 输出 Exodus 格式
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step01_diffusion
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step1.i

[Kernels]
[diff]
type = Diffusion
variable = u
[]
[]

[BCs]
[right]
type = NeumannBC
variable = u
boundary = 1
value = 1
[]
[]

[Kernels]
[adv_u]
implicit = false
type = ConservativeAdvection
variable = u
velocity = '1 0 0'
[]
[]

[Kernels]
[ffn]
type = BodyForce
variable = u
function = f_fn
[]
[]

[Variables]
[u]
order = FIRST
family = LAGRANGE
block = 0
[]
[v]
order = FIRST
family = LAGRANGE
block = 1
[]
[]

[Executioner]
solve_type = PJFNK

[Executioner]
type = Steady
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'

int x = 2;

{
int x = 5; // "遮蔽"了另一个 x - 不好的做法
assert(x == 5);
}

assert(x == 2);

a = b * (c - 4) / d++;

if (a && b && f()) { e = a; }

a = std::pow(r, 2); // 调用标准 pow() 函数
b = YourLib::pow(r, 2); // 调用 YourLib 命名空间或类中的 pow()

using std::pow; // 现在 "pow" 可以自动表示 "std::pow"
using YourLib::pow; // 或者它可以表示 "YourLib::pow"...

c = pow(r, 2); // 模棱两可，或者根据 r 的类型推断

t = my_obj.someFunction();
b = my_ptr->someFunction();

float pi = 3.14;

int approx_pi = static_cast<int>(pi);

float f = (float) "3.14"; // 无法编译

unsigned int huge_value = 4294967295; // ok
int i = huge_value; // 值被静默更改！
int j = static_cast<int>(huge_value); // 没用！

int i = cast_int<int>(huge_value); // 在非优化版本中会触发断言失败

for (int i=0; i<10; ++i) { foo(i); }
for (auto val : my_container) { foo(val); }
while (boolean-expression) { bar(); }
do { baz(); } while (boolean-expression);

if (boolean-expression) { }
else if (boolean-expression) { }
else { }

make

// "MyObject" 的 "myMethod" 函数的定义
void MyObject::myMethod()
{
std::cout << "Hello, World!\n";
}
MyNamespace::a = 2.718;
MyObject::myMethod();

int *p;

int a;
p = &a;

*p = 5; // 通过 "p" 设置 "a" 的值
std::cout << *p << "\n"; // 打印 5
std::cout << a << "\n"; // 打印 5

int a = 5;
int *p; // 声明一个指针
p = &a; // 将 'p' 设置为 'a' 的地址
*p = *p + 2; // 获取 p 指向的值，加 2，
// 将结果存储在同一位置
std::cout << a << "\n"; // 打印 7
std::cout << *p << "\n"; // 打印 7
std::cout << p << "\n"; // 打印一个地址 (0x7fff5fbfe95c)

p = &a;

p = p + 1000;

*p; // 访问 &a + 1000 处的内存

int a = 5;
int &r = a; // 定义并初始化一个引用
r = r + 2;
std::cout << a << "\n"; // 打印 7
std::cout << r << "\n"; // 打印 7
std::cout << &r << "\n"; // 打印 a 的地址

&r = &r + 1; // 无法编译
int &r; // 无法编译

std::vector<int> four_ints(4);
int &r = four_ints[0];
r = 5; // 有效: four_ints 现在是 {5,0,0,0}
four_ints.clear(); // four_ints 现在是 {}
r = 6; // 未定义行为，鼻中妖魔

Foo & Foo::operator= (const Foo & other) {
this->member1 = other.member1; // 必须拷贝所有东西
this->member2 = other.member2;
}

Foo & Foo::operator= (Foo && other) {
// std::move "将左值变为右值"
this->member1 = std::move(other.member1); // 可以廉价地“窃取”内存
this->member2 = std::move(other.member2);
}

int *iPtr; // 声明
iPtr = &c;
int a = b + *iPtr;

int &iRef = c; // 必须初始化
int a = b + iRef;

std::sqrt(a + b); // "a+b" 创建一个很快就会消失的对象

result = someFunction(a, b, my_shape);

void swap(int a, int b)
{
int temp = a;
a = b;
b = temp;
}
int i = 1;
int j = 2;
swap (i, j); // i 和 j 是参数
std::cout << i << " " << j; // 打印 1 2
// i 和 j 没有被交换

void swap(int &a, int &b)
{
int temp = a;
a = b;
b = temp;
}
int i = 1;
int j = 2;
swap (i, j); // i 和 j 是参数
std::cout << i << " " << j; // 打印 2 1
// i 和 j 被正确交换

{
int a = 4;
int *b = new int; // 动态分配；b 的值是什么？
auto c = std::make_unique<int>(6); // 动态分配
*b = 5;
int d = a + *b + *c;
std::cout << d; // 打印 15
}
// a, b, 和 c 不再在栈上
// *c 被 std::unique_ptr 从堆中删除
// *b 是泄漏的内存 - 我们忘记了 "delete b"，现在为时已晚！

int a;
auto b = std::make_unique<int>(); // 动态分配
int &r = *b; // 创建对新分配对象的引用
a = 4;
r = 5;
int c = a + r;
std::cout << c; // 打印 9

{
std::string name("myObject");
print(name);
}
void print(const std::string & name)
{
name = "MineNow"; // 编译时错误
const_cast<std::string &>(name) = "MineNow"; // 只是糟糕的代码
...
}

constexpr int factorial(int n)
{
if (n <= 1)
return 1;
return n * factorial(n-1);
}
{
constexpr int a = factorial(6); // 直接编译为 a = 720
int b = 6;
function_which_might_modify(b);
int c = factorial(b); // 在运行时计算
}

int foo(int value);
int foo(float value);
int foo(float value, bool is_initialized);
...

template <class T>
T getMax (T a, T b)
{
if (a > b)
return a;
else
return b;
}

template <class T>
T getMax (T a, T b)
{
return (a > b ? a : b); // "三元" 运算符
}
int i = 5, j = 6, k;
float x = 3.142; y = 2.718, z;
k = getMax(i, j); // 使用 int 版本
z = getMax(x, y); // 使用 float 版本
k = getMax<int>(i, j); // 显式调用 int 版本

template<class T>
void print(T value)
{
std::cout << value << std::endl;
}

template<>
void print<bool>(bool value)
{
if (value)
std::cout << "true\n";
else
std::cout << "false\n";
}

int main()
{
int a = 5;
bool b = true;
print(a); // 打印 5
print(b); // 打印 true
}

class Point
{
public:
Point(float x, float y); // 构造函数

// 访问器
float getX();
float getY();
void setX(float x);
void setY(float y);
private:
float _x, _y;
};

class Point
{
public:
Point(float x, float y);
float getX();
float getY();
void setX(float x);
void setY(float y);
private:
// 存储一个值向量而不是单独的标量
std::vector<float> _coords;
};

#include "Point.h"
Point::Point(float x, float y)
{
_coords.push_back(x);
_coords.push_back(y);
}

float Point::getX() { return _coords[0]; }
float Point::getY() { return _coords[1]; }
void Point::setX(float x) { _coords[0] = x; }
void Point::setY(float y) { _coords[1] = y; }

#include "Point.h"
int main()
{
Point p1(1, 2);
Point p2 = Point(3, 4);
Point p3; // 编译错误，没有默认构造函数
std::cout << p1.getX() << "," << p1.getY() << "\n"
<< p2.getX() << "," << p2.getY() << "\n";
}

Point a(1,2), b(3,4), c(5,6);
// 使用访问器将 c 赋值为 a + b
c.setX(a.getX() + b.getX());
c.setY(a.getY() + b.getY());

c = a + b;

Point Point::operator+(const Point & p)
{
return Point(_x + p._x, _y + p._y);
}

Point & Point::operator=(const Point & p)
{
_x = p._x;
_y = p._y;
return *this;
}

#include "Point.h"

int main()
{
Point p1(0, 0), p2(1, 2), p3(3, 4);
p1 = p2 + p3;
std::cout << p1.getX() << "," << p1.getY() << "\n";
}

class Shape {
public:
Shape(int x=0, int y=0): _x(x), _y(y) {} // 构造函数
virtual ~Shape() {} // 析构函数
virtual float area() const = 0; // 纯虚函数
void printPosition() const; // 主体在别处出现

protected:
// 形状质心的坐标
int _x;
int _y;
};

#include "Shape.h"
class Rectangle: public Shape
{
public:
Rectangle(int width, int height, int x=0, int y=0) :
Shape(x,y),
_width(width),
_height(height)
{}

virtual float area() const override { return _width * _height; }

protected:
int _width;
int _height;
};

Rectangle rectangle(3, 4);
Shape & s_ref = rectangle;
Shape * s_ptr = &rectangle;

// mesh 是一个指向 Mesh 对象的指针
Mesh * mesh;

// params 是一个对 InputParameters 对象的引用
InputParameters & params;

// 以下是相同的
// value 是一个对常量 Real 对象的引用
const Real & value;
Real const & value;

// 创建几个形状
Rectangle r(3, 4);
Circle c(3, 10, 10);
printInformation(r); // 将 Rectangle 传入 Shape 引用
printInformation(c); // 将 Circle 传入 Shape 引用
...
void printInformation(const Shape & shape)
{
shape.printPosition();
std::cout << shape.area() << '\n';
}
// (0, 0)
// 12
// (10, 10)
// 28.274

// 类模板
template<typename T> class A { };

git clang-format branch_name_here

namespace moose // 小写命名空间名称
{
// 不要为命名空间增加缩进级别

int // 返回类型应单独一行
junkFunction()
{
// 缩进两个空格！
if (name == "moose") // 控制关键字 "if" 后有空格
{
// '&' 和 '*' 两侧有空格
SomeClass & a_ref;
SomeClass * a_pointer;
}

// 对于 if 后的单个语句，省略大括号，语句必须在单独一行
if (name == "squirrel")
doStuff();
else
doOtherStuff();

// 单个语句的分支和循环没有大括号
for (unsigned int i = 0; i < some_number; ++i) // 控制关键字 "for" 后有空格
doSomething();

// 赋值运算符周围有空格
Real foo = 5.0;

switch (stuff) // 控制关键字 "switch" 后有空格
{
// 缩进 case 语句
case 2:
junk = 4;
break;
case 3:
{ // 仅当 case 语句中有声明时才使用大括号
int bar = 9;
junk = bar;
break;
}
default:
junk = 8;
}

while (--foo) // 控制关键字 "while" 后有空格
std::cout << "Count down " << foo;
}

// （短）函数定义在单行上
SomeClass::SomeFunc() {}

// 构造函数初始化列表可以全部在同一行上。
SomeClass::SomeClass() : member_a(2), member_b(3) { }

// 长（即不适合单行）初始化列表使用四空格缩进和每项一行。
SomeOtherClass::SomeOtherClass()
: member_a(2),
member_b(3),
member_c(4),
member_d(5),
member_e(6),
member_f(7),
member_g(8),
member_h(9)
{ // 大括号在单独的行上，因为函数定义已经超过 1 行
}

} // namespace moose

auto dof = elem->dof_number(0, 0, 0); // dof_id_type, MOOSE 中默认为 64 位
int dof2 = elem->dof_number(0, 1, 0); // 对于 2B+ DoFs，dof2 会出错
auto & obj = getSomeObject();
auto & elem_it = mesh.active_local_elements_begin();
auto & [key, value] = map.find(some_item);

// 这里不能使用引用
for (auto it = container.begin(); it != container.end(); ++it)
doSomething(*it);

// 这里使用引用
for (auto & obj : container)
doSomething(obj);

// 尽可能显式地在捕获列表中列出捕获的变量。
std::for_each(container.begin(), container.end(), [= &local_var](Foo & foo) {
foo.item = local_var.item;
foo.item2 = local_var.item2;
});

unsigned int i = 4; // 内置类型
SomeObject junk(17); // 对象
SomeObject * stuff = functionReturningPointer(); // 指针
auto & [key, value] = functionReturningKVPair(); // 单个元组成员
auto heap_memory = std::make_unique<HeapObject>(a,b); // 非容器动态分配
std::vector<int> v = {1, 2, 4, 8, 16}; // 简单容器

find . -name '.[Chi]' -or -name '.py' | xargs perl -pli -e 's/\s+$//'

// 前向声明
class Something;

#include "LocalInclude.h"
#include "libmesh/libmesh_include.h"

#include <system_library>

/**

Kernel 类负责计算各种物理的残差。
/
class Kernel
{
public:
/*
这个构造函数应该最常使用。它初始化了残差计算所需的所有内部引用。
@param system 这个变量所在的系统
@param var_name 这个 Kernel 将为其计算残差的变量。
*/
Kernel(System * system, std::string var_name);
/**

这个函数用于根据 junk 获取 stuff。
@param junk 你想获取的 stuff 的索引
@return 与你传入的 junk 相关联的 stuff
*/
int returnStuff(int junk);
protected:
/// 这是这个类持有的 stuff。
std::vector<int> stuff;
};

class MyClass:
def init(self):
self.public_member
self._protected_member
self.__private_member

#pragma once

#include "BaseObject.h"

class CustomObject : public BaseObject
{
public:
static InputParameters validParams();

CustomObject(const InputParameters & parameters);

protected:

virtual Real doSomething() override;

const Real & _scale;
};

#include "CustomObject.h"

registerMooseObject("CustomApp", CustomObject);

InputParameters
CustomObject::validParams()
{
InputParameters params = BaseObject::validParams();
params.addClassDescription("The CustomObject does something with a scale parameter.");
params.addParam<Real>("scale", 1, "A scale factor for use when doing something.");
return params;
}

CustomObject::CustomObject(const InputParameters & parameters) :
BaseObject(parameters),
_scale(getParam<Real>("scale"))
{
}

double
CustomObject::doSomething()
{
// Do some sort of import calculation here that needs a scale factor
return _scale;
}

public:
...
static InputParameters Convection::validParams();

InputParameters
Convection::validParams()
{
InputParameters params = Kernel::validParams(); // 从父类开始
params.addRequiredParam<RealVectorValue>("velocity", "Velocity Vector");
params.addParam<Real>("coefficient", "Diffusion coefficient");
return params;
}

#include "MooseEnum.h"

// 有效选项在以空格分隔的列表中指定。
// 您可以选择提供默认值作为第二个参数。
// MooseEnums 区分大小写但大小写不敏感。
MooseEnum option_enum("first=1 second fourth=4", "second");

// 在字符串上下文中使用
if (option_enum == "first")
doSomething();

// 在整数上下文中使用
switch (option_enum)
{
case 1: ... break;
case 2: ... break;
case 4: ... break;
default: ... ;
}

InputParameters MyObject::validParams()
{
InputParameters params = ParentObject::validParams();
MooseEnum component("X Y Z"); // 未提供默认值
params.addRequiredParam<MooseEnum>("component", component,
"The X, Y, or Z component");
return params;
}

InputParameters MyObject::validParams()
{
InputParameters params = ParentObject::validParams();
MultiMooseEnum transforms("scale rotate translate");
params.addRequiredParam<MultiMooseEnum>("transforms", transforms,
"The transforms to perform");
return params;
}

#pragma once

#include "ADKernelGrad.h"

class ADDiffusion : public ADKernelGrad
{
public:
static InputParameters validParams();

ADDiffusion(const InputParameters & parameters);

protected:
virtual ADRealVectorValue precomputeQpResidual() override;
};

#include "ADDiffusion.h"

registerMooseObject("MooseApp", ADDiffusion);

InputParameters
ADDiffusion::validParams()
{
auto params = ADKernelGrad::validParams();
params.addClassDescription("Same as Diffusion in terms of physics/residual, but the Jacobian "
"is computed using forward automatic differentiation");
return params;
}

ADDiffusion::ADDiffusion(const InputParameters & parameters) : ADKernelGrad(parameters) {}

ADRealVectorValue
ADDiffusion::precomputeQpResidual()
{
return _grad_u[_qp];
}

#pragma once

// 在这里包含 "ADKernel" Kernel 以便我们可以扩展它
#include "ADKernel.h"

/**

计算残差贡献：K / mu * grad_u * grad_phi.
*/
class DarcyPressure : public ADKernel
{
public:
static InputParameters validParams();
DarcyPressure(const InputParameters & parameters);

protected:
/// ADKernel 对象必须覆盖 precomputeQpResidual
virtual ADReal computeQpResidual() override;

/// 从输入文件设置的引用
const Real & _permeability;
const Real & _viscosity;
};

#include "DarcyPressure.h"

registerMooseObject("DarcyThermoMechApp", DarcyPressure);

InputParameters
DarcyPressure::validParams()
{
InputParameters params = ADKernel::validParams();
params.addClassDescription("Compute the diffusion term for Darcy pressure (
p
p
) equation: "
"
−
n
a
b
l
a
c
d
o
t
f
r
a
c
m
a
t
h
b
f
K
m
u
n
a
b
l
a
p
=
0
−
nabla
cdot
fracmathbfKmu
nablap=0
");

// 添加一个必需的参数。如果输入文件中未提供此参数，MOOSE 将会报错。
params.addRequiredParam<Real>("permeability", "The permeability (
m
a
t
h
r
m
K
mathrmK
) of the fluid.");

// 添加一个带默认值的参数；此值可以在输入文件中覆盖。
params.addParam<Real>(
"viscosity",
7.98e-4,
"The viscosity (
m
u
mu
) of the fluid in Pa, the default is for water at 30 degrees C.");
return params;
}

DarcyPressure::DarcyPressure(const InputParameters & parameters)
: ADKernel(parameters),

Generated code
// 从输入文件中获取参数
_permeability(getParam<Real>("permeability")),
_viscosity(getParam<Real>("viscosity"))
content_copy
download
Use code with caution.
{
}

ADReal
DarcyPressure::computeQpResidual()
{
return (_permeability / _viscosity) * _grad_test[_i][_qp] * _grad_u[_qp];
}

@@ -1,51 +1,52 @@
[Mesh]
[gmg]

type = GeneratedMeshGenerator # Can generate simple lines, rectangles and rectangular prisms
dim = 2 # Dimension of the mesh
nx = 100 # Number of elements in the x direction
ny = 10 # Number of elements in the y direction
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 10
xmax = 0.304
ymax = 0.0257
[]
coord_type = RZ # Axisymmetric RZ
rz_coord_axis = X # Which axis the symmetry is around
coord_type = RZ
rz_coord_axis = X
[]
[Variables/pressure]

Adds a Linear Lagrange variable by default
[]

-[Kernels/diffusion]
-- type = ADDiffusion # Laplacian operator using automatic differentiation
-- variable = pressure # Operate on the "pressure" variable from above
+[Kernels/darcy_pressure]

type = DarcyPressure
variable = pressure
permeability = 0.8451e-9 # (m^2) 1mm spheres.
[]
[BCs]
[inlet]

type = DirichletBC # Simple u=value BC
variable = pressure # Variable to be set
boundary = left # Name of a sideset in the mesh
value = 4000 # (Pa) From Figure 2 from paper. First data point for 1mm spheres.
type = DirichletBC
variable = pressure
boundary = left
value = 4000
[]
[outlet]
type = DirichletBC<<<{"description": "Imposes the essential boundary condition
u
=
g
u=g
, where
g
g
 is a constant, controllable value.", "href": "../source/bcs/DirichletBC.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
value<<<{"description": "Value of the BC"}>>> = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
type = DirichletBC
variable = pressure
boundary = right
value = 0
[]
[]
[Problem]

type = FEProblem # This is the "normal" type of Finite Element Problem in MOOSE
type = FEProblem
[]
[Executioner]

type = Steady # Steady state problem
solve_type = NEWTON # Perform a Newton solve, uses AD to compute Jacobian terms
petsc_options_iname = '-pc_type -pc_hypre_type' # PETSc option pairs with values below
type = Steady
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
[]
[Outputs]

exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true # Output Exodus format
exodus = true
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step02_darcy_pressure
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step2.i

[Mesh]
[fmg]
type = FileMeshGenerator
file = square.e
[]
[]

[Mesh]
[generated]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = 1
nx = 2
ymin = -2
ymax = 3
ny = 3
elem_type = 'TRI3'
[]
[]

[Mesh]
file = three_block.e

这些名称将动态应用于网格
以便它们可以在输入文件中使用
此外，它们将显示在输出文件中
block_id = '1 2 3'
block_name = 'wood steel copper'

boundary_id = '1 2'
boundary_name = 'left right'
[]

[BCs]
active = 'left right'

[left]
type = DirichletBC
variable = u
boundary = 'left'
value = 0
[]

[right]
type = DirichletBC
variable = u
boundary = 'right'
value = 1
[]
[]

[Materials]
active = empty

[empty]
type = MTMaterial
block = 'wood steel copper'
[]
[]

[Outputs]
exodus = true
[]

[Outputs]
[out]
type = Exodus
[]
[]

[Outputs]
[out]
type = Exodus
output_material_properties = true
# 从输出中移除一些量
hide = 'power_pp pressure_var'
[]
[]

[Outputs]
interval = 10 # 这是一个时间步间隔
[exo]
type = Exodus
interval = 1 # 覆盖顶层的 interval
[]
[cp]
type = Checkpoint # 使用顶层指定的 interval
[]
[]

[Outputs]
exodus = true # 创建 input_out.e
[other] # 创建 input_other.e
type = Exodus
interval = 2
[]
[base]
type = Exodus
file_base = out # 创建 out.e
[]
[]

addParam<MaterialPropertyName>("combination_property_name", 12345,
"The name of the property providing the luggage combination");

[Materials]
[block_1]
type = OutputTestMaterial
block = 1
output_properties = 'real_property tensor_property'
outputs = exodus
variable = u
[]
[block_2]
type = OutputTestMaterial
block = 2
output_properties = 'vector_property tensor_property'
outputs = exodus
variable = u
[]
[]

[Outputs]
exodus = true
[]

[Functions]
[sin_fn]
type = ParsedFunction
expression = sin(x)
[]
[cos_fn]
type = ParsedFunction
expression = cos(x)
[]
[fn]
type = ParsedFunction
expression = 's/c'
symbol_names = 's c'
symbol_values = 'sin_fn cos_fn'
[]
[]

// 添加一个带默认常量的 Function
params.addParam<FunctionName>("pressure_grad", "0.5", "The gradient of ...");

// 添加一个带默认解析函数的 Function
params.addParam<FunctionName>("power_history", "t+100*sin(y)", "The power history of ...");

#pragma once

#include "Material.h"

/**

Material 对象继承自 Material 并重写 computeQpProperties。
它们的工作是声明供其他对象（如 Kernels 和 BoundaryConditions）在计算中使用的属性。
*/
class PackedColumn : public Material
{
public:
static InputParameters validParams();
PackedColumn(const InputParameters & parameters);

protected:
/// 必要的重写。这里是计算属性值的地方。
virtual void computeQpProperties() override;

/// 柱中球体的半径
const Function & _input_radius;

/// 来自输入文件的粘度值
const Real & _input_viscosity;

/// 根据半径 (mm) 计算的渗透率 (K)
ADMaterialProperty<Real> & _permeability;

/// 流体的粘度 (mu)
ADMaterialProperty<Real> & _viscosity;
};

#include "PackedColumn.h"
#include "Function.h"

registerMooseObject("DarcyThermoMechApp", PackedColumn);

InputParameters
PackedColumn::validParams()
{
InputParameters params = Material::validParams();

// 用于插值渗透率的球体半径参数。
params.addParam<FunctionName>("radius",
"1.0",
"The radius of the steel spheres (mm) that are packed in the "
"column for computing permeability.");
params.addParam<Real>(
"viscosity",
7.98e-4,
"The viscosity (
m
u
mu
) of the fluid in Pa, the default is for water at 30 degrees C.");

return params;
}

PackedColumn::PackedColumn(const InputParameters & parameters)
: Material(parameters),

Generated code
// 从输入文件中获取参数
_input_radius(getFunction("radius")),
_input_viscosity(getParam<Real>("viscosity")),

// 此对象生产的材料属性
_permeability(declareADProperty<Real>("permeability")),
_viscosity(declareADProperty<Real>("viscosity"))
content_copy
download
Use code with caution.
{
}

void
PackedColumn::computeQpProperties()
{
// 来自论文：表 1
std::vector<Real> sphere_sizes = {1, 3};
std::vector<Real> permeability = {0.8451e-9, 8.968e-9};

Real value = _input_radius.value(_t, _q_point[_qp]);
mooseAssert(value >= 1 && value <= 3,
"The radius range must be in the range [1, 3], but " << value << " provided.");

_viscosity[_qp] = _input_viscosity;

// 我们将使用上面两点的简单线性插值来计算渗透率：
// y0 * (x1 - x) + y1 * (x - x0)
// y(x) = -------------------------------
// x1 - x0
_permeability[_qp] =
(permeability[0] * (sphere_sizes[1] - value) + permeability[1] * (value - sphere_sizes[0])) /
(sphere_sizes[1] - sphere_sizes[0]);
}

@@ -1,24 +1,32 @@

#pragma once

// 在这里包含 "ADKernel" Kernel 以便我们可以扩展它
#include "ADKernel.h"

/**

计算残差贡献：K / mu * grad_u * grad_phi.
*/
class DarcyPressure : public ADKernel
{
public:
static InputParameters validParams();
DarcyPressure(const InputParameters & parameters);

protected:
/// ADKernel 对象必须覆盖 precomputeQpResidual
virtual ADReal computeQpResidual() override;

/// 从输入文件设置的引用
const Real & _permeability;
const Real & _viscosity;
// 从 Material 系统设置的引用
/// 粘度。这被声明为 \p ADMaterialProperty，意味着任何来自生产 \p Material 的导数信息
/// 都将被保留，非线性求解的完整性也将同样被保留
const ADMaterialProperty<Real> & _permeability;
/// 粘度。这被声明为 \p ADMaterialProperty，意味着任何来自生产 \p Material 的导数信息
/// 都将被保留，非线性求解的完整性也将同样被保留
const ADMaterialProperty<Real> & _viscosity;
};

@@ -1,37 +1,26 @@

#include "DarcyPressure.h"

registerMooseObject("DarcyThermoMechApp", DarcyPressure);

InputParameters
DarcyPressure::validParams()
{
InputParameters params = ADKernel::validParams();
params.addClassDescription("Compute the diffusion term for Darcy pressure (
p
p
) equation: "
"
−
n
a
b
l
a
c
d
o
t
f
r
a
c
m
a
t
h
b
f
K
m
u
n
a
b
l
a
p
=
0
−
nabla
cdot
fracmathbfKmu
nablap=0
");
// 添加一个必需的参数。如果输入文件中未提供此参数，MOOSE 将会报错。
params.addRequiredParam<Real>("permeability", "The permeability (
m
a
t
h
r
m
K
mathrmK
) of the fluid.");
// 添加一个带默认值的参数；此值可以在输入文件中覆盖。
params.addParam<Real>(
Generated code
"viscosity",
content_copy
download
Use code with caution.
Generated code
7.98e-4,
content_copy
download
Use code with caution.
Generated code
"The viscosity ($\\mu$) of the fluid in Pa, the default is for water at 30 degrees C.");
content_copy
download
Use code with caution.
return params;
}
DarcyPressure::DarcyPressure(const InputParameters & parameters)
: ADKernel(parameters),
// 从输入文件中获取参数
_permeability(getParam<Real>("permeability")),
_viscosity(getParam<Real>("viscosity"))
_permeability(getADMaterialProperty<Real>("permeability")),
_viscosity(getADMaterialProperty<Real>("viscosity"))
{
}
ADReal
DarcyPressure::computeQpResidual()
{

return (_permeability / _viscosity) * _grad_test[_i][_qp] * _grad_u[_qp];
return (_permeability[_qp] / _viscosity[_qp]) * _grad_test[_i][_qp] * _grad_u[_qp];
}

@@ -1,52 +1,55 @@
[Mesh]
[gmg]

type = GeneratedMeshGenerator # Can generate simple lines, rectangles and rectangular prisms
dim = 2 # Dimension of the mesh
nx = 100 # Number of elements in the x direction
ny = 10 # Number of elements in the y direction
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 10
xmax = 0.304
ymax = 0.0257
[]
coord_type = RZ # Axisymmetric RZ
rz_coord_axis = X # Which axis the symmetry is around
coord_type = RZ
rz_coord_axis = X
[]
[Variables/pressure]

Adds a Linear Lagrange variable by default
[]

[Kernels/darcy_pressure]
type = DarcyPressure
variable = pressure
-- permeability = 0.8451e-9 # (m^2) 1mm spheres.
[]

[BCs]
[inlet]

type = DirichletBC # Simple u=value BC
variable = pressure # Variable to be set
boundary = left # Name of a sideset in the mesh
value = 4000 # (Pa) From Figure 2 from paper. First data point for 1mm spheres.
type = DirichletBC
variable = pressure
boundary = left
value = 4000
[]
[outlet]
type = DirichletBC<<<{"description": "Imposes the essential boundary condition
u
=
g
u=g
, where
g
g
 is a constant, controllable value.", "href": "../source/bcs/DirichletBC.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
value<<<{"description": "Value of the BC"}>>> = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
type = DirichletBC
variable = pressure
boundary = right
value = 0
[]
[]
+[Materials/column]

type = PackedColumn
+[]
[Problem]

type = FEProblem # This is the "normal" type of Finite Element Problem in MOOSE
type = FEProblem
[]
[Executioner]

type = Steady # Steady state problem
solve_type = NEWTON # Perform a Newton solve, uses AD to compute Jacobian terms
petsc_options_iname = '-pc_type -pc_hypre_type' # PETSc option pairs with values below
type = Steady
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
[]
[Outputs]

exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true # Output Exodus format
exodus = true
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step03_darcy_material
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step3.i

@@ -1,55 +1,57 @@
[Mesh]
[gmg]

type = GeneratedMeshGenerator # Can generate simple lines, rectangles and rectangular prisms
dim = 2 # Dimension of the mesh
nx = 100 # Number of elements in the x direction
ny = 10 # Number of elements in the y direction
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 10
xmax = 0.304
ymax = 0.0257
[]
coord_type = RZ # Axisymmetric RZ
rz_coord_axis = X # Which axis the symmetry is around
coord_type = RZ
rz_coord_axis = X
[]
[Variables/pressure]

Adds a Linear Lagrange variable by default
[]

[Kernels/darcy_pressure]
type = DarcyPressure
variable = pressure
[]

[BCs]
[inlet]

type = DirichletBC # Simple u=value BC
variable = pressure # Variable to be set
boundary = left # Name of a sideset in the mesh
value = 4000 # (Pa) From Figure 2 from paper. First data point for 1mm spheres.
type = DirichletBC
variable = pressure
boundary = left
value = 4000
[]
[outlet]
type = DirichletBC<<<{"description": "Imposes the essential boundary condition
u
=
g
u=g
, where
g
g
 is a constant, controllable value.", "href": "../source/bcs/DirichletBC.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
value<<<{"description": "Value of the BC"}>>> = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
type = DirichletBC
variable = pressure
boundary = right
value = 0
[]
[]
[Materials/column]
type = PackedColumn
-+ radius = '1 + 2/3.04x'
++ radius = '1 + 2/0.304x'

outputs = exodus
[]
[Problem]

type = FEProblem # This is the "normal" type of Finite Element Problem in MOOSE
type = FEProblem
[]
[Executioner]

type = Steady # Steady state problem
solve_type = NEWTON # Perform a Newton solve, uses AD to compute Jacobian terms
petsc_options_iname = '-pc_type -pc_hypre_type' # PETSc option pairs with values below
type = Steady
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
[]
[Outputs]

exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true # Output Exodus format
exodus = true
[]

tests/
kernels/
my_kernel_test/
my_kernel_test.e [输入网格]
my_kernel_test.i [输入文件]
tests [测试规范文件]
gold/ [用于验证解的黄金标准文件夹]
out.e [解]

[Tests]
[my_kernel_test]
type = Exodiff
input = my_kernel_test.i
exodiff = my_kernel_test_out.e
[]

[kernel_check_exception]
type = RunException
input = my_kernel_exception.i
expect_err = 'Bad stuff happened with variable \w+'
[]
[]

./run_tests -j 12

./run_tests -h

./run_tests --dump

[AuxVariables]
[aux]
order = FIRST
family = LAGRANGE_VEC
[]
[]
[AuxKernels]
[parsed]
type = ParsedVectorAux
variable = aux
expression_x = 'x + y'
expression_y = 'x - y'
[]
[]

#pragma once

#include "AuxKernel.h"

/**

辅助内核，负责在给定多种流体属性和压力梯度的情况下计算达西速度。
*/
class DarcyVelocity : public VectorAuxKernel
{
public:
static InputParameters validParams();
DarcyVelocity(const InputParameters & parameters);

protected:
/**

AuxKernels 必须重写 computeValue。computeValue() 在每个求积点上被调用。
对于节点辅助变量，这些求积点与节点重合。
*/
virtual RealVectorValue computeValue() override;
/// 耦合变量的梯度
const VariableGradient & _pressure_gradient;

/// 保存来自材料系统的渗透率和粘度
const ADMaterialProperty<Real> & _permeability;
const ADMaterialProperty<Real> & _viscosity;
};

#include "DarcyVelocity.h"

#include "metaphysicl/raw_type.h"

registerMooseObject("DarcyThermoMechApp", DarcyVelocity);

InputParameters
DarcyVelocity::validParams()
{
InputParameters params = VectorAuxKernel::validParams();

// 添加一个“耦合参数”以从输入文件中获取一个变量。
params.addRequiredCoupledVar("pressure", "The pressure field.");

return params;
}

DarcyVelocity::DarcyVelocity(const InputParameters & parameters)
: VectorAuxKernel(parameters),

Generated code
// 获取变量的梯度
_pressure_gradient(coupledGradient("pressure")),

// 设置对渗透率 MaterialProperty 的引用。
// 只有作用于单元辅助变量的 AuxKernels 才能这样做
_permeability(getADMaterialProperty<Real>("permeability")),

// 设置对粘度 MaterialProperty 的引用。
// 只有作用于单元辅助变量的 AuxKernels 才能这样做
_viscosity(getADMaterialProperty<Real>("viscosity"))
content_copy
download
Use code with caution.
{
}

RealVectorValue
DarcyVelocity::computeValue()
{
// 访问此求积点处的压力梯度，然后提取其请求的“分量”（x、y 或 z）。
// 注意，获取梯度的特定分量是使用括号运算符完成的。
return -MetaPhysicL::raw_value(_permeability[_qp] / _viscosity[_qp]) * _pressure_gradient[_qp];
}

@@ -1,55 +1,70 @@
[Mesh]
[gmg]

type = GeneratedMeshGenerator # Can generate simple lines, rectangles and rectangular prisms
dim = 2 # Dimension of the mesh
nx = 100 # Number of elements in the x direction
ny = 10 # Number of elements in the y direction
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 10
xmax = 0.304
ymax = 0.0257
[]
coord_type = RZ # Axisymmetric RZ
rz_coord_axis = X # Which axis the symmetry is around
coord_type = RZ
rz_coord_axis = X
[]
[Variables/pressure]

Adds a Linear Lagrange variable by default
[]

+[AuxVariables/velocity]

order = CONSTANT
family = MONOMIAL_VEC
+[]
[Kernels/darcy_pressure]
type = DarcyPressure
variable = pressure
+[]
+
+[AuxKernels]

[velocity]
type = DarcyVelocity
variable = velocity
execute_on = timestep_end
pressure = pressure
[]
[]
[BCs]
[inlet]

type = DirichletBC # Simple u=value BC
variable = pressure # Variable to be set
boundary = left # Name of a sideset in the mesh
value = 4000 # (Pa) From Figure 2 from paper. First data point for 1mm spheres.
type = DirichletBC
variable = pressure
boundary = left
value = 4000
[]
[outlet]
type = DirichletBC<<<{"description": "Imposes the essential boundary condition
u
=
g
u=g
, where
g
g
 is a constant, controllable value.", "href": "../source/bcs/DirichletBC.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
value<<<{"description": "Value of the BC"}>>> = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
type = DirichletBC
variable = pressure
boundary = right
value = 0
[]
[]
[Materials/column]
type = PackedColumn

radius = 1
[]
[Problem]

type = FEProblem # This is the "normal" type of Finite Element Problem in MOOSE
type = FEProblem
[]
[Executioner]

type = Steady # Steady state problem
solve_type = NEWTON # Perform a Newton solve, uses AD to compute Jacobian terms
petsc_options_iname = '-pc_type -pc_hypre_type' # PETSc option pairs with values below
type = Steady
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
[]
[Outputs]

exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true # Output Exodus format
exodus = true
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step04_velocity_aux
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step4.i

cd ~/projects/moose/tutorials/darcy_thermo_mech/step04_velocity_aux
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step4.i Executioner/nl_rel_tol=1e-12

#pragma once

#include "ADDiffusion.h"

class ADHeatConduction : public ADDiffusion
{
public:
static InputParameters validParams();

ADHeatConduction(const InputParameters & parameters);

protected:
virtual ADRealVectorValue precomputeQpResidual() override;

const ADMaterialProperty<Real> & _thermal_conductivity;
};

#include "ADHeatConduction.h"

registerMooseObject("HeatTransferApp", ADHeatConduction);

InputParameters
ADHeatConduction::validParams()
{
InputParameters params = ADDiffusion::validParams();
params.addParam<MaterialPropertyName>("thermal_conductivity",
"thermal_conductivity",
"the name of the thermal conductivity material property");
params.set<bool>("use_displaced_mesh") = true;
return params;
}

ADHeatConduction::ADHeatConduction(const InputParameters & parameters)
: ADDiffusion(parameters),
_thermal_conductivity(getADMaterialProperty<Real>("thermal_conductivity"))
{
}

ADRealVectorValue
ADHeatConduction::precomputeQpResidual()
{
return _thermal_conductivity[_qp] * ADDiffusion::precomputeQpResidual();
}

[Materials/material]
type = GenericConstantMaterial
prop_names = 'conductivity density'
prop_values = '0.01 200'
[]

[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 10
xmax = 0.304 # 测试室长度
ymax = 0.0257 # 测试室半径
[]
coord_type = RZ
rz_coord_axis = X
[]

[Variables]
[temperature]
[]
[]

[Kernels]
[heat_conduction]
type = ADHeatConduction
variable = temperature
[]
[]

[BCs]
[inlet_temperature]
type = DirichletBC
variable = temperature
boundary = left
value = 350 # (K)
[]
[outlet_temperature]
type = DirichletBC
variable = temperature
boundary = right
value = 300 # (K)
[]
[]

[Materials/steel]
type = ADGenericConstantMaterial
prop_names = thermal_conductivity
prop_values = 18 # K: (W/m*K) 来自维基百科 @296K
[]

[Problem]
type = FEProblem
[]

[Executioner]
type = Steady
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
exodus = true
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step5_heat_conduction
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step5a_steady.i

[Executioner]
type = Steady
[]

[Executioner]
type = Transient
scheme = 'implicit-euler'
solve_type = 'PJFNK'
start_time = 0.0
num_steps = 5
dt = 0.1
[]

[Executioner]
type = Transient
start_time = 0.0
num_steps = 6
dt = 0.1
[TimeIntegrator]
type = NewmarkBeta
beta = 0.4225
gamma = 0.8
[]
[]

[Executioner]
type = Transient
solve_type = NEWTON
start_time = 0.0
dtmin = 1.0
end_time = 10.0
[TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 1
linear_iteration_ratio = 1
dt = 5.0
[]
[]

[TimeSteppers]
[constant]
type = ConstantDT
dt = 0.2
[]
[hit_these_times]
type = TimeSequenceStepper
time_sequence = '0.5 1 1.5 2.1'
[]
[]

#pragma once

#include "ADTimeDerivative.h"

class ADHeatConductionTimeDerivative : public ADTimeDerivative
{
public:
static InputParameters validParams();

ADHeatConductionTimeDerivative(const InputParameters & parameters);

protected:
virtual ADReal precomputeQpResidual() override;

/// Specific heat material property
const ADMaterialProperty<Real> & _specific_heat;

/// Density material property
const ADMaterialProperty<Real> & _density;
};

#include "ADHeatConductionTimeDerivative.h"

registerMooseObject("HeatTransferApp", ADHeatConductionTimeDerivative);

InputParameters
ADHeatConductionTimeDerivative::validParams()
{
InputParameters params = ADTimeDerivative::validParams();
params.addClassDescription(
"AD Time derivative term
r
h
o
c
p
f
r
a
c
p
a
r
t
i
a
l
T
p
a
r
t
i
a
l
t
rhoc
p
​

fracpartialTpartialt
 of "
"the heat equation for quasi-constant specific heat
c
p
c
p
​

 and the density
r
h
o
rho
.");
params.set<bool>("use_displaced_mesh") = true;
params.addParam<MaterialPropertyName>(
"specific_heat", "specific_heat", "Property name of the specific heat material property");
params.addParam<MaterialPropertyName>(
"density_name", "density", "Property name of the density material property");
return params;
}

ADHeatConductionTimeDerivative::ADHeatConductionTimeDerivative(const InputParameters & parameters)
: ADTimeDerivative(parameters),
_specific_heat(getADMaterialProperty<Real>("specific_heat")),
_density(getADMaterialProperty<Real>("density_name"))
{
}

ADReal
ADHeatConductionTimeDerivative::precomputeQpResidual()
{
return _specific_heat[_qp] * _density[_qp] * ADTimeDerivative::precomputeQpResidual();
}

@@ -1,60 +1,68 @@
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 10
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
[]
coord_type = RZ
rz_coord_axis = X
[]

[Variables]
[temperature]

initial_condition = 300 # Start at room temperature
[]
[]
[Kernels]
[heat_conduction]

type = ADHeatConduction
variable = temperature
type = ADHeatConduction
variable = temperature
[]
[heat_conduction_time_derivative]
type = ADHeatConductionTimeDerivative
variable = temperature
specific_heat = specific_heat
density_name = density
[]
[]
[BCs]
[inlet_temperature]

type = DirichletBC
variable = temperature
boundary = left
value = 350 # (K)
type = DirichletBC
variable = temperature
boundary = left
value = 350
[]
[outlet_temperature]
type = DirichletBC
variable = temperature
boundary = right
value = 300 # (K)
type = DirichletBC
variable = temperature
boundary = right
value = 300
[]
[]
[Materials/steel]

type = ADGenericConstantMaterial
-- prop_names = thermal_conductivity
-- prop_values = 18 # K: (W/m*K) from wikipedia @296K
type = ADGenericConstantMaterial
prop_names = 'thermal_conductivity specific_heat density'
prop_values = '18 0.466 8000' # W/m*K, J/kg-K, kg/m^3 @ 296K
[]
[Problem]
type = FEProblem
[]

[Executioner]
-- type = Steady

type = Transient
num_steps = 10
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
[]
[Outputs]

exodus = true
exodus = true
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step5_heat_conduction
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step5b_transient.i

#pragma once

#include "DirichletBCBase.h"

/**

Boundary condition of a Dirichlet type
Sets the value in the node
*/
class DirichletBC : public DirichletBCBase
{
public:
static InputParameters validParams();
DirichletBC(const InputParameters & parameters);

protected:
virtual Real computeQpValue() override;

/// The value for this BC
const Real & _value;
};

#include "DirichletBC.h"

registerMooseObject("MooseApp", DirichletBC);

InputParameters
DirichletBC::validParams()
{
InputParameters params = DirichletBCBase::validParams();
params.addRequiredParam<Real>("value", "Value of the BC");
params.declareControllable("value");
params.addClassDescription("Imposes the essential boundary condition
u
=
g
u=g
, where
g
g
 "
"is a constant, controllable value.");
return params;
}

DirichletBC::DirichletBC(const InputParameters & parameters)
: DirichletBCBase(parameters), _value(getParam<Real>("value"))
{
}

Real
DirichletBC::computeQpValue()
{
return _value;
}

#pragma once

#include "GenericIntegratedBC.h"

/**

Implements a simple constant Neumann BC where grad(u)=value on the boundary.
Uses the term produced from integrating the diffusion operator by parts.
*/
template <bool is_ad>
class NeumannBCTempl : public GenericIntegratedBC<is_ad>
{
public:
static InputParameters validParams();
NeumannBCTempl(const InputParameters & parameters);

protected:
virtual GenericReal<is_ad> computeQpResidual() override;

/// Value of grad(u) on the boundary.
const Real & _value;

usingGenericIntegratedBCMembers;
};

typedef NeumannBCTempl<false> NeumannBC;
typedef NeumannBCTempl<true> ADNeumannBC;

#include "NeumannBC.h"

registerMooseObject("MooseApp", NeumannBC);
registerMooseObject("MooseApp", ADNeumannBC);

template <bool is_ad>
InputParameters
NeumannBCTempl<is_ad>::validParams()
{
InputParameters params = GenericIntegratedBC<is_ad>::validParams();
params.addParam<Real>("value",
0.0,
"For a Laplacian problem, the value of the gradient dotted with the "
"normals on the boundary.");
params.declareControllable("value");
params.addClassDescription("Imposes the integrated boundary condition "
"
f
r
a
c
p
a
r
t
i
a
l
u
p
a
r
t
i
a
l
n
=
h
fracpartialupartialn=h
, "
"where
h
h
 is a constant, controllable value.");
return params;
}

template <bool is_ad>
NeumannBCTempl<is_ad>::NeumannBCTempl(const InputParameters & parameters)
: GenericIntegratedBC<is_ad>(parameters), _value(this->template getParam<Real>("value"))
{
}

template <bool is_ad>
GenericReal<is_ad>
NeumannBCTempl<is_ad>::computeQpResidual()
{
return -_test[_i][_qp] * _value;
}

template class NeumannBCTempl<false>;
template class NeumannBCTempl<true>;

#pragma once

// 包含基类以便可以扩展它
#include "ADIntegratedBC.h"

/**

表示热传导“无 BC”边界条件的 IntegratedBC。
残差就是 -testkgrad_u*normal... 这是通过分部积分得到的项。
这是求解对流/扩散方程时截断较长域的标准技术。
另见：Griffiths, David F. "The 'no boundary condition' outflow
boundary condition.", International Journal for Numerical Methods
in Fluids, vol. 24, no. 4, 1997, pp. 393-411.
*/
class HeatConductionOutflow : public ADIntegratedBC
{
public:
static InputParameters validParams();
HeatConductionOutflow(const InputParameters & parameters);

protected:
/**

调用此函数以在边界上积分残差。
*/
virtual ADReal computeQpResidual() override;
/// 材料的导热系数
const ADMaterialProperty<Real> & _thermal_conductivity;
};

#include "HeatConductionOutflow.h"

registerMooseObject("DarcyThermoMechApp", HeatConductionOutflow);

InputParameters
HeatConductionOutflow::validParams()
{
InputParameters params = ADIntegratedBC::validParams();
params.addClassDescription("Compute the outflow boundary condition.");
return params;
}

HeatConductionOutflow::HeatConductionOutflow(const InputParameters & parameters)
: ADIntegratedBC(parameters),
_thermal_conductivity(getADMaterialProperty<Real>("thermal_conductivity"))
{
}

ADReal
HeatConductionOutflow::computeQpResidual()
{
return -_test[_i][_qp] * _thermal_conductivity[_qp] * _grad_u[_qp] * _normals[_qp];
}

@@ -1,68 +1,65 @@
[Mesh]
[gmg]

type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 10
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
type = GeneratedMeshGenerator
dim = 2
nx = 100
ny = 10
xmax = 0.304
ymax = 0.0257
[]
coord_type = RZ
rz_coord_axis = X
[]
[Variables]
[temperature]

initial_condition = 300 # Start at room temperature
initial_condition = 300
[]
[]
[Kernels]
[heat_conduction]

type = ADHeatConduction
variable = temperature
type = ADHeatConduction
variable = temperature
[]
[heat_conduction_time_derivative]
type = ADHeatConductionTimeDerivative
variable = temperature
-- specific_heat = specific_heat
-- density_name = density
type = ADHeatConductionTimeDerivative
variable = temperature
[]
[]
[BCs]
[inlet_temperature]

type = DirichletBC
variable = temperature
boundary = left
value = 350 # (K)
type = DirichletBC
variable = temperature
boundary = left
value = 350
[]
[outlet_temperature]
-- type = DirichletBC
-- variable = temperature
-- boundary = right
-- value = 300 # (K)
type = HeatConductionOutflow
variable = temperature
boundary = right
[]
[]
[Materials/steel]

type = ADGenericConstantMaterial
prop_names = 'thermal_conductivity specific_heat density'
prop_values = '18 0.466 8000' # W/m*K, J/kg-K, kg/m^3 @ 296K
type = ADGenericConstantMaterial
prop_names = 'thermal_conductivity specific_heat density'
prop_values = '18 0.466 8000'
[]
[Problem]
type = FEProblem
[]

[Executioner]
type = Transient
num_steps = 10
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
[]

[Outputs]

exodus = true
exodus = true
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step05_heat_conduction
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step5b_transient.i

#pragma once

#include "ADKernelValue.h"

/**

Kernel which implements the convective term in the transient heat
conduction equation, and provides coupling with the Darcy pressure
equation.
*/
class DarcyAdvection : public ADKernelValue
{
public:
static InputParameters validParams();
DarcyAdvection(const InputParameters & parameters);

protected:
/// ADKernelValue objects must override precomputeQpResidual
virtual ADReal precomputeQpResidual() override;

/// The gradient of pressure
const ADVariableGradient & _pressure_grad;

/// These references will be set by the initialization list so that
/// values can be pulled from the Material system.
const ADMaterialProperty<Real> & _permeability;
const ADMaterialProperty<Real> & _porosity;
const ADMaterialProperty<Real> & _viscosity;
const ADMaterialProperty<Real> & _density;
const ADMaterialProperty<Real> & _specific_heat;
};

#include "DarcyAdvection.h"

registerMooseObject("DarcyThermoMechApp", DarcyAdvection);

InputParameters
DarcyAdvection::validParams()
{
InputParameters params = ADKernelValue::validParams();
params.addRequiredCoupledVar("pressure", "The variable representing the pressure.");
return params;
}

DarcyAdvection::DarcyAdvection(const InputParameters & parameters)
: ADKernelValue(parameters),
// Couple to the gradient of the pressure
_pressure_grad(adCoupledGradient("pressure")),
// Grab necessary material properties
_permeability(getADMaterialProperty<Real>("permeability")),
_porosity(getADMaterialProperty<Real>("porosity")),
_viscosity(getADMaterialProperty<Real>("viscosity")),
_density(getADMaterialProperty<Real>("density")),
_specific_heat(getADMaterialProperty<Real>("specific_heat"))
{
}

ADReal
DarcyAdvection::precomputeQpResidual()
{
// See also: E. Majchrzak and L. Turchan, "The Finite Difference
// Method For Transient Convection Diffusion", Scientific Research
// of the Institute of Mathematics and Computer Science, vol. 1,
// no. 11, 2012, pp. 63-72.
// http://srimcs.im.pcz.pl/2012_1/art_07.pdf

// http://en.wikipedia.org/wiki/Superficial_velocity
ADRealVectorValue superficial_velocity =
_porosity[_qp] * -(_permeability[_qp] / _viscosity[_qp]) * _pressure_grad[_qp];

return _density[_qp] * _specific_heat[_qp] * superficial_velocity * _grad_u[_qp];
}

_permeability[_qp] =
(permeability[0] * (sphere_sizes[1] - value) + permeability[1] * (value - sphere_sizes[0])) /
(sphere_sizes[1] - sphere_sizes[0]);

_permeability_interpolation.setData(sphere_sizes, permeability);

_permeability[_qp] = _permeability_interpolation.sample(value);

@@ -1,34 +1,95 @@

#pragma once

#include "Material.h"
+
+// A helper class from MOOSE that linear interpolates x,y data
+#include "LinearInterpolation.h"

/**

Material objects inherit from Material and override computeQpProperties.
Their job is to declare properties for use by other objects in the
calculation such as Kernels and BoundaryConditions.
*/
class PackedColumn : public Material
{
public:
static InputParameters validParams();
PackedColumn(const InputParameters & parameters);

protected:
/// Necessary override. This is where the values of the properties are computed.
virtual void computeQpProperties() override;

/**
Helper function for reading CSV data for use in an interpolator object.
*/
bool initInputData(const std::string & param_name, ADLinearInterpolation & interp);
/// The radius of the spheres in the column
const Function & _input_radius;
/// Value of viscosity from the input file
const Real & _input_viscosity;
/// The input porosity
const Function & _input_porosity;
/// The permeability (K) computed based on the radius (mm)
/// Temperature
const ADVariableValue & _temperature;
/// Compute permeability based on the radius (mm)
LinearInterpolation _permeability_interpolation;
/// Fluid viscosity
bool _use_fluid_mu_interp;
const Real & _fluid_mu;
ADLinearInterpolation _fluid_mu_interpolation;
/// Fluid thermal conductivity
bool _use_fluid_k_interp = false;
const Real & _fluid_k;
ADLinearInterpolation _fluid_k_interpolation;
/// Fluid density
bool _use_fluid_rho_interp = false;
const Real & _fluid_rho;
ADLinearInterpolation _fluid_rho_interpolation;
/// Fluid specific heat
bool _use_fluid_cp_interp;
const Real & _fluid_cp;
ADLinearInterpolation _fluid_cp_interpolation;
/// Solid thermal conductivity
bool _use_solid_k_interp = false;
const Real & _solid_k;
ADLinearInterpolation _solid_k_interpolation;
/// Solid density
bool _use_solid_rho_interp = false;
const Real & _solid_rho;
ADLinearInterpolation _solid_rho_interpolation;
/// Solid specific heat
bool _use_solid_cp_interp;
const Real & _solid_cp;
ADLinearInterpolation _solid_cp_interpolation;
/// The permeability (K)
ADMaterialProperty<Real> & _permeability;
/// The viscosity of the fluid (mu)
ADMaterialProperty<Real> & _viscosity;
/// The porosity (eps)
ADMaterialProperty<Real> & _porosity;
/// The bulk thermal conductivity
ADMaterialProperty<Real> & _thermal_conductivity;
/// The bulk heat capacity
ADMaterialProperty<Real> & _specific_heat;
/// The bulk density
ADMaterialProperty<Real> & _density;
};

@@ -1,58 +1,177 @@

#include "PackedColumn.h"
#include "Function.h"
+#include "DelimitedFileReader.h"

registerMooseObject("DarcyThermoMechApp", PackedColumn);

InputParameters
PackedColumn::validParams()
{
InputParameters params = Material::validParams();

params.addRequiredCoupledVar("temperature", "The temperature (C) of the fluid.");
// Parameter for radius of the spheres used to interpolate permeability.
params.addParam<FunctionName>("radius",
"1.0",
"The radius of the steel spheres (mm) that are packed in the "
"column for computing permeability.");
// http://en.wikipedia.org/wiki/Close-packing_of_equal_spheres
params.addParam<FunctionName>(
Generated code
"porosity", 0.25952, "Porosity of porous media, default is for closed packed spheres.");
content_copy
download
Use code with caution.
// Fluid properties
params.addParam<Real>(
Generated code
"viscosity",
content_copy
download
Use code with caution.
Generated code
7.98e-4,
content_copy
download
Use code with caution.
Generated code
"The viscosity ($\\mu$) of the fluid in Pa, the default is for water at 30 degrees C.");
content_copy
download
Use code with caution.
Generated code
"fluid_viscosity", 1.002e-3, "Fluid viscosity (Pa s); default is for water at 20C).");
content_copy
download
Use code with caution.
params.addParam<FileName>(
Generated code
"fluid_viscosity_file",
content_copy
download
Use code with caution.
Generated code
"The name of a file containing the fluid viscosity (Pa-s) as a function of temperature "
content_copy
download
Use code with caution.
Generated code
"(C); if provided the constant value is ignored.");
content_copy
download
Use code with caution.
params.addParam<Real>("fluid_thermal_conductivity",
Generated code
0.59803,
content_copy
download
Use code with caution.
Generated code
"Fluid thermal conductivity (W/(mK); default is for water at 20C).");
content_copy
download
Use code with caution.
params.addParam<FileName>(
Generated code
"fluid_thermal_conductivity_file",
content_copy
download
Use code with caution.
Generated code
"The name of a file containing fluid thermal conductivity (W/(mK)) as a function of "
content_copy
download
Use code with caution.
Generated code
"temperature (C); if provided the constant value is ignored.");
content_copy
download
Use code with caution.
params.addParam<Real>(
Generated code
"fluid_density", 998.21, "Fluid density (kg/m^3); default is for water at 20C).");
content_copy
download
Use code with caution.
params.addParam<FileName>("fluid_density_file",
Generated code
"The name of a file containing fluid density (kg/m^3) as a function "
content_copy
download
Use code with caution.
Generated code
"of temperature (C); if provided the constant value is ignored.");
content_copy
download
Use code with caution.
params.addParam<Real>(
Generated code
"fluid_specific_heat", 4157.0, "Fluid specific heat (J/(kgK); default is for water at 20C).");
content_copy
download
Use code with caution.
params.addParam<FileName>(
Generated code
"fluid_specific_heat_file",
content_copy
download
Use code with caution.
Generated code
"The name of a file containing fluid specific heat (J/(kgK) as a function of temperature "
content_copy
download
Use code with caution.
Generated code
"(C); if provided the constant value is ignored.");
content_copy
download
Use code with caution.
// Solid properties
// https://en.wikipedia.org/wiki/Stainless_steel#Properties
params.addParam<Real>("solid_thermal_conductivity",
Generated code
15.0,
content_copy
download
Use code with caution.
Generated code
"Solid thermal conductivity (W/(mK); default is for AISI/ASTIM 304 "
content_copy
download
Use code with caution.
Generated code
"stainless steel at 20C).");
content_copy
download
Use code with caution.
params.addParam<FileName>(
Generated code
"solid_thermal_conductivity_file",
content_copy
download
Use code with caution.
Generated code
"The name of a file containing solid thermal conductivity (W/(mK)) as a function of "
content_copy
download
Use code with caution.
Generated code
"temperature (C); if provided the constant value is ignored.");
content_copy
download
Use code with caution.
params.addParam<Real>(
Generated code
"solid_density",
content_copy
download
Use code with caution.
Generated code
7900,
content_copy
download
Use code with caution.
Generated code
"Solid density (kg/m^3); default is for AISI/ASTIM 304 stainless steel at 20C).");
content_copy
download
Use code with caution.
params.addParam<FileName>("solid_density_file",
Generated code
"The name of a file containing solid density (kg/m^3) as a function "
content_copy
download
Use code with caution.
Generated code
"of temperature (C); if provided the constant value is ignored.");
content_copy
download
Use code with caution.
params.addParam<Real>(
Generated code
"solid_specific_heat",
content_copy
download
Use code with caution.
Generated code
500,
content_copy
download
Use code with caution.
Generated code
"Solid specific heat (J/(kgK); default is for AISI/ASTIM 304 stainless steel at 20C).");
content_copy
download
Use code with caution.
params.addParam<FileName>(
Generated code
"solid_specific_heat_file",
content_copy
download
Use code with caution.
Generated code
"The name of a file containing solid specific heat (J/(kgK) as a function of temperature "
content_copy
download
Use code with caution.
Generated code
"(C); if provided the constant value is ignored.");
content_copy
download
Use code with caution.
return params;
}
PackedColumn::PackedColumn(const InputParameters & parameters)
: Material(parameters),

Generated code
// Get the parameters from the input file
 _input_radius(getFunction("radius")),
content_copy
download
Use code with caution.
_input_viscosity(getParam<Real>("viscosity")),
_input_porosity(getFunction("porosity")),
_temperature(adCoupledValue("temperature")),
// Fluid
_fluid_mu(getParam<Real>("fluid_viscosity")),
_fluid_k(getParam<Real>("fluid_thermal_conductivity")),
_fluid_rho(getParam<Real>("fluid_density")),
_fluid_cp(getParam<Real>("fluid_specific_heat")),
// Solid
_solid_k(getParam<Real>("solid_thermal_conductivity")),
_solid_rho(getParam<Real>("solid_density")),
_solid_cp(getParam<Real>("solid_specific_heat")),
// Material Properties being produced by this object
_permeability(declareADProperty<Real>("permeability")),
_viscosity(declareADProperty<Real>("viscosity"))
_viscosity(declareADProperty<Real>("viscosity")),
_porosity(declareADProperty<Real>("porosity")),
_thermal_conductivity(declareADProperty<Real>("thermal_conductivity")),
_specific_heat(declareADProperty<Real>("specific_heat")),
_density(declareADProperty<Real>("density"))
{
// Set data for permeability interpolation
std::vector<Real> sphere_sizes = {1, 3};
std::vector<Real> permeability = {0.8451e-9, 8.968e-9};
_permeability_interpolation.setData(sphere_sizes, permeability);
// Fluid viscosity, thermal conductivity, density, and specific heat
_use_fluid_mu_interp = initInputData("fluid_viscosity_file", _fluid_mu_interpolation);
_use_fluid_k_interp = initInputData("fluid_thermal_conductivity_file", _fluid_k_interpolation);
_use_fluid_rho_interp = initInputData("fluid_density_file", _fluid_rho_interpolation);
_use_fluid_cp_interp = initInputData("fluid_specific_heat_file", _fluid_cp_interpolation);
// Solid thermal conductivity, density, and specific heat
_use_solid_k_interp = initInputData("solid_thermal_conductivity_file", _solid_k_interpolation);
_use_solid_rho_interp = initInputData("solid_density_file", _solid_rho_interpolation);
_use_solid_cp_interp = initInputData("solid_specific_heat_file", _solid_cp_interpolation);
}
void
PackedColumn::computeQpProperties()
{

// From the paper: Table 1
std::vector<Real> sphere_sizes = {1, 3};
std::vector<Real> permeability = {0.8451e-9, 8.968e-9};
// Current temperature
ADReal temp = _temperature[_qp] - 273.15;
Real value = _input_radius.value(_t, _q_point[_qp]);
mooseAssert(value >= 1 && value <= 3,
Generated code
"The radius range must be in the range [1, 3], but " << value << " provided.");
content_copy
download
Use code with caution.
// Permeability
Real radius_value = _input_radius.value(_t, _q_point[_qp]);
mooseAssert(radius_value >= 1 && radius_value <= 3,
Generated code
"The radius range must be in the range [1, 3], but " << radius_value << " provided.");
content_copy
download
Use code with caution.
_permeability[_qp] = _permeability_interpolation.sample(radius_value);
_viscosity[_qp] = _input_viscosity;
// Porosity
Real porosity_value = _input_porosity.value(_t, _q_point[_qp]);
mooseAssert(porosity_value > 0 && porosity_value <= 1,
Generated code
"The porosity range must be in the range (0, 1], but " << porosity_value
content_copy
download
Use code with caution.
Generated code
<< " provided.");
content_copy
download
Use code with caution.
_porosity[_qp] = porosity_value;
// We'll calculate permeability using a simple linear interpolation of the two points above:
// y0 * (x1 - x) + y1 * (x - x0)
// y(x) = -------------------------------
// x1 - x0
_permeability[_qp] =
Generated code
(permeability[0] * (sphere_sizes[1] - value) + permeability[1] * (value - sphere_sizes[0])) /
content_copy
download
Use code with caution.
Generated code
(sphere_sizes[1] - sphere_sizes[0]);
content_copy
download
Use code with caution.
// Fluid properties
_viscosity[_qp] = _use_fluid_mu_interp ? _fluid_mu_interpolation.sample(temp) : _fluid_mu;
ADReal fluid_k = _use_fluid_k_interp ? _fluid_k_interpolation.sample(temp) : _fluid_k;
ADReal fluid_rho = _use_fluid_rho_interp ? _fluid_rho_interpolation.sample(temp) : _fluid_rho;
ADReal fluid_cp = _use_fluid_cp_interp ? _fluid_cp_interpolation.sample(temp) : _fluid_cp;
// Solid properties
ADReal solid_k = _use_solid_k_interp ? _solid_k_interpolation.sample(temp) : _solid_k;
ADReal solid_rho = _use_solid_rho_interp ? _solid_rho_interpolation.sample(temp) : _solid_rho;
ADReal solid_cp = _use_solid_cp_interp ? _solid_cp_interpolation.sample(temp) : _solid_cp;
// Compute the heat conduction material properties as a linear combination of
// the material properties for fluid and steel.
_thermal_conductivity[_qp] = _porosity[_qp] * fluid_k + (1.0 - _porosity[_qp]) * solid_k;
_density[_qp] = _porosity[_qp] * fluid_rho + (1.0 - _porosity[_qp]) * solid_rho;
_specific_heat[_qp] = _porosity[_qp] * fluid_cp + (1.0 - _porosity[_qp]) * solid_cp;
}
+bool
+PackedColumn::initInputData(const std::string & param_name, ADLinearInterpolation & interp)
+{

if (isParamValid(param_name))
{
const std::string & filename = getParam<FileName>(param_name);
MooseUtils::DelimitedFileReader reader(filename, &_communicator);
reader.setComment("#");
reader.read();
interp.setData(reader.getData(0), reader.getData(1));
return true;
}
return false;
+}

#pragma once

// Include the base class so it can be extended
#include "ADIntegratedBC.h"

/**

An IntegratedBC representing the "No BC" boundary condition for
heat conduction.
The residual is simply -testkgrad_u*normal... the term you get
from integration by parts. This is a standard technique for
truncating longer domains when solving the convection/diffusion
equation.
See also: Griffiths, David F. "The 'no boundary condition' outflow
boundary condition.", International Journal for Numerical Methods
in Fluids, vol. 24, no. 4, 1997, pp. 393-411.
*/
class HeatConductionOutflow : public ADIntegratedBC
{
public:
static InputParameters validParams();
HeatConductionOutflow(const InputParameters & parameters);

protected:
/**

This is called to integrate the residual across the boundary.
*/
virtual ADReal computeQpResidual() override;
/// Thermal conductivity of the material
const ADMaterialProperty<Real> & _thermal_conductivity;
};

#include "HeatConductionOutflow.h"

registerMooseObject("DarcyThermoMechApp", HeatConductionOutflow);

InputParameters
HeatConductionOutflow::validParams()
{
InputParameters params = ADIntegratedBC::validParams();
params.addClassDescription("Compute the outflow boundary condition.");
return params;
}

HeatConductionOutflow::HeatConductionOutflow(const InputParameters & parameters)
: ADIntegratedBC(parameters),
_thermal_conductivity(getADMaterialProperty<Real>("thermal_conductivity"))
{
}

ADReal
HeatConductionOutflow::computeQpResidual()
{
return -_test[_i][_qp] * _thermal_conductivity[_qp] * _grad_u[_qp] * _normals[_qp];
}

@@ -1,65 +1,112 @@
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2

nx = 100
nx = 200
ny = 10
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
[]
coord_type = RZ
rz_coord_axis = X
[]
[Variables]

[pressure]
[]
[temperature]
initial_condition = 300 # Start at room temperature
initial_condition = 300
[]
[]
[Kernels]

[darcy_pressure]
type = DarcyPressure
variable = pressure
[]
[heat_conduction]
type = ADHeatConduction
variable = temperature
type = ADHeatConduction
variable = temperature
[]
[heat_conduction_time_derivative]
type = ADHeatConductionTimeDerivative
variable = temperature
type = ADHeatConductionTimeDerivative
variable = temperature
[]
[heat_convection]
type = DarcyAdvection
variable = temperature
pressure = pressure
[]
[]
[BCs]
[inlet_temperature]

type = FunctionDirichletBC
variable = temperature
boundary = left
function = 'if(t<0,350+50*t,350)'
type = FunctionDirichletBC
variable = temperature
boundary = left
function = 'if(t<0,350+50*t,350)'
[]
[outlet_temperature]
type = HeatConductionOutflow
variable = temperature
boundary = right
[]
[inlet]
type = DirichletBC
variable = pressure
boundary = left
value = 4000
[]
[outlet]
type = DirichletBC
variable = pressure
boundary = right
value = 0
[]
[]
-[Materials/steel]
-- type = ADGenericConstantMaterial
-- prop_names = 'thermal_conductivity specific_heat density'
-- prop_values = '18 0.466 8000' # W/m*K, J/kg-K, kg/m^3 @ 296K
+[Materials/column]

type = PackedColumn
temperature = temperature
radius = 1
+[]
+[AuxVariables/velocity]

order = CONSTANT
family = MONOMIAL_VEC
+[]
+[AuxKernels/velocity]

type = DarcyVelocity
variable = velocity
execute_on = timestep_end
pressure = pressure
[]
[Problem]
type = FEProblem
[]

[Executioner]
type = Transient

num_steps = 10
end_time = 100
dt = 0.25
start_time = -1
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
automatic_scaling = true
steady_state_tolerance = 1e-5
steady_state_detection = true
[TimeStepper]
type = FunctionDT
function = 'if(t<0,0.1,0.25)'
[]
[]
[Outputs]

exodus = true
exodus = true
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step06_coupled_darcy_heat_conduction
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step6a_coupled.i

@@ -1,112 +1,117 @@
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2
nx = 200
ny = 10
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
[]
coord_type = RZ
rz_coord_axis = X
[]

[Variables]
[pressure]
[]
[temperature]
initial_condition = 300 # Start at room temperature
[]
[]

[Kernels]
[darcy_pressure]
type = DarcyPressure
variable = pressure
[]
[heat_conduction]
type = ADHeatConduction
variable = temperature
[]
[heat_conduction_time_derivative]
type = ADHeatConductionTimeDerivative
variable = temperature
[]
[heat_convection]
type = DarcyAdvection
variable = temperature
pressure = pressure
[]
[]

[BCs]
[inlet_temperature]
type = FunctionDirichletBC
variable = temperature
boundary = left
function = 'if(t<0,350+50*t,350)'
[]
[outlet_temperature]
type = HeatConductionOutflow
variable = temperature
boundary = right
[]
[inlet]

type = DirichletBC
type = FunctionDirichletBC
variable = pressure
boundary = left
value = 4000 # (Pa) From Figure 2 from paper. First data point for 1mm spheres.
function = 2000sin(0.466pi*t) # Inlet signal from Fig. 3
[]
[outlet]
type = DirichletBC
type = FunctionDirichletBC
variable = pressure
boundary = right
value = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
function = 2000cos(0.466pi*t) # Outlet signal from Fig. 3
[]
[]
[Materials/column]
type = PackedColumn

radius = 1
temperature = temperature
radius = 1
fluid_viscosity_file = data/water_viscosity.csv
fluid_density_file = data/water_density.csv
fluid_thermal_conductivity_file = data/water_thermal_conductivity.csv
fluid_specific_heat_file = data/water_specific_heat.csv
outputs = exodus
[]
[AuxVariables/velocity]
order = CONSTANT
family = MONOMIAL_VEC
[]

[AuxKernels/velocity]
type = DarcyVelocity
variable = velocity
execute_on = timestep_end
pressure = pressure
[]

[Problem]
type = FEProblem
[]

[Executioner]
type = Transient
end_time = 100
dt = 0.25
start_time = -1

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

automatic_scaling = true
steady_state_tolerance = 1e-5
steady_state_detection = true

[TimeStepper]
type = FunctionDT

function = 'if(t<0,0.1,0.25)'
function = 'if(t<0,0.1,(2pi/(0.466pi))/16)' # dt to always hit the peaks of sine/cosine BC
[]
[]
[Outputs]

exodus = true
exodus = true
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step06_coupled_darcy_heat_conduction
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step6b_transient_inflow.i

@@ -1,112 +1,82 @@
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2
nx = 200
ny = 10
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
[]
coord_type = RZ
rz_coord_axis = X
[]

[Variables]
-- [pressure]
-- []
[temperature]

initial_condition = 300 # Start at room temperature
initial_condition = 300
[]
[]
+[AuxVariables/pressure]
+[]
+
+[AuxKernels/pressure]

type = FunctionAux
variable = pressure
function = '4000 - 3000 * x - 3000 * txx*y'
execute_on = timestep_end
+[]
[Kernels]
-- [darcy_pressure]
-- type = DarcyPressure
-- variable = pressure
-- []
[heat_conduction]

type = ADHeatConduction
variable = temperature
type = ADHeatConduction
variable = temperature
[]
[heat_conduction_time_derivative]
type = ADHeatConductionTimeDerivative
variable = temperature
type = ADHeatConductionTimeDerivative
variable = temperature
[]
[heat_convection]
type = DarcyAdvection
variable = temperature
pressure = pressure
[]
[]
[BCs]
[inlet_temperature]
-- type = FunctionDirichletBC
-+ type = DirichletBC

variable = temperature
boundary = left
-- function = 'if(t<0,350+50*t,350)'
-+ value = 350
type = DirichletBC
variable = temperature
boundary = left
value = 350
[]
[outlet_temperature]
type = HeatConductionOutflow
variable = temperature
boundary = right
[]
-- [inlet]
-- type = DirichletBC
-- variable = pressure
-- boundary = left
-- value = 4000 # (Pa) From Figure 2 from paper. First data point for 1mm spheres.
-- []
-- [outlet]
-- type = DirichletBC
-- variable = pressure
-- boundary = right
-- value = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
-- []
[]
[Materials/column]
type = PackedColumn
-- temperature = temperature
-+ temperature = 293.15 # 20C

temperature = 293.15
radius = 1
--[]
--
-[AuxVariables/velocity]
-- order = CONSTANT
-- family = MONOMIAL_VEC
--[]
--
-[AuxKernels/velocity]
-- type = DarcyVelocity
-- variable = velocity
-- execute_on = timestep_end
-- pressure = pressure
[]
[Problem]
type = FEProblem
[]

[Executioner]
type = Transient
-- end_time = 100
-- dt = 0.25
-- start_time = -1

num_steps = 300
dt = 0.1
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
--
-- automatic_scaling = true
-- steady_state_tolerance = 1e-5
-- steady_state_detection = true
--
-- [TimeStepper]
-- type = FunctionDT
-- function = 'if(t<0,0.1,0.25)'
-- []
[]
[Outputs]

exodus = true
exodus = true
[]

[Outputs]
execute_on = TIMESTEP_END

[Outputs]
[exo] # 文件名后缀
type = Exodus
execute_on = 'LINEAR TIMESTEP_END'
output_material_properties = true
[]
[]

[Debug]
show_material_props = true
[]

[Debug]
show_actions = true
[]

[Debug]
show_execution_order = ALWAYS
[]

[Debug]
show_var_residual_norms = true
[]

[Outputs]
exodus = true
execute_on = INITIAL
[]

[Adaptivity]
initial_steps = 2
cycles_per_step = 2
marker = marker
initial_marker = marker
max_h_level = 2
[Indicators]
[indicator]
type = GradientJumpIndicator
variable = u
[]
[]
[Markers]
[marker]
type = ErrorFractionMarker
indicator = indicator
coarsen = 0.1
refine = 0.7
[]
[]
[]

@@ -1,112 +1,112 @@
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2

nx = 200
ny = 10
nx = 30
ny = 3
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
[]
coord_type = RZ
rz_coord_axis = X
[]
[Variables]
[pressure]
[]
[temperature]

initial_condition = 300 # Start at room temperature
initial_condition = 300
[]
[]
[Kernels]
[darcy_pressure]
type = DarcyPressure
variable = pressure
[]
[heat_conduction]

type = ADHeatConduction
variable = temperature
type = ADHeatConduction
variable = temperature
[]
[heat_conduction_time_derivative]
type = ADHeatConductionTimeDerivative
variable = temperature
type = ADHeatConductionTimeDerivative
variable = temperature
[]
[heat_convection]
type = DarcyAdvection
variable = temperature
pressure = pressure
[]
[]
[BCs]
[inlet_temperature]

type = FunctionDirichletBC
variable = temperature
boundary = left
function = 'if(t<0,350+50*t,350)'
type = FunctionDirichletBC
variable = temperature
boundary = left
function = 'if(t<0,350+50*t,350)'
[]
[outlet_temperature]
type = HeatConductionOutflow
variable = temperature
boundary = right
[]
[inlet]
type = DirichletBC
variable = pressure
boundary = left
value = 4000 # (Pa) From Figure 2 from paper. First data point for 1mm spheres.
type = DirichletBC
variable = pressure
boundary = left
value = 4000
[]
[outlet]
type = DirichletBC
variable = pressure
boundary = right
value = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
type = DirichletBC
variable = pressure
boundary = right
value = 0
[]
[]
[Materials/column]
type = PackedColumn
temperature = temperature
radius = 1
[]

[AuxVariables/velocity]
order = CONSTANT
family = MONOMIAL_VEC
[]

[AuxKernels/velocity]
type = DarcyVelocity
variable = velocity
execute_on = timestep_end
pressure = pressure
[]

[Problem]
type = FEProblem
[]

[Executioner]
type = Transient
end_time = 100
dt = 0.25
start_time = -1

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

automatic_scaling = true
steady_state_tolerance = 1e-5
steady_state_detection = true

[TimeStepper]
type = FunctionDT
function = 'if(t<0,0.1,0.25)'
[]
[]

[Outputs]

exodus = true
exodus = true
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step07_adaptivity
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step7a_coarse.i

@@ -1,112 +1,113 @@
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2
nx = 30
ny = 3
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
[]
coord_type = RZ
rz_coord_axis = X

uniform_refine = 3
[]
[Variables]
[pressure]
[]
[temperature]

initial_condition = 300 # Start at room temperature
initial_condition = 300
[]
[]
[Kernels]
[darcy_pressure]
type = DarcyPressure
variable = pressure
[]
[heat_conduction]

type = ADHeatConduction
variable = temperature
type = ADHeatConduction
variable = temperature
[]
[heat_conduction_time_derivative]
type = ADHeatConductionTimeDerivative
variable = temperature
type = ADHeatConductionTimeDerivative
variable = temperature
[]
[heat_convection]
type = DarcyAdvection
variable = temperature
pressure = pressure
[]
[]
[BCs]
[inlet_temperature]

type = FunctionDirichletBC
variable = temperature
boundary = left
function = 'if(t<0,350+50*t,350)'
type = FunctionDirichletBC
variable = temperature
boundary = left
function = 'if(t<0,350+50*t,350)'
[]
[outlet_temperature]
type = HeatConductionOutflow
variable = temperature
boundary = right
[]
[inlet]
type = DirichletBC
variable = pressure
boundary = left
value = 4000 # (Pa) From Figure 2 from paper. First data point for 1mm spheres.
type = DirichletBC
variable = pressure
boundary = left
value = 4000
[]
[outlet]
type = DirichletBC
variable = pressure
boundary = right
value = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
type = DirichletBC
variable = pressure
boundary = right
value = 0
[]
[]
[Materials/column]
type = PackedColumn
temperature = temperature
radius = 1
[]

[AuxVariables/velocity]
order = CONSTANT
family = MONOMIAL_VEC
[]

[AuxKernels/velocity]
type = DarcyVelocity
variable = velocity
execute_on = timestep_end
pressure = pressure
[]

[Problem]
type = FEProblem
[]

[Executioner]
type = Transient
end_time = 100
dt = 0.25
start_time = -1

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

automatic_scaling = true
steady_state_tolerance = 1e-5
steady_state_detection = true

[TimeStepper]
type = FunctionDT
function = 'if(t<0,0.1,0.25)'
[]
[]

[Outputs]

exodus = true
exodus = true
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step07_adaptivity
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step7b_fine.i

@@ -1,113 +1,129 @@
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2
nx = 30
ny = 3
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
[]
coord_type = RZ
rz_coord_axis = X

uniform_refine = 3
uniform_refine = 3 # This is now the max_h_level, not a uniform refinement
[]
[Variables]
[pressure]
[]
[temperature]

initial_condition = 300 # Start at room temperature
initial_condition = 300
[]
[]
[Kernels]
[darcy_pressure]
type = DarcyPressure
variable = pressure
[]
[heat_conduction]

type = ADHeatConduction
variable = temperature
type = ADHeatConduction
variable = temperature
[]
[heat_conduction_time_derivative]
type = ADHeatConductionTimeDerivative
variable = temperature
type = ADHeatConductionTimeDerivative
variable = temperature
[]
[heat_convection]
type = DarcyAdvection
variable = temperature
pressure = pressure
[]
[]
[BCs]
[inlet_temperature]

type = FunctionDirichletBC
variable = temperature
boundary = left
function = 'if(t<0,350+50*t,350)'
type = FunctionDirichletBC
variable = temperature
boundary = left
function = 'if(t<0,350+50*t,350)'
[]
[outlet_temperature]
type = HeatConductionOutflow
variable = temperature
boundary = right
[]
[inlet]
type = DirichletBC
variable = pressure
boundary = left
value = 4000 # (Pa) From Figure 2 from paper. First data point for 1mm spheres.
type = DirichletBC
variable = pressure
boundary = left
value = 4000
[]
[outlet]
type = DirichletBC
variable = pressure
boundary = right
value = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
type = DirichletBC
variable = pressure
boundary = right
value = 0
[]
[]
[Materials/column]
type = PackedColumn
temperature = temperature
radius = 1
[]

[AuxVariables/velocity]
order = CONSTANT
family = MONOMIAL_VEC
[]

[AuxKernels/velocity]
type = DarcyVelocity
variable = velocity
execute_on = timestep_end
pressure = pressure
[]

[Problem]
type = FEProblem
[]

[Executioner]
type = Transient
end_time = 100
dt = 0.25
start_time = -1

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

automatic_scaling = true
steady_state_tolerance = 1e-5
steady_state_detection = true

[TimeStepper]
type = FunctionDT
function = 'if(t<0,0.1,0.25)'
[]
[]

[Outputs]

exodus = true
exodus = true
+[]
+[Adaptivity]

marker = error_frac
max_h_level = 3
[Indicators/temperature_jump]
type = GradientJumpIndicator
variable = temperature
scale_by_flux_faces = true
[]
[Markers/error_frac]
type = ErrorFractionMarker
coarsen = 0.15
indicator = temperature_jump
refine = 0.7
[]
[]

cd ~/projects/moose/tutorials/darcy_thermo_mech/step07_adaptivity
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step7a_adapt.i

@@ -1,129 +1,161 @@
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2

nx = 30
ny = 3
nx = 40
ny = 4
xmax = 0.304 # Length of test chamber
ymax = 0.0257 # Test chamber radius
[]
[bottom]
type = SubdomainBoundingBoxGenerator
input = gmg
location = inside
bottom_left = '0 0 0'
top_right = '0.304 0.01285 0'
block_id = 1
[]
coord_type = RZ
rz_coord_axis = X
uniform_refine = 3
uniform_refine = 3 # This is now the max_h_level, not a uniform refinement
[]
[Variables]
[pressure]
[]
[temperature]

initial_condition = 300 # Start at room temperature
initial_condition = 300
[]
[]
[Kernels]
[darcy_pressure]
type = DarcyPressure
variable = pressure
[]
[heat_conduction]

type = ADHeatConduction
variable = temperature
type = ADHeatConduction
variable = temperature
[]
[heat_conduction_time_
Generated org
derivative]
-    type = ADHeatConductionTimeDerivative
-    variable = temperature
+    type = ADHeatConductionTimeDerivative
+    variable = temperature
   []
   [heat_convection]
     type = DarcyAdvection
     variable = temperature
     pressure = pressure
   []
 []

 [BCs]
   [inlet_temperature]
-    type = FunctionDirichletBC
-    variable = temperature
-    boundary = left
-    function = 'if(t<0,350+50*t,350)'
+    type = FunctionDirichletBC
+    variable = temperature
+    boundary = left
+    function = 'if(t<0,350+50*t,350)'
   []
   [outlet_temperature]
     type = HeatConductionOutflow
     variable = temperature
     boundary = right
   []
   [inlet]
-    type = DirichletBC
-    variable = pressure
-    boundary = left
-    value = 4000 # (Pa) From Figure 2 from paper.  First data point for 1mm spheres.
+    type = DirichletBC
+    variable = pressure
+    boundary = left
+    value = 4000
   []
   [outlet]
-    type = DirichletBC
-    variable = pressure
-    boundary = right
-    value = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
+    type = DirichletBC
+    variable = pressure
+    boundary = right
+    value = 0
   []
 []

-[Materials/column]
--  type = PackedColumn
--  temperature = temperature
--  radius = 1
+[Materials]
+  viscosity_file = data/water_viscosity.csv
+  density_file = data/water_density.csv
+  thermal_conductivity_file = data/water_thermal_conductivity.csv
+  specific_heat_file = data/water_specific_heat.csv
+
+  [column_bottom]
+    type = PackedColumn
+    block = 1
+    radius = 1.15
+    temperature = temperature
+    fluid_viscosity_file = ${viscosity_file}
+    fluid_density_file = ${density_file}
+    fluid_thermal_conductivity_file = ${thermal_conductivity_file}
+    fluid_specific_heat_file = ${specific_heat_file}
+  []
+  [column_top]
+    type = PackedColumn
+    block = 0
+    radius = 1
+    temperature = temperature
+    porosity = '0.25952 + 0.7*x/0.304'
+    fluid_viscosity_file = ${viscosity_file}
+    fluid_density_file = ${density_file}
+    fluid_thermal_conductivity_file = ${thermal_conductivity_file}
+    fluid_specific_heat_file = ${specific_heat_file}
+  []
 []

 [AuxVariables/velocity]
   order = CONSTANT
   family = MONOMIAL_VEC
 []

 [AuxKernels/velocity]
   type = DarcyVelocity
   variable = velocity
   execute_on = timestep_end
   pressure = pressure
 []

 [Problem]
   type = FEProblem
 []

 [Executioner]
   type = Transient
   end_time = 100
   dt = 0.25
   start_time = -1

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

   automatic_scaling = true
   steady_state_tolerance = 1e-5
   steady_state_detection = true

   [TimeStepper]
     type = FunctionDT
     function = 'if(t<0,0.1,0.25)'
   []
 []

-[Outputs]
--  exodus = true
+[Outputs/out]
+  type = Exodus
+  output_material_properties = true
 []

 [Adaptivity]
   marker = error_frac
   max_h_level = 3
   [Indicators/temperature_jump]
-    type = GradientJumpIndicator
-    variable = temperature
-    scale_by_flux_faces = true
+    type = GradientJumpIndicator
+    variable = temperature
+    scale_by_flux_faces = true
   []
   [Markers/error_frac]
-    type = ErrorFractionMarker
--    coarsen = 0.15
--    indicator = temperature_jump
--    refine = 0.7
+    type = ErrorFractionMarker
+    coarsen = 0.025
+    indicator = temperature_jump
+    refine = 0.9
   []
 []

Real averageValue(SubdomainID block) const;

InputParameters
BlockAverageDiffusionMaterial::validParams()
{
  InputParameters params = Material::validParams();
  params.addRequiredParam<UserObjectName>("block_average_userobject", "Computes the ...");
  return params;
}

const BlockAverageValue & _block_average_value;

BlockAverageDiffusionMaterial::BlockAverageDiffusionMaterial(const InputParameters & parameters) :
    Material(parameters),
    _block_average_value(getUserObject<BlockAverageValue>("block_average_userobject"))
{}

_diffusivity[_qp] = 0.5 * _block_average_value.averageValue(_current_elem->subdomain_id());

[Output]
  csv = true
[]

const PostprocessorValue & _pps_value;

_pps_value(getPostprocessorValue("postprocessor")),

params.addParam<PostprocessorName>("postprocessor", 1.2345, "Doc String");

VectorPostprocessorValue & _pid;

_pid(declareVector("pid")),

[Output]
  csv = true
[]

const VectorPostprocessorValue & _x_values;

_x_values(getVectorPostprocessorValue("vectorpostprocessor", _x_name)),

 @@ -1,112 +1,137 @@
 [Mesh]
   [gmg]
     type = GeneratedMeshGenerator
     dim = 2
-    nx = 200
-    ny = 10
+    nx = 30
+    ny = 3
     xmax = 0.304 # Length of test chamber
     ymax = 0.0257 # Test chamber radius
   []
   coord_type = RZ
   rz_coord_axis = X
+  uniform_refine = 2
 []

 [Variables]
   [pressure]
   []
   [temperature]
-    initial_condition = 300 # Start at room temperature
+    initial_condition = 300
   []
 []

 [Kernels]
   [darcy_pressure]
     type = DarcyPressure
     variable = pressure
   []
   [heat_conduction]
-    type = ADHeatConduction
-    variable = temperature
+    type = ADHeatConduction
+    variable = temperature
   []
   [heat_conduction_time_derivative]
-    type = ADHeatConductionTimeDerivative
-    variable = temperature
+    type = ADHeatConductionTimeDerivative
+    variable = temperature
   []
   [heat_convection]
     type = DarcyAdvection
     variable = temperature
     pressure = pressure
   []
 []

 [BCs]
   [inlet_temperature]
-    type = FunctionDirichletBC
-    variable = temperature
-    boundary = left
-    function = 'if(t<0,350+50*t,350)'
+    type = FunctionDirichletBC
+    variable = temperature
+    boundary = left
+    function = 'if(t<0,350+50*t,350)'
   []
   [outlet_temperature]
     type = HeatConductionOutflow
     variable = temperature
     boundary = right
   []
   [inlet]
-    type = DirichletBC
-    variable = pressure
-    boundary = left
-    value = 4000 # (Pa) From Figure 2 from paper.  First data point for 1mm spheres.
+    type = DirichletBC
+    variable = pressure
+    boundary = left
+    value = 4000
   []
   [outlet]
-    type = DirichletBC
-    variable = pressure
-    boundary = right
-    value = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
+    type = DirichletBC
+    variable = pressure
+    boundary = right
+    value = 0
   []
 []

 [Materials/column]
   type = PackedColumn
   temperature = temperature
   radius = 1
+  porosity = '0.25952 + 0.7*y/0.0257'
+[]
+
+[Postprocessors]
+  [average_temperature]
+    type = ElementAverageValue
+    variable = temperature
+  []
+  [outlet_heat_flux]
+    type = ADSideDiffusiveFluxIntegral
+    variable = temperature
+    boundary = right
+    diffusivity = thermal_conductivity
+  []
+[]
+
+[VectorPostprocessors/temperature_sample]
+  type = LineValueSampler
+  num_points = 500
+  start_point = '0.1 0      0'
+  end_point =   '0.1 0.0257 0'
+  variable = temperature
+  sort_by = y
 []

 [AuxVariables/velocity]
   order = CONSTANT
   family = MONOMIAL_VEC
 []

 [AuxKernels/velocity]
   type = DarcyVelocity
   variable = velocity
   execute_on = timestep_end
   pressure = pressure
 []

 [Problem]
   type = FEProblem
 []

 [Executioner]
   type = Transient
   end_time = 100
   dt = 0.25
   start_time = -1

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

   automatic_scaling = true
   steady_state_tolerance = 1e-5
   steady_state_detection = true

   [TimeStepper]
     type = FunctionDT
     function = 'if(t<0,0.1,0.25)'
   []
 []

 [Outputs]
-  exodus = true
+  exodus = true
+  csv = true
 []

cd ~/projects/moose/tutorials/darcy_thermo_mech/step08_postprocessors
make -j 12 # 使用您系统的处理器数量
cd problems
../darcy_thermo_mech-opt -i step8.i