发布于
RB 学习笔记
从二维热方程出发，用 Fortran 逐步实现 FGT 热传播、自适应四叉树、无滑移 Stokes、Navier–Stokes 与 Rayleigh–Bénard 对流求解器。
Chapter 1：混合边界热方程的直接传播
1. 目标
考虑单位方盒上的温度扰动 \(\theta=T-(1-z)\)：
\[
\theta_t=D\Delta\theta,\qquad
\theta(x,0,t)=\theta(x,1,t)=0,\qquad
\theta(x+1,z,t)=\theta(x,z,t).
\]
2. 完整代码
module heat_direct
use iso_fortran_env, only: real64
implicit none
private
integer, parameter, public :: rk = real64
real(rk), parameter :: pi = acos(-1.0_rk)
public :: heat_step
contains
function gaussian(r, delta) result(g)
real(rk), intent(in) :: r, delta
real(rk) :: g
g = exp(-r*r/delta) / sqrt(pi * delta)
end function gaussian
subroutine heat_step(source, x, z, D, h, images, u)
real(rk), intent(in) :: source(:, :)
real(rk), intent(in) :: x(:), z(:), D, h
integer, intent(in) :: images
real(rk), intent(out) :: u(size(x), size(z))
real(rk) :: kx(size(x), size(source, 1))
real(rk) :: kz(size(z), size(source, 2))
real(rk) :: delta, q, shift
integer :: n, i, j, p
n = size(source, 1)
if (n < 1 .or. size(source, 2) /= n) error stop "source must be a square array"
if (D <= 0.0_rk .or. h <= 0.0_rk .or. images < 0) error stop "invalid heat parameters"
delta = 4.0_rk * D * h
kx = 0.0_rk
kz = 0.0_rk
do i = 1, size(x)
do j = 1, n
q = (real(j, rk) - 0.5_rk) / real(n, rk)
do p = -images, images
shift = real(p, rk)
kx(i, j) = kx(i, j) + gaussian(x(i) - q + shift, delta)
end do
end do
end do
do i = 1, size(z)
do j = 1, n
q = (real(j, rk) - 0.5_rk) / real(n, rk)
do p = -images, images
shift = 2.0_rk * real(p, rk)
kz(i, j) = kz(i, j) + gaussian(z(i) - q + shift, delta)-gaussian(z(i) + q + shift, delta)
end do
end do
end do
u = matmul(kx, matmul(source, transpose(kz))) / real(n, rk)**2
end subroutine heat_step
end module heat_direct
3. 传播公式
二维热方程经过时间 \(h\) 的自由空间热核为
\[
G_h(x-\xi,z-\eta)
=
\frac{1}{4\pi Dh}
\exp\left[-\frac{(x-\xi)^2+(z-\eta)^2}{4Dh}\right].
\]
4. 从积分公式到矩阵乘法
4.1 中点采样与接口约定
混合边界下的精确传播写成
\[
\theta(x,z,t+h)
=
\int_0^1\int_0^1
K_x(x,\xi)\,
\theta(\xi,\eta,t)\,
K_z(z,\eta)
\,d\eta\,d\xi.
\]
5. 模块中的执行流程
source、x、z、D、h、images
→ 检查源数组形状和参数范围
→ 计算 delta = 4Dh
→ 累加水平周期核 kx
→ 累加竖直奇反射核 kz
→ 两次矩阵乘法并乘面积权重
→ 返回目标位置上的温度扰动 u
当 \(N_x,N_z\) 与 \(n\) 同量级时，两次稠密矩阵乘法的算术量为 \(O(n^3)\)；直接枚举全部二维源与目标则为 \(O(n^4)\)。核矩阵构造还需要 \(O((2M+1)n(N_x+N_z))\) 次一维核计算。
Chapter 2：Chebyshev 节点与盒内插值
1. 目标
将一个正方形盒子中的场表示为张量积多项式。每个方向使用 \(k\) 个 Chebyshev 根节点，保存 \(k\times k\) 个函数值，再通过插值计算任意盒内位置的值。
cheb_box 提供三个操作：生成节点 cheb_nodes、构造一维插值矩阵 cheb_matrix、计算二维插值 evaluate_box。后续的源场表示、网格细化和粗化都以这一表示为基础。
2. 完整代码
module cheb_box
use iso_fortran_env, only: real64
implicit none
private
integer, parameter, public :: rk = real64
real(rk), parameter :: pi = acos(-1.0_rk)
public :: cheb_nodes, cheb_matrix, evaluate_box
contains
subroutine cheb_nodes(k, nodes)
integer, intent(in) :: k
real(rk), intent(out) :: nodes(k)
integer :: j
if (k < 2) error stop "k must be at least 2"
do j = 1, k
nodes(j) = cos(pi * (real(j, rk) - 0.5_rk) / real(k, rk))
end do
end subroutine cheb_nodes
subroutine cheb_matrix(k, target, matrix)
integer, intent(in) :: k
real(rk), intent(in) :: target(:)
real(rk), intent(out) :: matrix(size(target), k)
real(rk) :: nodes(k), weights(k), angle
integer :: i, j
call cheb_nodes(k, nodes)
do j = 1, k
angle = pi * (real(j, rk) - 0.5_rk) / real(k, rk)
weights(j) = (-1.0_rk)**(j - 1) * sin(angle)
end do
do i = 1, size(target)
j = minloc(abs(target(i) - nodes), dim=1)
if (abs(target(i) - nodes(j)) <= 8.0_rk*epsilon(1.0_rk)) then
matrix(i, :) = 0.0_rk
matrix(i, j) = 1.0_rk
else
matrix(i, :) = weights / (target(i) - nodes)
matrix(i, :) = matrix(i, :) / sum(matrix(i, :))
end if
end do
end subroutine cheb_matrix
subroutine evaluate_box(values, center, width, x, z, u)
real(rk), intent(in) :: values(:, :)
real(rk), intent(in) :: center(2), width, x(:), z(:)
real(rk), intent(out) :: u(size(x), size(z))
real(rk) :: bx(size(x), size(values, 1))
real(rk) :: bz(size(z), size(values, 2))
integer :: k
k = size(values, 1)
if (size(values, 2) /= k .or. width <= 0.0_rk) error stop "invalid box"
call cheb_matrix(k, 2.0_rk*(x - center(1))/width, bx)
call cheb_matrix(k, 2.0_rk*(z - center(2))/width, bz)
u = matmul(bx, matmul(values, transpose(bz)))
end subroutine evaluate_box
end module cheb_box
Chapter 3：单盒 Gaussian 积分与 FGT 接口
1. 目标
把一个盒子上的连续密度交给FGT库，计算
\[
P(\boldsymbol x)
=
\int_B
e^{-|\boldsymbol x-\boldsymbol y|^2/\delta}
\,\sigma(\boldsymbol y)\,d\boldsymbol y.
\]
gaussian_box 负责准备数据并调用 boxfgt。输入是盒内多项式密度的节点值，输出是同一盒节点上的 Gaussian 体积势。这个核尚未包含热方程的 \(1/(\pi\delta)\) 归一化因子。
2. 完整代码
module fgt_single_box
use cheb_box, only: rk
use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
implicit none
private
public :: gaussian_box, boxfgt
interface
subroutine boxfgt(nd, d, delta, eps, ipoly, iperiod, k, np, nb, nlev, ltree, itree, iptr, centers, boxsize, fvals, ifpgh, pot, grad, hess, ifnewtree, nt, targs, ifpghtarg, pote, grade, hesse, timeinfo)
import rk
integer :: nd, d, ipoly, k, np, nb, nlev, ltree
integer, intent(inout) :: iperiod
integer :: itree(ltree), iptr(8), ifpgh, ifnewtree, nt, ifpghtarg
real(rk) :: delta, eps, centers(d,nb), boxsize(0:nlev)
real(rk) :: fvals(nd,np,nb), pot(nd,np,nb)
real(rk) :: grad(nd,d,np,*), hess(nd,d*(d+1)/2,np,*)
real(rk) :: targs(d,nt), pote(nd,nt)
real(rk) :: grade(nd,d,*), hesse(nd,d*(d+1)/2,*), timeinfo(*)
end subroutine boxfgt
end interface
contains
subroutine gaussian_box(k, values, center, width, delta, tolerance, potential)
integer, intent(in) :: k
real(rk), intent(in) :: values(k,k), center(2), width, delta, tolerance
real(rk), intent(out) :: potential(k,k)
integer :: tree(20), ptr(8), nb, nlev, ltree, np, iperiod
real(rk) :: centers(2,1), boxsize(0:0), dd, eps
real(rk) :: f(1,k*k,1), p(1,k*k,1)
real(rk) :: grad(1,2,k*k,1), hess(1,3,k*k,1)
real(rk) :: target(2,1), pt(1,1), gt(1,2,1), ht(1,3,1), times(100)
if (k < 2 .or. k > 32) error stop "unsupported order"
if (.not. all(ieee_is_finite([center, width, delta, tolerance]))) error stop "nonfinite parameter"
if (width <= 0.0_rk .or. delta <= 0.0_rk) error stop "invalid geometry or delta"
if (tolerance < 1.0e-14_rk .or. tolerance >= 0.1_rk) error stop "use a tolerance between 1e-14 and 0.1"
if (.not. all(ieee_is_finite(values))) error stop "nonfinite density"
nb = 1
nlev = 0
ltree = 20
np = k*k
iperiod = 0
ptr = [1, 3, 4, 5, 6, 10, 11, 20]
tree = -1
tree(1:2) = [1, 1]
tree(ptr(2)) = 0
tree(ptr(4)) = 0
tree(ptr(6)) = 1
tree(ptr(7)) = 1
centers(:,1) = center
boxsize(0) = width
dd = delta
eps = tolerance
f(1,:,1) = reshape(values, [np])
p = 0.0_rk
grad = 0.0_rk
hess = 0.0_rk
target = 0.0_rk
pt = 0.0_rk
gt = 0.0_rk
ht = 0.0_rk
times = 0.0_rk
call boxfgt(1, 2, dd, eps, 1, iperiod, k, np, nb, nlev, ltree, tree, ptr, centers, boxsize, f, 1, p, grad, hess, 0, 0, target, 0, pt, gt, ht, times)
if (.not. all(ieee_is_finite(p))) error stop "nonfinite FGT output"
potential = reshape(p(1,:,1), [k,k])
end subroutine gaussian_box
end module fgt_single_box
Chapter 4：固定树上的混合边界热传播
1. 目标
将单位物理盒延拓成 \([0,2]\times[0,2]\) 上的周期场，用一次周期 FGT 实现水平周期、上下零温扰边界的热传播。
这一步把 Chapter 1 中显式构造的混合边界核，改写成“密度延拓 → 周期 Gaussian 积分 → 限制回物理域”。
2. 完整代码
module heat_fgt_mixed
use cheb_box, only: rk
use fgt_single_box, only: boxfgt
use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
implicit none
private
public :: heat_step_fgt
contains
subroutine heat_step_fgt(k, source, diffusion, step, tolerance, output)
integer, intent(in) :: k
real(rk), intent(in) :: source(k,k), diffusion, step, tolerance
real(rk), intent(out) :: output(k,k)
integer :: tree(90), ptr(8), nb, nlev, ltree, np, iperiod
real(rk), parameter :: pi = acos(-1.0_rk)
real(rk) :: centers(2,5), boxsize(0:1), delta, eps
real(rk) :: f(1,k*k,5), p(1,k*k,5)
real(rk) :: grad(1,2,k*k,5), hess(1,3,k*k,5)
real(rk) :: target(2,1), pt(1,1), gt(1,2,1), ht(1,3,1), times(100)
if (k < 2 .or. k > 32) error stop "unsupported order"
if (.not. all(ieee_is_finite([diffusion, step, tolerance]))) error stop "nonfinite heat parameter"
if (diffusion <= 0.0_rk .or. step <= 0.0_rk) error stop "diffusion and step must be positive"
if (tolerance < 1.0e-14_rk .or. tolerance >= 0.1_rk) error stop "invalid tolerance"
if (.not. all(ieee_is_finite(source))) error stop "nonfinite source"
nb = 5
nlev = 1
ltree = 90
np = k*k
iperiod = 1
delta = 4.0_rk*diffusion*step
eps = tolerance
if (.not. ieee_is_finite(delta) .or. delta <= 0.0_rk) error stop "invalid delta"
ptr = [1, 5, 10, 15, 20, 40, 45, 90]
tree = -1
tree(1:4) = [1, 1, 2, 5]
tree(ptr(2):ptr(2)+4) = [0, 1, 1, 1, 1]
tree(ptr(3):ptr(3)+4) = [-1, 1, 1, 1, 1]
tree(ptr(4):ptr(4)+4) = [4, 0, 0, 0, 0]
tree(ptr(5):ptr(5)+3) = [2, 3, 4, 5]
tree(ptr(6):ptr(6)+4) = 0
centers(:,1) = [1.0_rk, 1.0_rk]
centers(:,2) = [0.5_rk, 0.5_rk]
centers(:,3) = [1.5_rk, 0.5_rk]
centers(:,4) = [0.5_rk, 1.5_rk]
centers(:,5) = [1.5_rk, 1.5_rk]
boxsize = [2.0_rk, 1.0_rk]
f = 0.0_rk
f(1,:,2) = reshape(source, [np])
f(1,:,3) = f(1,:,2)
f(1,:,4) = reshape(-source(:,k:1:-1), [np])
f(1,:,5) = f(1,:,4)
p = 0.0_rk
grad = 0.0_rk
hess = 0.0_rk
target = 0.0_rk
pt = 0.0_rk
gt = 0.0_rk
ht = 0.0_rk
times = 0.0_rk
call boxfgt(1, 2, delta, eps, 1, iperiod, k, np, nb, nlev, ltree, tree, ptr, centers, boxsize, f, 1, p, grad, hess, 0, 0, target, 0, pt, gt, ht, times)
output = reshape(p(1,:,2), [k,k])/(pi*delta)
if (.not. all(ieee_is_finite(output))) error stop "nonfinite output"
end subroutine heat_step_fgt
end module heat_fgt_mixed
Chapter 5：带源热方程的二阶时间推进
1. 目标
在扩散方程中加入已知源项：
\[
\theta_t=D\Delta\theta+f(x,z,t).
\]
forced_step 接收当前状态及时间区间两端的源值，调用已有热传播算子完成一步更新。state 是温度扰动，force0、force1 才是这里的 PDE 源项。
2. 完整代码
module heat_forced
use cheb_box, only: rk
use heat_fgt_mixed, only: heat_step_fgt
use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
implicit none
private
public :: forced_step
contains
subroutine forced_step(k, state, force0, force1, diffusion, step, tolerance, next)
integer, intent(in) :: k
real(rk), intent(in) :: state(k,k), force0(k,k), force1(k,k)
real(rk), intent(in) :: diffusion, step, tolerance
real(rk), intent(out) :: next(k,k)
real(rk) :: work(k,k)
if (.not. ieee_is_finite(step) .or. step <= 0.0_rk) error stop "invalid step"
if (.not. all(ieee_is_finite(force0)) .or. .not. all(ieee_is_finite(force1))) error stop "nonfinite forcing"
work = state + 0.5_rk*step*force0
call heat_step_fgt(k, work, diffusion, step, tolerance, next)
next = next + 0.5_rk*step*force1
if (.not. all(ieee_is_finite(next))) error stop "nonfinite forced output"
end subroutine forced_step
end module heat_forced
4. 三行实现对应三个操作
| 实现 | 数学操作 |
|---|---|
work = state + 0.5_rk*step*force0 |
形成 \(\theta^n+\tfrac h2f^n\) |
call heat_step_fgt(...,work,...,next) |
对整个组合场施加 \(H(h)\) |
next = next + 0.5_rk*step*force1 |
加入 \(\tfrac h2f^{n+1}\) |
左端源项需要经历整个时间间隔的扩散，右端源项在终点加入。因此不能把两端源值都放到传播器外面相加。
work 是局部临时数组，state 不被覆盖；调用成功后由外部把 next 提交为新状态。三组输入数组使用相同节点与形状。
Chapter 6：自适应四叉树与周期 2:1 平衡
1. 目标
用多个不同大小的盒子表示单位域中的场：每个叶盒仍保存 \(k\times k\) 个 Chebyshev 节点值，误差过大的区域继续细分。
build_field 根据函数回调生成树，evaluate_field 在树上定位并插值，mark_unbalanced 标记破坏 2:1 平衡的粗盒。启用 balance=.true. 时，建树同时补齐水平方向周期接缝处的平衡。
2. 完整代码
module adaptive_boxes
use cheb_box, only: rk, cheb_nodes, evaluate_box
use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
implicit none
private
public :: adaptive_field, build_field, evaluate_field, mark_unbalanced
type :: adaptive_field
integer :: k = 0, nbox = 0
real(rk), allocatable :: center(:,:), width(:)
real(rk), allocatable :: values(:,:,:), indicator(:)
integer, allocatable :: level(:), parent(:), child(:,:)
end type adaptive_field
abstract interface
function scalar_field(x, z) result(f)
import rk
real(rk), intent(in) :: x, z
real(rk) :: f
end function scalar_field
end interface
contains
subroutine build_field(fun, k, tolerance, max_level, max_boxes, field, balance)
procedure(scalar_field) :: fun
integer, intent(in) :: k, max_level, max_boxes
real(rk), intent(in) :: tolerance
type(adaptive_field), intent(out) :: field
logical, intent(in), optional :: balance
logical :: enforce_balance, marked(max_boxes)
real(rk) :: nodes(k), xs(k), zs(k)
real(rk) :: xp(2*k+1), zp(2*k+1)
real(rk) :: check_values(2*k+1,2*k+1)
real(rk) :: truth(2*k+1,2*k+1), r(2*k+1), offset(2)
integer :: b, i, j, c, id, m
if (k < 2 .or. max_level < 0 .or. max_boxes < 1) error stop "invalid tree settings"
if (.not. ieee_is_finite(tolerance) .or. tolerance <= 0.0_rk) error stop "invalid tolerance"
enforce_balance = .false.
if (present(balance)) enforce_balance = balance
marked = .false.
field%k = k
field%nbox = 1
allocate(field%center(2,max_boxes), field%width(max_boxes))
allocate(field%values(k,k,max_boxes), field%indicator(max_boxes))
allocate(field%level(max_boxes), field%parent(max_boxes))
allocate(field%child(4,max_boxes))
field%center = 0.0_rk
field%width = 0.0_rk
field%values = 0.0_rk
field%indicator = 0.0_rk
field%level = 0
field%parent = 0
field%child = 0
field%center(:,1) = [0.5_rk,0.5_rk]
field%width(1) = 1.0_rk
call cheb_nodes(k, nodes)
m = 2*k+1
do i = 1, m
r(i) = -1.0_rk + 2.0_rk*real(i-1,rk)/real(m-1,rk)
end do
b = 1
do
do while (b <= field%nbox)
if (field%child(1,b) /= 0) then
b = b+1
cycle
end if
xs = field%center(1,b) + field%width(b)*nodes/2.0_rk
zs = field%center(2,b) + field%width(b)*nodes/2.0_rk
do j = 1, k
do i = 1, k
field%values(i,j,b) = fun(xs(i), zs(j))
end do
end do
xp = field%center(1,b) + field%width(b)*r/2.0_rk
zp = field%center(2,b) + field%width(b)*r/2.0_rk
call evaluate_box(field%values(:,:,b), field%center(:,b), field%width(b), xp, zp, check_values)
do j = 1, m
do i = 1, m
truth(i,j) = fun(xp(i), zp(j))
end do
end do
if (.not. all(ieee_is_finite(check_values)) .or. .not. all(ieee_is_finite(truth))) error stop "nonfinite field"
field%indicator(b) = maxval(abs(check_values-truth))
if (field%indicator(b) > tolerance .or. marked(b)) then
if (field%level(b) >= max_level) error stop "maximum level reached"
if (field%nbox+4 > max_boxes) error stop "maximum box count reached"
do c = 1, 4
id = field%nbox+c
offset = [real(2*mod(c-1,2)-1,rk), real(2*((c-1)/2)-1,rk)]
field%center(:,id) = field%center(:,b) + field%width(b)*offset/4.0_rk
field%width(id) = field%width(b)/2.0_rk
field%level(id) = field%level(b)+1
field%parent(id) = b
field%child(c,b) = id
end do
field%nbox = field%nbox+4
end if
b = b+1
end do
if (.not. enforce_balance) exit
call mark_unbalanced(field, marked)
if (.not. any(marked(1:field%nbox))) exit
b = 1
end do
end subroutine build_field
subroutine evaluate_field(field, x, z, u)
type(adaptive_field), intent(in) :: field
real(rk), intent(in) :: x(:), z(:)
real(rk), intent(out) :: u(size(x),size(z))
real(rk) :: value(1,1)
integer :: i, j, b, c
if (field%nbox < 1) error stop "empty field"
if (.not. all(ieee_is_finite(x)) .or. .not. all(ieee_is_finite(z))) error stop "nonfinite target"
if (any(x < 0.0_rk) .or. any(x > 1.0_rk) .or. any(z < 0.0_rk) .or. any(z > 1.0_rk)) error stop "target outside unit box"
do j = 1, size(z)
do i = 1, size(x)
b = 1
do while (field%child(1,b) /= 0)
c = 1
if (x(i) >= field%center(1,b)) c = c+1
if (z(j) >= field%center(2,b)) c = c+2
b = field%child(c,b)
end do
call evaluate_box(field%values(:,:,b), field%center(:,b), field%width(b), [x(i)], [z(j)], value)
u(i,j) = value(1,1)
end do
end do
end subroutine evaluate_field
subroutine mark_unbalanced(field, marked)
type(adaptive_field), intent(in) :: field
logical, intent(out) :: marked(:)
real(rk), parameter :: geometry_tol = 32.0_rk*epsilon(1.0_rk)
real(rk) :: distance(2), reach
integer :: a, b
marked = .false.
do a = 1, field%nbox
if (field%child(1,a) /= 0) cycle
do b = a+1, field%nbox
if (field%child(1,b) /= 0) cycle
if (abs(field%level(a)-field%level(b)) <= 1) cycle
distance = abs(field%center(:,a)-field%center(:,b))
distance(1) = min(distance(1), 1.0_rk-distance(1))
reach = (field%width(a)+field%width(b))/2.0_rk
if (any(distance > reach+geometry_tol)) cycle
if (field%level(a) < field%level(b)) then
marked(a) = .true.
else
marked(b) = .true.
end if
end do
end do
end subroutine mark_unbalanced
end module adaptive_boxes
Chapter 7：将四叉树转换为 FGT 数据布局
1. 目标
由于 FGT 接口要求同层盒子连续排列，并用一个整数数组保存拓扑。
故这里 pack_tree 完成三件事：按层重新编号、转换父子关系、将叶盒值按库的节点顺序打包。它不改变盒子几何，也不负责误差细化或平衡。
2. 完整代码
module fgt_tree_layout
use cheb_box, only: rk
use adaptive_boxes, only: adaptive_field
implicit none
private
public :: fgt_layout, pack_tree
type :: fgt_layout
integer :: k = 0, nbox = 0, nlevels = 0, npbox = 0, ltree = 0
integer :: ptr(8) = 0
integer, allocatable :: tree(:), new_to_old(:), old_to_new(:)
real(rk), allocatable :: centers(:,:), boxsize(:), density(:,:,:)
end type fgt_layout
contains
subroutine pack_tree(field, packed)
type(adaptive_field), intent(in) :: field
type(fgt_layout), intent(out) :: packed
integer :: nb, lev, first, b, old, c, base
if (field%nbox < 1 .or. field%k < 2 .or. field%k > 32) error stop "invalid field for FGT"
nb = field%nbox
packed%k = field%k
packed%nbox = nb
packed%nlevels = maxval(field%level(1:nb))
packed%npbox = field%k**2
packed%ptr(1) = 1
packed%ptr(2) = packed%ptr(1)+2*(packed%nlevels+1)
packed%ptr(3) = packed%ptr(2)+nb
packed%ptr(4) = packed%ptr(3)+nb
packed%ptr(5) = packed%ptr(4)+nb
packed%ptr(6) = packed%ptr(5)+4*nb
packed%ptr(7) = packed%ptr(6)+nb
packed%ptr(8) = packed%ptr(7)+9*nb
packed%ltree = packed%ptr(8)
allocate(packed%tree(packed%ltree))
allocate(packed%new_to_old(nb), packed%old_to_new(nb))
allocate(packed%centers(2,nb), packed%boxsize(0:packed%nlevels))
allocate(packed%density(1,packed%npbox,nb))
packed%tree = -1
packed%density = 0.0_rk
packed%tree(packed%ptr(6):packed%ptr(6)+nb-1) = 0
b = 0
do lev = 0, packed%nlevels
first = b+1
do old = 1, nb
if (field%level(old) /= lev) cycle
b = b+1
packed%new_to_old(b) = old
packed%old_to_new(old) = b
end do
if (b < first) error stop "missing tree level"
base = packed%ptr(1)+2*lev
packed%tree(base:base+1) = [first,b]
packed%boxsize(lev) = field%width(1)/2.0_rk**lev
end do
if (b /= nb) error stop "invalid box levels"
do b = 1, nb
old = packed%new_to_old(b)
packed%centers(:,b) = field%center(:,old)
packed%tree(packed%ptr(2)+b-1) = field%level(old)
if (field%parent(old) /= 0) then
packed%tree(packed%ptr(3)+b-1) = packed%old_to_new(field%parent(old))
end if
packed%tree(packed%ptr(4)+b-1) = 0
if (field%child(1,old) == 0) then
packed%density(1,:,b) = reshape(field%values(field%k:1:-1,field%k:1:-1,old), [packed%npbox])
else
packed%tree(packed%ptr(4)+b-1) = 4
base = packed%ptr(5)+4*(b-1)
do c = 1, 4
packed%tree(base+c-1) = packed%old_to_new(field%child(c,old))
end do
end if
end do
end subroutine pack_tree
end module fgt_tree_layout
Chapter 8：自适应树上的 Gaussian 积分
1. 目标
把已有的 adaptive_field 交给 FGT，并将结果恢复到我们的树编号与节点顺序。
gaussian_tree 返回形状为 (k,k,nbox) 的势值数组。默认计算自由空间 Gaussian 积分；periodic=.true. 时，计算周期长度等于根盒边长的周期积分。本层仍不做热核归一化。
2. 完整代码
module fgt_adaptive
use cheb_box, only: rk
use adaptive_boxes, only: adaptive_field
use fgt_tree_layout, only: fgt_layout, pack_tree
use fgt_single_box, only: boxfgt
use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
implicit none
private
public :: gaussian_tree
contains
subroutine gaussian_tree(field, delta, tolerance, potential, periodic)
type(adaptive_field), intent(in) :: field
real(rk), intent(in) :: delta, tolerance
real(rk), allocatable, intent(out) :: potential(:,:,:)
logical, intent(in), optional :: periodic
type(fgt_layout), allocatable :: packed
integer :: k, np, nb, nlev, ltree, iperiod, b, old
real(rk) :: dd, eps, values(field%k,field%k)
real(rk), allocatable :: p(:,:,:), grad(:,:,:,:), hess(:,:,:,:)
real(rk) :: target(2,1), pt(1,1), gt(1,2,1), ht(1,3,1), times(100)
if (.not. all(ieee_is_finite([delta,tolerance]))) error stop "nonfinite FGT parameter"
if (delta <= 0.0_rk) error stop "delta must be positive"
if (tolerance < 1.0e-14_rk .or. tolerance >= 0.1_rk) error stop "invalid FGT tolerance"
allocate(packed)
call pack_tree(field, packed)
if (.not. all(ieee_is_finite(packed%density))) error stop "nonfinite density"
k = packed%k
np = packed%npbox
nb = packed%nbox
nlev = packed%nlevels
ltree = packed%ltree
iperiod = 0
if (present(periodic)) then
if (periodic) iperiod = 1
end if
dd = delta
eps = tolerance
allocate(p(1,np,nb), grad(1,2,np,nb), hess(1,3,np,nb))
allocate(potential(k,k,nb))
p = 0.0_rk
grad = 0.0_rk
hess = 0.0_rk
potential = 0.0_rk
target = 0.0_rk
pt = 0.0_rk
gt = 0.0_rk
ht = 0.0_rk
times = 0.0_rk
call boxfgt(1, 2, dd, eps, 1, iperiod, k, np, nb, nlev, ltree, packed%tree, packed%ptr, packed%centers, packed%boxsize, packed%density, 1, p, grad, hess, 0, 0, target, 0, pt, gt, ht, times)
if (nb /= field%nbox .or. nlev /= packed%nlevels .or. ltree /= packed%ltree) error stop "unexpected tree change"
do b = 1, nb
old = packed%new_to_old(b)
if (field%child(1,old) /= 0) cycle
if (.not. all(ieee_is_finite(p(1,:,b)))) error stop "nonfinite FGT output"
values = reshape(p(1,:,b), [k,k])
potential(:,:,old) = values(k:1:-1,k:1:-1)
end do
end subroutine gaussian_tree
end module fgt_adaptive
Chapter 9：自适应树上的混合边界热传播
1. 目标
将固定四叶盒的混合边界延拓推广到任意有效的物理自适应树。
extend_mixed 创建四份平移或反射的子树；heat_step_adaptive 调用周期 FGT，取回物理域中的叶节点结果，并完成热核归一化。
2. 完整代码
module heat_adaptive_mixed
use cheb_box, only: rk
use adaptive_boxes, only: adaptive_field
use fgt_adaptive, only: gaussian_tree
use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
implicit none
private
public :: extend_mixed, heat_step_adaptive
contains
subroutine extend_mixed(field, extended)
type(adaptive_field), intent(in) :: field
type(adaptive_field), intent(out) :: extended
integer, parameter :: reflected_child(4) = [3,4,1,2]
integer :: nb, k, q, b, c, old_child, offset, id
if (field%nbox < 1) error stop "empty input tree"
if (maxval(abs(field%center(:,1)-0.5_rk)) > 1.0e-12_rk .or. abs(field%width(1)-1.0_rk) > 1.0e-12_rk) &
error stop "expected a unit physical box"
nb = field%nbox
k = field%k
extended%k = k
extended%nbox = 1+4*nb
allocate(extended%center(2,extended%nbox), extended%width(extended%nbox))
allocate(extended%values(k,k,extended%nbox), extended%indicator(extended%nbox))
allocate(extended%level(extended%nbox), extended%parent(extended%nbox))
allocate(extended%child(4,extended%nbox))
extended%center = 0.0_rk
extended%width = 0.0_rk
extended%values = 0.0_rk
extended%indicator = 0.0_rk
extended%level = 0
extended%parent = 0
extended%child = 0
extended%center(:,1) = [1.0_rk,1.0_rk]
extended%width(1) = 2.0_rk
do q = 1, 4
offset = 1+(q-1)*nb
extended%child(q,1) = offset+1
do b = 1, nb
id = offset+b
extended%center(:,id) = field%center(:,b)
extended%center(1,id) = field%center(1,b)+real(mod(q-1,2),rk)
extended%width(id) = field%width(b)
extended%level(id) = field%level(b)+1
extended%indicator(id) = field%indicator(b)
extended%values(:,:,id) = field%values(:,:,b)
if (b == 1) then
extended%parent(id) = 1
else
extended%parent(id) = offset+field%parent(b)
end if
if (q >= 3) then
extended%center(2,id) = 2.0_rk-field%center(2,b)
extended%values(:,:,id) = -field%values(:,k:1:-1,b)
end if
if (field%child(1,b) == 0) cycle
do c = 1, 4
old_child = c
if (q >= 3) old_child = reflected_child(c)
extended%child(c,id) = offset+field%child(old_child,b)
end do
end do
end do
end subroutine extend_mixed
subroutine heat_step_adaptive(field, diffusion, step, tolerance, potential)
type(adaptive_field), intent(in) :: field
real(rk), intent(in) :: diffusion, step, tolerance
real(rk), allocatable, intent(out) :: potential(:,:,:)
type(adaptive_field), allocatable :: extended
real(rk), allocatable :: all_potential(:,:,:)
real(rk), parameter :: pi = acos(-1.0_rk)
real(rk) :: delta
if (.not. all(ieee_is_finite([diffusion,step]))) error stop "nonfinite heat parameter"
if (diffusion <= 0.0_rk .or. step <= 0.0_rk) error stop "diffusion and step must be positive"
delta = 4.0_rk*diffusion*step
allocate(extended)
call extend_mixed(field, extended)
call gaussian_tree(extended, delta, tolerance, all_potential, periodic=.true.)
allocate(potential(field%k,field%k,field%nbox))
potential = all_potential(:,:,2:field%nbox+1)/(pi*delta)
end subroutine heat_step_adaptive
end module heat_adaptive_mixed
Chapter 10：输出分辨率检查与重新传播
1. 目标
热传播后的场可能进入原先很粗的区域，即使源场已被准确表示，原树也未必能准确表示输出。
checked_heat_step 比较父盒输出插值与细化节点上的重新传播结果。超出容差时细化正式树、恢复平衡，再从同一个数值源重新传播。
2. 完整代码
module heat_output_control
use cheb_box, only: rk, cheb_nodes, evaluate_box
use adaptive_boxes, only: adaptive_field, mark_unbalanced
use heat_adaptive_mixed, only: heat_step_adaptive
use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
implicit none
private
public :: refine_leaves, checked_heat_step
contains
subroutine refine_leaves(field, marked, max_level)
type(adaptive_field), intent(inout) :: field
logical, intent(in) :: marked(:)
integer, intent(in) :: max_level
real(rk) :: nodes(field%k), x(field%k), z(field%k), offset(2)
integer :: b, c, id, old_count
old_count = field%nbox
if (size(marked) < old_count) error stop "mark array too short"
call cheb_nodes(field%k, nodes)
do b = 1, old_count
if (.not. marked(b) .or. field%child(1,b) /= 0) cycle
if (field%level(b) >= max_level) error stop "maximum level reached"
if (field%nbox+4 > size(field%width)) error stop "maximum box count reached"
do c = 1, 4
id = field%nbox+c
offset = [real(2*mod(c-1,2)-1,rk), real(2*((c-1)/2)-1,rk)]
field%center(:,id) = field%center(:,b)+field%width(b)*offset/4.0_rk
field%width(id) = field%width(b)/2.0_rk
field%level(id) = field%level(b)+1
field%parent(id) = b
field%child(:,id) = 0
field%child(c,b) = id
x = field%center(1,id)+field%width(id)*nodes/2.0_rk
z = field%center(2,id)+field%width(id)*nodes/2.0_rk
call evaluate_box(field%values(:,:,b), field%center(:,b), field%width(b), x, z, field%values(:,:,id))
field%indicator(id) = -1.0_rk
end do
field%indicator(b) = -1.0_rk
field%nbox = field%nbox+4
end do
end subroutine refine_leaves
subroutine checked_heat_step(source, diffusion, step, fgt_eps, output_tol, &
max_level, max_passes, output, passes, first_error)
type(adaptive_field), intent(in) :: source
real(rk), intent(in) :: diffusion, step, fgt_eps, output_tol
integer, intent(in) :: max_level, max_passes
type(adaptive_field), intent(out) :: output
integer, intent(out) :: passes
real(rk), intent(out) :: first_error
type(adaptive_field) :: work, probe
real(rk), allocatable :: coarse(:,:,:), fine(:,:,:), eta(:)
logical, allocatable :: marked(:)
real(rk) :: nodes(source%k), x(source%k), z(source%k)
real(rk) :: predicted(source%k,source%k)
integer :: b, c, id
if (.not. ieee_is_finite(output_tol) .or. output_tol <= 0.0_rk) error stop "invalid output tolerance"
if (max_passes < 1 .or. max_level < 0) error stop "invalid refinement limit"
work = source
allocate(marked(size(work%width)), eta(size(work%width)))
call mark_unbalanced(work, marked)
if (any(marked)) error stop "input tree must be balanced"
call cheb_nodes(source%k, nodes)
passes = 0
first_error = -1.0_rk
do
passes = passes+1
call heat_step_adaptive(work, diffusion, step, fgt_eps, coarse)
probe = work
marked = .false.
marked(1:work%nbox) = work%child(1,1:work%nbox) == 0
call refine_leaves(probe, marked, max_level+1)
call heat_step_adaptive(probe, diffusion, step, fgt_eps, fine)
eta = 0.0_rk
do b = 1, work%nbox
if (work%child(1,b) /= 0) cycle
do c = 1, 4
id = probe%child(c,b)
x = probe%center(1,id)+probe%width(id)*nodes/2.0_rk
z = probe%center(2,id)+probe%width(id)*nodes/2.0_rk
call evaluate_box(coarse(:,:,b), work%center(:,b), work%width(b), x, z, predicted)
eta(b) = max(eta(b), maxval(abs(predicted-fine(:,:,id))))
end do
end do
if (.not. all(ieee_is_finite(eta))) error stop "nonfinite output indicator"
if (passes == 1) first_error = maxval(eta)
if (maxval(eta) <= output_tol) then
output = work
output%values(:,:,1:work%nbox) = coarse
output%indicator = -1.0_rk
do b = 1, work%nbox
if (work%child(1,b) == 0) output%indicator(b) = eta(b)
end do
return
end if
if (passes >= max_passes) error stop "output resolution did not converge"
marked = eta > output_tol
call refine_leaves(work, marked, max_level)
do
call mark_unbalanced(work, marked)
if (.not. any(marked)) exit
call refine_leaves(work, marked, max_level)
end do
end do
end subroutine checked_heat_step
end module heat_output_control
Chapter 11：受控粗化与重网格
1. 目标
扩散使场逐渐平滑后，部分细盒不再必要。coarsen_field 尝试用父盒多项式替代原子树，并在控制场转移误差的同时保持周期 2:1 平衡。
粗化不推进时间。输入 reference 和输出 result 表示同一时刻的场，且必须使用不同变量。
2. 完整代码
module field_regrid
use cheb_box, only: rk, cheb_nodes, evaluate_box
use adaptive_boxes, only: adaptive_field, evaluate_field, mark_unbalanced
use heat_output_control, only: refine_leaves
use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
implicit none
private
public :: coarsen_field
contains
subroutine transfer_indicator(reference, center, width, values, error)
type(adaptive_field), intent(in) :: reference
real(rk), intent(in) :: center(2), width, values(:,:)
real(rk), intent(out) :: error
integer :: a, i, m
real(rk) :: r(2*reference%k+1), x(2*reference%k+1), z(2*reference%k+1)
real(rk) :: old_values(2*reference%k+1,2*reference%k+1)
real(rk) :: new_values(2*reference%k+1,2*reference%k+1)
real(rk) :: lower(2), upper(2)
logical :: covered
m = size(r)
r = [(real(i-1,rk)/real(m-1,rk), i=1,m)]
error = 0.0_rk
covered = .false.
do a = 1, reference%nbox
if (reference%child(1,a) /= 0) cycle
lower = max(center-width/2.0_rk, &
reference%center(:,a)-reference%width(a)/2.0_rk)
upper = min(center+width/2.0_rk, &
reference%center(:,a)+reference%width(a)/2.0_rk)
if (any(upper <= lower)) cycle
covered = .true.
x = lower(1)+(upper(1)-lower(1))*r
z = lower(2)+(upper(2)-lower(2))*r
call evaluate_box(reference%values(:,:,a), reference%center(:,a), &
reference%width(a), x, z, old_values)
call evaluate_box(values, center, width, x, z, new_values)
if (.not. all(ieee_is_finite(old_values)) .or. &
.not. all(ieee_is_finite(new_values))) error stop "nonfinite transfer"
error = max(error, maxval(abs(new_values-old_values)))
end do
if (.not. covered) error stop "box outside reference tree"
end subroutine transfer_indicator
subroutine coarsen_field(reference, tolerance, result, monitor)
type(adaptive_field), intent(in) :: reference
real(rk), intent(in) :: tolerance
type(adaptive_field), intent(out) :: result
real(rk), intent(out) :: monitor
logical, allocatable :: marked(:)
real(rk) :: nodes(reference%k), error
integer :: b, max_level
if (.not. ieee_is_finite(tolerance) .or. tolerance <= 0.0_rk) &
error stop "invalid transfer tolerance"
if (reference%nbox < 1) error stop "empty reference tree"
allocate(marked(size(reference%width)))
call mark_unbalanced(reference, marked)
if (any(marked)) error stop "reference tree must be balanced"
max_level = maxval(reference%level(1:reference%nbox))
call cheb_nodes(reference%k, nodes)
result = reference
result%nbox = 1
result%center = 0.0_rk
result%width = 0.0_rk
result%values = 0.0_rk
result%indicator = -1.0_rk
result%level = 0
result%parent = 0
result%child = 0
call visit(1,1)
do
call mark_unbalanced(result, marked)
if (.not. any(marked)) exit
call refine_leaves(result, marked, max_level)
end do
monitor = 0.0_rk
do b = 1, result%nbox
if (result%child(1,b) /= 0) cycle
call transfer_indicator(reference, result%center(:,b), &
result%width(b), result%values(:,:,b), error)
monitor = max(monitor,error)
end do
if (monitor > tolerance) error stop "final transfer check failed"
result%indicator = -1.0_rk
contains
recursive subroutine visit(old, id)
integer, intent(in) :: old, id
real(rk) :: x(reference%k), z(reference%k)
real(rk) :: candidate(reference%k,reference%k), local_error
integer :: c, first
result%center(:,id) = reference%center(:,old)
result%width(id) = reference%width(old)
result%level(id) = reference%level(old)
if (reference%child(1,old) == 0) then
result%values(:,:,id) = reference%values(:,:,old)
return
end if
x = reference%center(1,old)+reference%width(old)*nodes/2.0_rk
z = reference%center(2,old)+reference%width(old)*nodes/2.0_rk
call evaluate_field(reference, x, z, candidate)
call transfer_indicator(reference, reference%center(:,old), &
reference%width(old), candidate, local_error)
if (local_error <= tolerance) then
result%values(:,:,id) = candidate
return
end if
if (result%nbox+4 > size(result%width)) error stop "maximum box count reached"
first = result%nbox+1
result%nbox = result%nbox+4
do c = 1, 4
result%child(c,id) = first+c-1
result%parent(first+c-1) = id
end do
do c = 1, 4
call visit(reference%child(c,old),first+c-1)
end do
end subroutine visit
end subroutine coarsen_field
end module field_regrid
Chapter 12：连续自适应热推进
1. 目标
把已有的热传播、输出细化和粗化串成时间步。现在传入的是上一时刻的数值场，以及可以在任意位置、任意时刻计算的已知源项。
2. 完整代码
module heat_time_driver
use cheb_box, only: rk
use adaptive_boxes, only: adaptive_field, build_field, evaluate_field
use heat_output_control, only: checked_heat_step
use field_regrid, only: coarsen_field
use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
implicit none
private
public :: advance_known_heat
abstract interface
real(kind=kind(1.0d0)) function known_source(x,z,t)
import rk
real(rk), intent(in) :: x,z,t
end function known_source
end interface
contains
subroutine advance_known_heat(state,time,h,diffusion,forcing,space_tol,next,passes,transfer_error)
type(adaptive_field), intent(in) :: state
real(rk), intent(in) :: time,h,diffusion,space_tol
procedure(known_source) :: forcing
type(adaptive_field), intent(out) :: next
integer, intent(out) :: passes
real(rk), intent(out) :: transfer_error
type(adaptive_field) :: combined, evolved
real(rk) :: first_error
if (h <= 0 .or. diffusion <= 0) error stop "positive timestep and diffusion required"
call build_field(combined_value,state%k,space_tol,9,size(state%width),combined,.true.)
call checked_heat_step(combined,diffusion,h,1.0e-12_rk,space_tol,9,6,evolved,passes,first_error)
! Rebuild after adding the endpoint source, which can contain new spatial structure.
call build_field(endpoint_value,state%k,space_tol,9,size(state%width),combined,.true.)
call coarsen_field(combined,space_tol,next,transfer_error)
if (.not. all(ieee_is_finite(next%values(:,:,1:next%nbox)))) error stop "invalid heat state"
contains
real(rk) function combined_value(x,z) result(v)
real(rk), intent(in) :: x,z
real(rk) :: old(1,1)
call evaluate_field(state,[x],[z],old)
v = old(1,1)+h*forcing(x,z,time)/2
end function combined_value
real(rk) function endpoint_value(x,z) result(v)
real(rk), intent(in) :: x,z
real(rk) :: old(1,1)
call evaluate_field(evolved,[x],[z],old)
v = old(1,1)+h*forcing(x,z,time+h)/2
end function endpoint_value
end subroutine advance_known_heat
end module heat_time_driver
Chapter 13：小型数值工具
1. 目标
后面需要解小型椭圆系统、计算积分和导数。把这些重复操作集中到一个模块，主算法只保留物理步骤。
2. 完整代码
module numerical_tools
use cheb_box, only: rk
implicit none
private
public :: pi, inverse_matrix, gauss_rule, interpolation_matrix, differentiation_matrix
real(rk), parameter :: pi = acos(-1.0_rk)
contains
function inverse_matrix(a) result(b)
real(rk), intent(in) :: a(:,:)
real(rk) :: b(size(a,1),size(a,1)), work(size(a,1),size(a,1)), row(size(a,1)), pivot, factor
integer :: n, i, j, p
n = size(a,1)
if (size(a,2) /= n) error stop "inverse requires a square matrix"
work = a
b = 0
do i = 1, n
b(i,i) = 1
end do
do j = 1, n
p = j-1+maxloc(abs(work(j:n,j)),dim=1)
if (abs(work(p,j)) < tiny(1.0_rk)) error stop "singular matrix"
row = work(j,:)
work(j,:) = work(p,:)
work(p,:) = row
row = b(j,:)
b(j,:) = b(p,:)
b(p,:) = row
pivot = work(j,j)
work(j,:) = work(j,:)/pivot
b(j,:) = b(j,:)/pivot
do i = 1, n
if (i == j) cycle
factor = work(i,j)
work(i,:) = work(i,:)-factor*work(j,:)
b(i,:) = b(i,:)-factor*b(j,:)
end do
end do
end function inverse_matrix
subroutine gauss_rule(n,x,w)
integer, intent(in) :: n
real(rk), intent(out) :: x(n), w(n)
real(rk) :: z, previous, p0, p1, p2, derivative
integer :: i, j, it
do i = 1, (n+1)/2
z = cos(pi*(real(i,rk)-0.25_rk)/(real(n,rk)+0.5_rk))
do it = 1, 100
p0 = 1
p1 = z
do j = 2, n
p2 = ((2*j-1)*z*p1-(j-1)*p0)/j
p0 = p1
p1 = p2
end do
derivative = n*(z*p1-p0)/(z*z-1)
previous = z
z = z-p1/derivative
if (abs(z-previous) < 4*epsilon(z)) exit
end do
if (it > 100) error stop "Gauss rule iteration failed"
x(i) = -z
x(n+1-i) = z
w(i) = 2/((1-z*z)*derivative**2)
w(n+1-i) = w(i)
end do
end subroutine gauss_rule
function interpolation_matrix(nodes,weights,targets) result(b)
real(rk), intent(in) :: nodes(:), weights(:), targets(:)
real(rk) :: b(size(targets),size(nodes)), difference(size(nodes))
integer :: i, j
do i = 1, size(targets)
difference = targets(i)-nodes
j = minloc(abs(difference),dim=1)
if (abs(difference(j)) < 8*epsilon(1.0_rk)) then
b(i,:) = 0
b(i,j) = 1
else
b(i,:) = weights/difference
b(i,:) = b(i,:)/sum(b(i,:))
end if
end do
end function interpolation_matrix
function differentiation_matrix(nodes,weights) result(d)
real(rk), intent(in) :: nodes(:), weights(:)
real(rk) :: d(size(nodes),size(nodes))
integer :: i, j
d = 0
do i = 1, size(nodes)
do j = 1, size(nodes)
if (i /= j) d(i,j) = weights(j)/(weights(i)*(nodes(i)-nodes(j)))
end do
d(i,i) = -sum(d(i,:))
end do
end function differentiation_matrix
end module numerical_tools
Chapter 14：周期 Poisson、Hodge 分解与周期 Stokes
1. 目标
先在两个方向都周期的单位方盒上建立椭圆逆与无散分解。该模块接收前面的自适应树，Poisson 逆由真实 FGT 热传播积分组成。
2. 完整代码
module periodic_operators
use cheb_box, only: rk, cheb_nodes
use adaptive_boxes, only: adaptive_field
use heat_output_control, only: refine_leaves
use fgt_adaptive, only: gaussian_tree
use numerical_tools, only: pi, differentiation_matrix, gauss_rule
implicit none
private
public :: periodic_heat, periodic_poisson, hodge_project, periodic_stokes_step, field_derivative, field_mean, balance_periodic
contains
real(rk) function field_mean(field) result(mean)
type(adaptive_field), intent(in) :: field
real(rk) :: weights(field%k),angle(field%k)
integer :: i,n,b
angle = [(pi*(i-0.5_rk)/field%k,i=1,field%k)]
weights = 1
do n = 2, field%k-1, 2
weights = weights+2*cos(n*angle)/(1.0_rk-n*n)
end do
weights = weights/field%k
mean = 0
do b = 1, field%nbox
if (field%child(1,b) /= 0) cycle
mean = mean+field%width(b)**2*dot_product(weights,matmul(field%values(:,:,b),weights))
end do
mean = mean/field%width(1)**2
end function field_mean
function field_derivative(field,axis) result(derivative)
type(adaptive_field), intent(in) :: field
integer, intent(in) :: axis
type(adaptive_field) :: derivative
real(rk) :: nodes(field%k),weights(field%k),d(field%k,field%k)
integer :: j,b
if (axis < 1 .or. axis > 2) error stop "derivative axis must be 1 or 2"
call cheb_nodes(field%k,nodes)
weights = [((-1.0_rk)**(j-1)*sin(pi*(j-0.5_rk)/field%k),j=1,field%k)]
d = differentiation_matrix(nodes,weights)
derivative = field
derivative%values = 0
derivative%indicator = -1
do b = 1, field%nbox
if (field%child(1,b) /= 0) cycle
if (axis == 1) then
derivative%values(:,:,b) = 2*matmul(d,field%values(:,:,b))/field%width(b)
else
derivative%values(:,:,b) = 2*matmul(field%values(:,:,b),transpose(d))/field%width(b)
end if
end do
end function field_derivative
subroutine balance_periodic(field)
type(adaptive_field), intent(inout) :: field
logical :: marked(size(field%width))
real(rk) :: distance(2),reach
integer :: a,b
do
marked = .false.
do a = 1, field%nbox
if (field%child(1,a) /= 0) cycle
do b = a+1, field%nbox
if (field%child(1,b) /= 0 .or. abs(field%level(a)-field%level(b)) <= 1) cycle
distance = abs(field%center(:,a)-field%center(:,b))
distance = min(distance,field%width(1)-distance)
reach = (field%width(a)+field%width(b))/2+32*epsilon(1.0_rk)
if (any(distance > reach)) cycle
if (field%level(a) < field%level(b)) then
marked(a) = .true.
else
marked(b) = .true.
end if
end do
end do
if (.not. any(marked)) exit
call refine_leaves(field,marked,10)
end do
end subroutine balance_periodic
function periodic_heat(field,diffusion,h,eps) result(output)
type(adaptive_field), intent(in) :: field
real(rk), intent(in) :: diffusion,h,eps
type(adaptive_field) :: output
real(rk), allocatable :: potential(:,:,:)
if (h <= 0 .or. diffusion <= 0) error stop "positive heat lag required"
output = field
call gaussian_tree(field,4*diffusion*h,eps,potential,.true.)
output%values(:,:,1:field%nbox) = potential/(4*pi*diffusion*h)
output%indicator = -1
end function periodic_heat
subroutine periodic_poisson(rho,tolerance,phi,monitor,heat_calls)
type(adaptive_field), intent(in) :: rho
real(rk), intent(in) :: tolerance
type(adaptive_field), intent(out) :: phi
real(rk), intent(out) :: monitor
integer, intent(out) :: heat_calls
type(adaptive_field) :: source,coarse,fine
real(rk) :: mean,bound,left_tail,right_tail,start,finish,gauge
integer :: panels,b
if (abs(rho%width(1)-1) > epsilon(1.0_rk)) error stop "Poisson requires unit periodic square"
if (tolerance < 1.0e-11_rk .or. tolerance > 1.0e-3_rk) error stop "invalid Poisson tolerance"
mean = field_mean(rho)
if (abs(mean) > tolerance/8) error stop "non-neutral periodic Poisson source"
source = rho
bound = 0
do b = 1, rho%nbox
if (rho%child(1,b) /= 0) cycle
source%values(:,:,b) = rho%values(:,:,b)-mean
bound = max(bound,maxval(abs(source%values(:,:,b))))
end do
heat_calls = 0
monitor = 0
phi = source
phi%values = 0
if (bound < tiny(1.0_rk)) return
start = min(1.0e-5_rk,tolerance/(128*bound))
finish = max(0.25_rk,log(max(1.0_rk,128*bound/tolerance))/(4*pi*pi))
left_tail = start*bound
right_tail = bound*exp(-4*pi*pi*finish)/(pi*pi*(1-exp(-4*pi*pi*finish))**2)
panels = ceiling(log(finish/start)/2)
call integrate_heat(8,coarse)
call integrate_heat(16,fine)
monitor = maxval(abs(coarse%values(:,:,1:rho%nbox)-fine%values(:,:,1:rho%nbox)))+left_tail+right_tail
if (monitor > tolerance/2) error stop "heat-integral quadrature unresolved"
phi = fine
gauge = field_mean(phi)
do b = 1, phi%nbox
if (phi%child(1,b) == 0) phi%values(:,:,b) = phi%values(:,:,b)-gauge
end do
contains
subroutine integrate_heat(order,result)
integer, intent(in) :: order
type(adaptive_field), intent(out) :: result
type(adaptive_field) :: evolved
real(rk) :: nodes(order),weights(order),lo,hi,mid,half,lag,weight
integer :: panel,j
result = source
result%values = 0
call gauss_rule(order,nodes,weights)
do panel = 1, panels
lo = log(start)+log(finish/start)*(panel-1)/panels
hi = log(start)+log(finish/start)*panel/panels
mid = (lo+hi)/2
half = (hi-lo)/2
do j = 1, order
lag = exp(mid+half*nodes(j))
weight = half*weights(j)*lag
evolved = periodic_heat(source,1.0_rk,lag,1.0e-12_rk)
result%values(:,:,1:source%nbox) = result%values(:,:,1:source%nbox)-weight*evolved%values(:,:,1:source%nbox)
heat_calls = heat_calls+1
end do
end do
end subroutine integrate_heat
end subroutine periodic_poisson
subroutine hodge_project(fx,fz,tolerance,ux,uz,pressure,monitor,calls)
type(adaptive_field), intent(in) :: fx,fz
real(rk), intent(in) :: tolerance
type(adaptive_field), intent(out) :: ux,uz,pressure
real(rk), intent(out) :: monitor
integer, intent(out) :: calls
type(adaptive_field) :: gx,gz,rho
if (fx%nbox /= fz%nbox .or. fx%k /= fz%k) error stop "Hodge requires a common tree"
if (any(fx%child(:,1:fx%nbox) /= fz%child(:,1:fz%nbox))) error stop "Hodge tree mismatch"
gx = field_derivative(fx,1)
gz = field_derivative(fz,2)
rho = gx
rho%values = gx%values+gz%values
call periodic_poisson(rho,tolerance,pressure,monitor,calls)
gx = field_derivative(pressure,1)
gz = field_derivative(pressure,2)
ux = fx
uz = fz
ux%values = fx%values-gx%values
uz%values = fz%values-gz%values
end subroutine hodge_project
function periodic_stokes_step(velocity,force0,force1,viscosity,h) result(next)
type(adaptive_field), intent(in) :: velocity,force0,force1
real(rk), intent(in) :: viscosity,h
type(adaptive_field) :: combined,next
if (velocity%nbox /= force0%nbox .or. velocity%nbox /= force1%nbox) error stop "Stokes source tree mismatch"
combined = velocity
combined%values = velocity%values+h*force0%values/2
next = periodic_heat(combined,viscosity,h,1.0e-12_rk)
next%values = next%values+h*force1%values/2
end function periodic_stokes_step
end module periodic_operators
Chapter 15：板间问题的辅助网格
1. 目标
为真实上下壁准备速度求导、Poisson 逆、密输出和去混叠所需的公共表示。水平使用 Fourier，竖直使用包含两壁的 Chebyshev–Lobatto 节点。
2. 完整代码
module slab_grid
use cheb_box, only: rk
use numerical_tools
implicit none
private
public :: auxiliary_grid, make_grid, fourier, physical, dx, dz, laplacian, poisson_dirichlet
public :: grid_coefficients, evaluate_grid, fine_values, restrict_product, integral, inner_product
type :: auxiliary_grid
integer :: nx = 0, nz = 0, nxf = 0, nzf = 0
real(rk), allocatable :: x(:), z(:), frequency(:), d(:,:), d2(:,:), v(:,:), iv(:,:), inverse(:,:,:)
real(rk), allocatable :: lift(:,:,:), xf(:), zf(:), vf(:,:), ivf(:,:), gram(:,:), moments(:)
complex(rk), allocatable :: transform(:,:), synthesis(:,:), fine_synthesis(:,:), fine_transform(:,:)
end type auxiliary_grid
contains
subroutine make_grid(g,nx,nz)
type(auxiliary_grid), intent(out) :: g
integer, intent(in) :: nx, nz
real(rk) :: weights(nz), matrix(nz-2,nz-2), a
integer :: i, j, m, mode
if (nx < 8 .or. mod(nx,2) /= 0 .or. nz < 8) error stop "invalid auxiliary grid"
g%nx = nx
g%nz = nz
g%nxf = 3*nx/2
g%nzf = 2*nz-1
allocate(g%x(nx),g%z(nz),g%frequency(nx),g%d(nz,nz),g%d2(nz,nz),g%v(nz,nz),g%iv(nz,nz))
allocate(g%inverse(nz-2,nz-2,nx),g%lift(nz,2,nx),g%gram(nz,nz),g%moments(nz))
allocate(g%xf(g%nxf),g%zf(g%nzf),g%vf(g%nzf,g%nzf),g%ivf(g%nzf,g%nzf))
allocate(g%transform(nx,nx),g%synthesis(nx,nx),g%fine_synthesis(g%nxf,nx),g%fine_transform(nx,g%nxf))
g%x = [(real(i-1,rk)/nx,i=1,nx)]
g%z = [(0.5_rk*(1-cos(pi*real(i-1,rk)/(nz-1))),i=1,nz)]
weights = [((-1.0_rk)**(i-1),i=1,nz)]
weights(1) = weights(1)/2
weights(nz) = weights(nz)/2
g%d = differentiation_matrix(g%z,weights)
g%d2 = matmul(g%d,g%d)
g%xf = [(real(i-1,rk)/g%nxf,i=1,g%nxf)]
g%zf = [(0.5_rk*(1-cos(pi*real(i-1,rk)/(g%nzf-1))),i=1,g%nzf)]
do j = 1, nz
g%v(:,j) = cos((j-1)*acos(2*g%z-1))
g%moments(j) = moment(j-1)
do i = 1, nz
g%gram(i,j) = (moment(i+j-2)+moment(abs(i-j)))/2
end do
end do
do j = 1, g%nzf
g%vf(:,j) = cos((j-1)*acos(2*g%zf-1))
end do
g%iv = inverse_matrix(g%v)
g%ivf = inverse_matrix(g%vf)
do m = 1, nx
mode = m-1
if (mode >= nx/2) mode = mode-nx
g%frequency(m) = 2*pi*mode
g%synthesis(:,m) = exp(cmplx(0.0_rk,g%frequency(m)*g%x,rk))
g%transform(m,:) = conjg(g%synthesis(:,m))/nx
g%fine_synthesis(:,m) = exp(cmplx(0.0_rk,g%frequency(m)*g%xf,rk))
g%fine_transform(m,:) = conjg(g%fine_synthesis(:,m))/g%nxf
a = abs(g%frequency(m))
matrix = -g%d2(2:nz-1,2:nz-1)
do i = 1, nz-2
matrix(i,i) = matrix(i,i)+a*a
end do
g%inverse(:,:,m) = inverse_matrix(matrix)
if (m == 1) then
g%lift(:,1,m) = 1-g%z
g%lift(:,2,m) = g%z
else
g%lift(:,1,m) = (exp(-a*g%z)-exp(-a*(2-g%z)))/(1-exp(-2*a))
g%lift(:,2,m) = (exp(-a*(1-g%z))-exp(-a*(1+g%z)))/(1-exp(-2*a))
end if
end do
end subroutine make_grid
pure real(rk) function moment(n)
integer, intent(in) :: n
moment = 0
if (mod(n,2) == 0) moment = 1.0_rk/(1.0_rk-real(n,rk)**2)
end function moment
function fourier(g,v) result(c)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: v(:,:)
complex(rk) :: c(g%nx,g%nz)
c = matmul(g%transform,cmplx(v,0.0_rk,rk))
end function fourier
function physical(g,c) result(v)
type(auxiliary_grid), intent(in) :: g
complex(rk), intent(in) :: c(:,:)
real(rk) :: v(g%nx,g%nz)
v = real(matmul(g%synthesis,c),rk)
end function physical
function dx(g,v) result(result)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: v(:,:)
real(rk) :: result(g%nx,g%nz)
complex(rk) :: c(g%nx,g%nz)
integer :: m
c = fourier(g,v)
do m = 1, g%nx
c(m,:) = cmplx(0.0_rk,g%frequency(m),rk)*c(m,:)
end do
result = physical(g,c)
end function dx
function dz(g,v) result(result)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: v(:,:)
real(rk) :: result(g%nx,g%nz)
result = matmul(v,transpose(g%d))
end function dz
function laplacian(g,v) result(result)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: v(:,:)
real(rk) :: result(g%nx,g%nz)
complex(rk) :: c(g%nx,g%nz)
integer :: m
c = fourier(g,v)
do m = 1, g%nx
c(m,:) = -g%frequency(m)**2*c(m,:)
end do
result = physical(g,c)+matmul(v,transpose(g%d2))
end function laplacian
function poisson_dirichlet(g,rho) result(p)
type(auxiliary_grid), intent(in) :: g
complex(rk), intent(in) :: rho(:,:)
complex(rk) :: p(g%nx,g%nz)
integer :: m
p = 0
do m = 1, g%nx
p(m,2:g%nz-1) = matmul(g%inverse(:,:,m),rho(m,2:g%nz-1))
end do
end function poisson_dirichlet
function grid_coefficients(g,v) result(c)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: v(:,:)
complex(rk) :: c(g%nx,g%nz)
c = matmul(fourier(g,v),transpose(g%iv))
end function grid_coefficients
function evaluate_grid(g,c,x,z) result(v)
type(auxiliary_grid), intent(in) :: g
complex(rk), intent(in) :: c(:,:)
real(rk), intent(in) :: x(:), z(:)
real(rk) :: v(size(x),size(z)), vz(size(z),g%nz)
complex(rk) :: ex(size(x),g%nx)
integer :: j
do j = 1, g%nx
ex(:,j) = exp(cmplx(0.0_rk,g%frequency(j)*x,rk))
end do
do j = 1, g%nz
vz(:,j) = cos((j-1)*acos(max(-1.0_rk,min(1.0_rk,2*z-1))))
end do
v = real(matmul(matmul(ex,c),transpose(vz)),rk)
end function evaluate_grid
function fine_values(g,v) result(f)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: v(:,:)
real(rk) :: f(g%nxf,g%nzf)
f = real(matmul(matmul(g%fine_synthesis,grid_coefficients(g,v)),transpose(g%vf(:,1:g%nz))),rk)
end function fine_values
function restrict_product(g,f) result(v)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: f(:,:)
real(rk) :: v(g%nx,g%nz)
complex(rk) :: c(g%nx,g%nzf)
c = matmul(matmul(g%fine_transform,cmplx(f,0.0_rk,rk)),transpose(g%ivf))
c(g%nx/2+1,:) = 0
v = physical(g,matmul(c(:,1:g%nz),transpose(g%v)))
end function restrict_product
real(rk) function integral(g,v) result(value)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: v(:,:)
value = dot_product(matmul(g%iv,sum(v,dim=1)/g%nx),g%moments)
end function integral
real(rk) function inner_product(g,v,w) result(value)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: v(:,:), w(:,:)
complex(rk) :: a(g%nx,g%nz), b(g%nx,g%nz)
integer :: m
a = grid_coefficients(g,v)
b = grid_coefficients(g,w)
value = 0
do m = 1, g%nx
if (m == g%nx/2+1) then
value = value+real(dot_product(a(m,:),matmul(g%gram,b(m,:))),rk)/2
else
value = value+real(dot_product(a(m,:),matmul(g%gram,b(m,:))),rk)
end if
end do
end function inner_product
end module slab_grid
Chapter 16：辅助表示与自适应 FGT 的连接
1. 目标
让速度涡量、平均流和温度的体积扩散都调用现有自适应混合边界热算子，同时把源表示和输出表示的误差分开检查。
2. 完整代码
module slab_heat
use cheb_box, only: rk, cheb_nodes, evaluate_box
use adaptive_boxes, only: adaptive_field, build_field, evaluate_field, mark_unbalanced
use heat_output_control, only: refine_leaves, checked_heat_step
use slab_grid, only: auxiliary_grid, grid_coefficients, evaluate_grid
implicit none
private
public :: sample_grid_field, heat_on_grid, heat_statistics
type :: heat_statistics
integer :: calls = 0, min_leaves = huge(1), max_leaves = 0, max_level = 0, refinements = 0
real(rk) :: source_monitor = 0, output_monitor = 0
end type heat_statistics
contains
subroutine sample_grid_field(g,values,k,tol,min_level,field)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: values(:,:), tol
integer, intent(in) :: k, min_level
type(adaptive_field), intent(out) :: field
complex(rk) :: coefficients(g%nx,g%nz)
real(rk) :: nodes(k), x(k), z(k), probes(2*k+1), xp(2*k+1), zp(2*k+1), predicted(2*k+1,2*k+1)
logical :: marked(4096), balance_marks(4096)
integer :: b,i,pass
coefficients = grid_coefficients(g,values)
call build_field(zero,k,tol,0,4096,field)
call cheb_nodes(k,nodes)
probes = [(-1.0_rk+real(i-1,rk)/k,i=1,2*k+1)]
do pass = 1, 9
marked = .false.
do b = 1, field%nbox
if (field%child(1,b) /= 0) cycle
x = field%center(1,b)+field%width(b)*nodes/2
z = field%center(2,b)+field%width(b)*nodes/2
field%values(:,:,b) = evaluate_grid(g,coefficients,x,z)
xp = field%center(1,b)+field%width(b)*probes/2
zp = field%center(2,b)+field%width(b)*probes/2
call evaluate_box(field%values(:,:,b),field%center(:,b),field%width(b),xp,zp,predicted)
field%indicator(b) = maxval(abs(predicted-evaluate_grid(g,coefficients,xp,zp)))
marked(b) = field%indicator(b) > tol .or. field%level(b) < min_level
end do
call mark_unbalanced(field,balance_marks)
marked = marked .or. balance_marks
if (.not. any(marked)) return
call refine_leaves(field,marked,7)
end do
error stop "source sampling did not converge"
contains
real(rk) function zero(x,z)
real(rk), intent(in) :: x,z
zero = 0*(x+z)
end function zero
end subroutine sample_grid_field
function heat_on_grid(g,values,diffusion,h,k,tol,min_level,stats) result(output)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: values(:,:), diffusion,h,tol
integer, intent(in) :: k,min_level
type(heat_statistics), intent(inout) :: stats
real(rk) :: output(g%nx,g%nz), first_error
type(adaptive_field) :: source, evolved
integer :: passes, leaves
if (h <= 0 .or. diffusion <= 0) error stop "invalid heat parameters"
if (maxval(abs(values)) < 1.0e-27_rk) then
output = 0
return
end if
call sample_grid_field(g,values,k,tol,min_level,source)
stats%source_monitor = max(stats%source_monitor,maxval(source%indicator(1:source%nbox)))
call checked_heat_step(source,diffusion,h,1.0e-12_rk,8*tol,7,5,evolved,passes,first_error)
call evaluate_field(evolved,g%x,g%z,output)
stats%calls = stats%calls+2*passes
leaves = count(evolved%child(1,1:evolved%nbox) == 0)
stats%min_leaves = min(stats%min_leaves,leaves)
stats%max_leaves = max(stats%max_leaves,leaves)
stats%max_level = max(stats%max_level,maxval(evolved%level(1:evolved%nbox)))
stats%refinements = stats%refinements+passes-1
stats%output_monitor = max(stats%output_monitor,maxval(evolved%indicator(1:evolved%nbox)))
end function heat_on_grid
end module slab_heat
Chapter 17：真实无滑移 Stokes 壁面
1. 目标
同时满足不可压条件和两壁的两个速度分量为零。采用流函数–涡量表示，以两列壁响应修正 FGT 得到的体积热势。
2. 完整代码
module wall_stokes
use cheb_box, only: rk
use numerical_tools, only: pi, inverse_matrix
use slab_grid
implicit none
private
public :: wall_operator, make_wall_operator, apply_lift, impose_noslip, velocity_from_vorticity, recover_pressure
type :: wall_operator
real(rk), allocatable :: j(:,:,:), b(:,:,:), influence_inverse(:,:,:)
end type wall_operator
contains
subroutine make_wall_operator(g,diffusion,h,wall)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: diffusion,h
type(wall_operator), intent(out) :: wall
real(rk) :: heat(g%nz,2), psi(g%nz,2), influence(2,2), rate, coefficient, a
integer :: m,n,c,nmax
if (diffusion <= 0 .or. h <= 0) error stop "invalid wall response parameters"
allocate(wall%j(g%nz,2,g%nx),wall%b(g%nz,2,g%nx),wall%influence_inverse(2,2,g%nx))
wall%j = 0
wall%influence_inverse = 0
nmax = ceiling(sqrt(40/(diffusion*h))/pi)+1
do m = 1, g%nx
a = abs(g%frequency(m))
heat = 0
! Only two boundary lifting functions are propagated analytically.
do n = 1, nmax
rate = a*a+(pi*n)**2
coefficient = 2*pi*n*exp(-diffusion*h*rate)/rate
heat(:,1) = heat(:,1)+coefficient*sin(pi*n*g%z)
heat(:,2) = heat(:,2)-(-1.0_rk)**n*coefficient*sin(pi*n*g%z)
end do
do c = 1, 2
wall%j(2:g%nz-1,c,m) = matmul(g%inverse(:,:,m),g%lift(2:g%nz-1,c,m)-heat(2:g%nz-1,c))/diffusion
end do
wall%b(:,:,m) = g%lift(:,:,m)-wall%j(:,:,m)/h
if (m == 1) cycle
psi = 0
psi(2:g%nz-1,:) = matmul(g%inverse(:,:,m),wall%b(2:g%nz-1,:,m))
influence(1,:) = matmul(g%d(1,:),psi)
influence(2,:) = matmul(g%d(g%nz,:),psi)
wall%influence_inverse(:,:,m) = inverse_matrix(influence)
end do
end subroutine make_wall_operator
function apply_lift(g,matrix,boundary) result(v)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: matrix(:,:,:)
complex(rk), intent(in) :: boundary(:,:)
complex(rk) :: v(g%nx,g%nz)
integer :: m
do m = 1, g%nx
v(m,:) = matmul(matrix(:,:,m),boundary(m,:))
end do
end function apply_lift
function impose_noslip(g,wall,rhs) result(omega)
type(auxiliary_grid), intent(in) :: g
type(wall_operator), intent(in) :: wall
complex(rk), intent(in) :: rhs(:,:)
complex(rk) :: omega(g%nx,g%nz), psi(g%nx,g%nz), boundary(2), new_boundary(2)
integer :: m
psi = poisson_dirichlet(g,rhs)
omega = rhs
omega(1,:) = 0
do m = 2, g%nx
boundary(1) = sum(g%d(1,:)*psi(m,:))
boundary(2) = sum(g%d(g%nz,:)*psi(m,:))
new_boundary = -matmul(wall%influence_inverse(:,:,m),boundary)
omega(m,:) = rhs(m,:)+matmul(wall%b(:,:,m),new_boundary)
end do
omega(g%nx/2+1,:) = 0
end function impose_noslip
subroutine velocity_from_vorticity(g,omega,mean_u,u,w)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: omega(:,:),mean_u(:)
real(rk), intent(out) :: u(g%nx,g%nz),w(g%nx,g%nz)
complex(rk) :: psi(g%nx,g%nz), c(g%nx,g%nz)
integer :: m
psi = poisson_dirichlet(g,fourier(g,omega))
psi(1,:) = 0
u = physical(g,matmul(psi,transpose(g%d)))+spread(mean_u,1,g%nx)
do m = 1, g%nx
c(m,:) = -cmplx(0.0_rk,g%frequency(m),rk)*psi(m,:)
end do
w = physical(g,c)
end subroutine velocity_from_vorticity
function recover_pressure(g,u,w,fx,fz,diffusion) result(p)
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: u(:,:),w(:,:),fx(:,:),fz(:,:),diffusion
real(rk) :: p(g%nx,g%nz), matrix(g%nz,g%nz), c(g%nz), antiderivative(g%nz+1), values(g%nz)
complex(rk) :: rho(g%nx,g%nz), wall(g%nx,g%nz), phat(g%nx,g%nz), rhs(g%nz)
integer :: m,i,n
if (any(shape(u) /= shape(w))) error stop "velocity shape mismatch"
rho = fourier(g,dx(g,fx)+dz(g,fz))
wall = fourier(g,diffusion*laplacian(g,w)+fz)
phat = 0
do m = 2, g%nx
matrix = g%d2
do i = 1, g%nz
matrix(i,i) = matrix(i,i)-g%frequency(m)**2
end do
matrix(1,:) = g%d(1,:)
matrix(g%nz,:) = g%d(g%nz,:)
rhs = rho(m,:)
rhs(1) = wall(m,1)
rhs(g%nz) = wall(m,g%nz)
phat(m,:) = matmul(inverse_matrix(matrix),rhs)
end do
! The horizontal mean obeys p_z = mean(F_z); integrate its Chebyshev series.
c = matmul(g%iv,sum(fz,dim=1)/g%nx)
antiderivative = 0
antiderivative(2) = c(1)
antiderivative(3) = c(2)/4
do n = 2, g%nz-1
antiderivative(n+2) = antiderivative(n+2)+c(n+1)/(2*(n+1))
antiderivative(n) = antiderivative(n)-c(n+1)/(2*(n-1))
end do
values = 0
do n = 0, g%nz
values = values+antiderivative(n+1)*cos(n*acos(2*g%z-1))/2
end do
phat(1,:) = cmplx(values,0.0_rk,rk)
p = physical(g,phat)
p = p-integral(g,p)
end function recover_pressure
end module wall_stokes
Chapter 18：从 Stokes 到 Navier–Stokes 与 Boussinesq
1. 目标
把体积扩散、壁修正、对流、浮力和温度方程放入同一个二阶预测校正时间步，并固定状态、诊断和检查点的语义。
2. 完整代码
module boussinesq_solver
use cheb_box, only: rk
use slab_grid
use slab_heat
use wall_stokes
use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
implicit none
private
public :: flow_state, flow_solver, configure_solver, zero_state, set_velocity, get_velocity, advance_flow
public :: flow_forcing, diagnostics, flow_pressure, save_checkpoint, load_checkpoint, external_force
type :: flow_state
real(rk), allocatable :: omega(:,:),mean_u(:),theta(:,:)
real(rk) :: time = 0
end type flow_state
type :: flow_solver
type(auxiliary_grid) :: grid
type(wall_operator) :: wall
type(heat_statistics) :: stats
real(rk) :: pr = 1, ra = 0, dt = 0.001_rk, tolerance = 1.0e-10_rk
integer :: k = 12, min_level = 1, steps = 0
logical :: nonlinear = .true., thermal_feedback = .true.
end type flow_solver
abstract interface
subroutine external_force(g,t,fx,fz,ft)
import rk, auxiliary_grid
type(auxiliary_grid), intent(in) :: g
real(rk), intent(in) :: t
real(rk), intent(out) :: fx(g%nx,g%nz),fz(g%nx,g%nz),ft(g%nx,g%nz)
end subroutine external_force
end interface
contains
subroutine configure_solver(s,nx,nz,k,dt,pr,ra,tolerance,min_level)
type(flow_solver), intent(out) :: s
integer, intent(in) :: nx,nz,k,min_level
real(rk), intent(in) :: dt,pr,ra,tolerance
if (.not. all(ieee_is_finite([dt,pr,ra,tolerance]))) error stop "nonfinite solver parameters"
if (dt <= 0 .or. pr <= 0 .or. tolerance <= 0 .or. k < 4 .or. k > 16) error stop "invalid solver parameters"
s%pr = pr
s%ra = ra
s%dt = dt
s%tolerance = tolerance
s%k = k
s%min_level = min_level
call make_grid(s%grid,nx,nz)
call make_wall_operator(s%grid,pr,dt,s%wall)
end subroutine configure_solver
subroutine zero_state(s,state)
type(flow_solver), intent(in) :: s
type(flow_state), intent(out) :: state
allocate(state%omega(s%grid%nx,s%grid%nz),state%theta(s%grid%nx,s%grid%nz),state%mean_u(s%grid%nz))
state%omega = 0
state%theta = 0
state%mean_u = 0
state%time = 0
end subroutine zero_state
subroutine set_velocity(s,state,u,w)
type(flow_solver), intent(in) :: s
type(flow_state), intent(inout) :: state
real(rk), intent(in) :: u(:,:),w(:,:)
state%mean_u = sum(u,dim=1)/s%grid%nx
state%omega = dx(s%grid,w)-dz(s%grid,u)
state%omega = state%omega-spread(sum(state%omega,dim=1)/s%grid%nx,1,s%grid%nx)
end subroutine set_velocity
subroutine get_velocity(s,state,u,w)
type(flow_solver), intent(in) :: s
type(flow_state), intent(in) :: state
real(rk), intent(out) :: u(s%grid%nx,s%grid%nz),w(s%grid%nx,s%grid%nz)
call velocity_from_vorticity(s%grid,state%omega,state%mean_u,u,w)
end subroutine get_velocity
subroutine physical_forcing(s,state,fx,fz,ft,external)
type(flow_solver), intent(in) :: s
type(flow_state), intent(in) :: state
real(rk), intent(out) :: fx(:,:),fz(:,:),ft(:,:)
procedure(external_force), optional :: external
real(rk) :: u(s%grid%nx,s%grid%nz),w(s%grid%nx,s%grid%nz)
real(rk) :: uf(s%grid%nxf,s%grid%nzf),wf(s%grid%nxf,s%grid%nzf)
real(rk) :: ex(size(fx,1),size(fx,2)),ez(size(fx,1),size(fx,2)),et(size(fx,1),size(fx,2))
associate(g => s%grid)
call get_velocity(s,state,u,w)
fx = 0
fz = s%ra*s%pr*state%theta
ft = 0
if (s%thermal_feedback) ft = w
if (s%nonlinear) then
uf = fine_values(g,u)
wf = fine_values(g,w)
fx = -restrict_product(g,uf*fine_values(g,dx(g,u))+wf*fine_values(g,dz(g,u)))
fz = fz-restrict_product(g,uf*fine_values(g,dx(g,w))+wf*fine_values(g,dz(g,w)))
ft = ft-restrict_product(g,uf*fine_values(g,dx(g,state%theta))+wf*fine_values(g,dz(g,state%theta)))
end if
if (present(external)) then
call external(g,state%time,ex,ez,et)
fx = fx+ex
fz = fz+ez
ft = ft+et
end if
end associate
end subroutine physical_forcing
function flow_forcing(s,state,external) result(source)
type(flow_solver), intent(in) :: s
type(flow_state), intent(in) :: state
procedure(external_force), optional :: external
real(rk) :: source(s%grid%nx,s%grid%nz,3),fx(s%grid%nx,s%grid%nz),fz(s%grid%nx,s%grid%nz)
call physical_forcing(s,state,fx,fz,source(:,:,3),external)
source(:,:,1) = dx(s%grid,fz)-dz(s%grid,fx)
source(:,:,1) = source(:,:,1)-spread(sum(source(:,:,1),dim=1)/s%grid%nx,1,s%grid%nx)
source(:,:,2) = spread(sum(fx,dim=1)/s%grid%nx,1,s%grid%nx)
end function flow_forcing
function source_potential(s,state,external) result(r)
type(flow_solver), intent(in) :: s
type(flow_state), intent(in) :: state
procedure(external_force), optional :: external
complex(rk) :: r(s%grid%nx,s%grid%nz,3)
real(rk) :: source(s%grid%nx,s%grid%nz,3)
integer :: c
source = flow_forcing(s,state,external)
do c = 1, 3
r(:,:,c) = poisson_dirichlet(s%grid,fourier(s%grid,source(:,:,c)))
if (c < 3) r(:,:,c) = r(:,:,c)/s%pr
end do
end function source_potential
subroutine advance_flow(s,state,next,external,corrections)
type(flow_solver), intent(inout) :: s
type(flow_state), intent(in) :: state
type(flow_state), intent(out) :: next
procedure(external_force), optional :: external
integer, intent(in), optional :: corrections
complex(rk) :: r0(s%grid%nx,s%grid%nz,3),r1(s%grid%nx,s%grid%nz,3),slope(s%grid%nx,s%grid%nz,3)
complex(rk) :: om(s%grid%nx,s%grid%nz),boundary(s%grid%nx,2),rhs(s%grid%nx,s%grid%nz)
real(rk) :: v(s%grid%nx,s%grid%nz),evolved(s%grid%nx,s%grid%nz),u(s%grid%nx,s%grid%nz),w(s%grid%nx,s%grid%nz),spacing,cfl
integer :: iteration,ncorrect,j,c
ncorrect = 2
if (present(corrections)) ncorrect = corrections
if (ncorrect < 1) error stop "at least one endpoint correction required"
call get_velocity(s,state,u,w)
cfl = 0
do j = 2, s%grid%nz-1
spacing = min(s%grid%z(j)-s%grid%z(j-1),s%grid%z(j+1)-s%grid%z(j))
cfl = max(cfl,s%dt*maxval(s%grid%nx*abs(u(:,j))+abs(w(:,j))/spacing))
end do
if (s%nonlinear .and. cfl > 0.9_rk) error stop "advection CFL too large; reduce dt"
if (s%dt*sqrt(abs(s%ra*s%pr)) > 0.5_rk) error stop "coupling timestep too large"
call zero_state(s,next)
associate(g => s%grid)
om = fourier(g,state%omega)
boundary(:,1) = om(:,1)
boundary(:,2) = om(:,g%nz)
r0 = source_potential(s,state,external)
r1 = r0
do iteration = 0, ncorrect
! Integrate the affine endpoint source exactly: q = (-D Lap_D)^(-1)(r1-r0)/h.
do c = 1, 3
slope(:,:,c) = poisson_dirichlet(g,(r1(:,:,c)-r0(:,:,c))/s%dt)
if (c < 3) slope(:,:,c) = slope(:,:,c)/s%pr
end do
v = physical(g,om-apply_lift(g,g%lift,boundary)-r0(:,:,1)+slope(:,:,1))
evolved = heat_on_grid(g,v,s%pr,s%dt,s%k,s%tolerance,s%min_level,s%stats)
rhs = fourier(g,evolved)+r1(:,:,1)-slope(:,:,1)+apply_lift(g,s%wall%j,boundary/s%dt)
next%omega = physical(g,impose_noslip(g,s%wall,rhs))
v = spread(state%mean_u,1,g%nx)-physical(g,r0(:,:,2)-slope(:,:,2))
evolved = heat_on_grid(g,v,s%pr,s%dt,s%k,s%tolerance,s%min_level,s%stats)+physical(g,r1(:,:,2)-slope(:,:,2))
next%mean_u = sum(evolved,dim=1)/g%nx
v = state%theta-physical(g,r0(:,:,3)-slope(:,:,3))
next%theta = heat_on_grid(g,v,1.0_rk,s%dt,s%k,s%tolerance,s%min_level,s%stats)+physical(g,r1(:,:,3)-slope(:,:,3))
next%time = state%time+s%dt
if (iteration < ncorrect) r1 = source_potential(s,next,external)
end do
end associate
if (.not. all(ieee_is_finite(next%omega)) .or. .not. all(ieee_is_finite(next%theta))) error stop "nonfinite flow state"
s%steps = s%steps+1
end subroutine advance_flow
function diagnostics(s,state) result(d)
type(flow_solver), intent(in) :: s
type(flow_state), intent(in) :: state
real(rk) :: d(14),u(s%grid%nx,s%grid%nz),w(s%grid%nx,s%grid%nz),tz(s%grid%nx,s%grid%nz),v(s%grid%nx,s%grid%nz)
call get_velocity(s,state,u,w)
associate(g => s%grid)
tz = dz(g,state%theta)
d(1) = state%time
d(2) = (inner_product(g,u,u)+inner_product(g,w,w))/2
d(3) = 1+inner_product(g,w,state%theta)
d(4) = 1-sum(tz(:,1))/g%nx
d(5) = 1-sum(tz(:,g%nz))/g%nx
d(6) = integral(g,state%theta)
d(7) = integral(g,state%theta*spread(1-g%z,1,g%nx))
d(8) = integral(g,state%theta*spread(g%z,1,g%nx))
d(9) = max(maxval(abs(u(:,[1,g%nz]))),maxval(abs(w(:,[1,g%nz]))))
d(10) = maxval(abs(state%theta(:,[1,g%nz])))
d(11) = maxval(abs(dx(g,u)+dz(g,w)))
d(12) = inner_product(g,state%theta,state%theta)/2
v = dx(g,u)
d(13) = inner_product(g,v,v)
v = dz(g,u)
d(13) = d(13)+inner_product(g,v,v)
v = dx(g,w)
d(13) = d(13)+inner_product(g,v,v)
v = dz(g,w)
d(13) = s%pr*(d(13)+inner_product(g,v,v))
v = dx(g,state%theta)
d(14) = inner_product(g,v,v)+inner_product(g,tz,tz)
end associate
end function diagnostics
function flow_pressure(s,state,external) result(p)
type(flow_solver), intent(in) :: s
type(flow_state), intent(in) :: state
procedure(external_force), optional :: external
real(rk) :: p(s%grid%nx,s%grid%nz),fx(s%grid%nx,s%grid%nz),fz(s%grid%nx,s%grid%nz),ft(s%grid%nx,s%grid%nz),u(s%grid%nx,s%grid%nz),w(s%grid%nx,s%grid%nz)
call get_velocity(s,state,u,w)
call physical_forcing(s,state,fx,fz,ft,external)
p = recover_pressure(s%grid,u,w,fx,fz,s%pr)
end function flow_pressure
subroutine save_checkpoint(s,state,path)
type(flow_solver), intent(in) :: s
type(flow_state), intent(in) :: state
character(*), intent(in) :: path
integer :: unit
open(newunit=unit,file=path,access="stream",form="unformatted",status="replace")
write(unit) "RBFGT002",s%grid%nx,s%grid%nz,s%k,s%min_level,s%dt,s%pr,s%ra,s%tolerance,s%nonlinear,s%thermal_feedback
write(unit) state%time,state%omega,state%mean_u,state%theta
close(unit)
end subroutine save_checkpoint
subroutine load_checkpoint(s,state,path)
type(flow_solver), intent(in) :: s
type(flow_state), intent(out) :: state
character(*), intent(in) :: path
character(8) :: magic
integer :: unit,nx,nz,k,level,ios
real(rk) :: dt,pr,ra,tol
logical :: nonlinear,thermal
open(newunit=unit,file=path,access="stream",form="unformatted",status="old")
read(unit,iostat=ios) magic,nx,nz,k,level,dt,pr,ra,tol,nonlinear,thermal
if (ios /= 0 .or. magic /= "RBFGT002") error stop "invalid checkpoint header"
if (nx /= s%grid%nx .or. nz /= s%grid%nz .or. k /= s%k .or. level /= s%min_level) error stop "checkpoint grid mismatch"
if (abs(dt-s%dt)+abs(pr-s%pr)+abs(ra-s%ra)+abs(tol-s%tolerance) > 1.0e-15_rk) error stop "checkpoint parameter mismatch"
if ((nonlinear .neqv. s%nonlinear) .or. (thermal .neqv. s%thermal_feedback)) error stop "checkpoint mode mismatch"
call zero_state(s,state)
read(unit,iostat=ios) state%time,state%omega,state%mean_u,state%theta
close(unit)
if (ios /= 0) error stop "truncated checkpoint"
if (.not. all(ieee_is_finite(state%omega)) .or. .not. all(ieee_is_finite(state%theta)) .or. .not. all(ieee_is_finite(state%mean_u))) error stop "nonfinite checkpoint"
end subroutine load_checkpoint
end module boussinesq_solver
Chapter 19：单位周期盒的独立线性参考
1. 目标
根据当前几何计算临界 Rayleigh 数和近临界增长率。单位宽周期盒只允许 \(a_m=2\pi m\)，不能直接采用连续最优波数对应的经典阈值。
2. 完整代码
module rb_linear_reference
use cheb_box, only: rk
use numerical_tools, only: pi
implicit none
private
public :: neutral_rayleigh, reference_growth, reference_mode
contains
subroutine characteristic_roots(ra,a,sigma,qr,qc)
real(rk), intent(in) :: ra,a,sigma
complex(rk), intent(out) :: qr,qc
real(rk) :: p,q,discriminant,u,v
! Pr=1: q^3 - 2*sigma*q^2 + sigma^2*q + Ra*a^2 = 0.
p = -sigma*sigma/3
q = 2*sigma**3/27+ra*a*a
discriminant = (q/2)**2+(p/3)**3
if (discriminant <= 0) error stop "reference outside one-real-root branch"
v = sign(abs(-q/2-sqrt(discriminant))**(1.0_rk/3),-q/2-sqrt(discriminant))
u = -p/(3*v)
qr = cmplx(u+v+2*sigma/3,0.0_rk,rk)
qc = cmplx(-(u+v)/2+2*sigma/3,sqrt(3.0_rk)*(u-v)/2,rk)
end subroutine characteristic_roots
function basis(ra,a,sigma,z,derivative,thermal,odd) result(v)
real(rk), intent(in) :: ra,a,sigma,z
integer, intent(in) :: derivative
logical, intent(in) :: thermal
logical, intent(in), optional :: odd
real(rk) :: v(3)
complex(rk) :: qr,qc,sr,sc,vr,vc
logical :: even_derivative
call characteristic_roots(ra,a,sigma,qr,qc)
sr = sqrt(a*a+qr)
sc = sqrt(a*a+qc)
even_derivative = mod(derivative,2) == 0
if (present(odd)) then
if (odd) even_derivative = .not. even_derivative
end if
if (even_derivative) then
vr = sr**derivative*cosh(sr*(z-0.5_rk))
vc = sc**derivative*cosh(sc*(z-0.5_rk))
else
vr = sr**derivative*sinh(sr*(z-0.5_rk))
vc = sc**derivative*sinh(sc*(z-0.5_rk))
end if
if (present(odd)) then
if (odd) then
if (abs(sr) > 1.0e-12_rk) then
vr = vr/sr
else
vr = 0
if (derivative == 0) vr = z-0.5_rk
if (derivative == 1) vr = 1
end if
end if
end if
if (thermal) then
vr = vr/(sigma-qr)
vc = vc/(sigma-qc)
end if
v = [real(vc,rk),aimag(vc),real(vr,rk)]
end function basis
function cross(a,b) result(c)
real(rk), intent(in) :: a(3),b(3)
real(rk) :: c(3)
c = [a(2)*b(3)-a(3)*b(2),a(3)*b(1)-a(1)*b(3),a(1)*b(2)-a(2)*b(1)]
end function cross
real(rk) function determinant(ra,a,sigma,odd) result(value)
real(rk), intent(in) :: ra,a,sigma
logical, intent(in), optional :: odd
real(rk) :: r1(3),r2(3),r3(3)
r1 = basis(ra,a,sigma,1.0_rk,0,.false.,odd)
r2 = basis(ra,a,sigma,1.0_rk,1,.false.,odd)
r3 = basis(ra,a,sigma,1.0_rk,0,.true.,odd)
r1 = r1/sqrt(sum(r1*r1))
r2 = r2/sqrt(sum(r2*r2))
r3 = r3/sqrt(sum(r3*r3))
value = dot_product(cross(r1,r2),r3)
end function determinant
real(rk) function neutral_rayleigh(mode,odd) result(ra)
integer, intent(in) :: mode
logical, intent(in), optional :: odd
real(rk) :: a,lo,hi,mid,flo,fhi,fmid
integer :: i
if (mode < 1 .or. mode > 8) error stop "reference mode outside scan range"
a = 2*pi*mode
lo = 100
flo = determinant(lo,a,0.0_rk,odd)
do i = 1, 200
hi = lo*1.1_rk
fhi = determinant(hi,a,0.0_rk,odd)
if (flo*fhi < 0) exit
lo = hi
flo = fhi
end do
if (i > 200) error stop "neutral branch not bracketed"
do i = 1, 100
mid = (lo+hi)/2
fmid = determinant(mid,a,0.0_rk,odd)
if (flo*fmid <= 0) then
hi = mid
else
lo = mid
flo = fmid
end if
if (hi-lo < 1.0e-9_rk) exit
end do
ra = (lo+hi)/2
end function neutral_rayleigh
real(rk) function reference_growth(ra,mode) result(sigma)
real(rk), intent(in) :: ra
integer, intent(in) :: mode
real(rk) :: a,lo,hi,mid,flo,fhi,fmid
integer :: i
a = 2*pi*mode
lo = -10
hi = 10
flo = determinant(ra,a,lo)
fhi = determinant(ra,a,hi)
if (flo*fhi >= 0) error stop "growth rate not bracketed on near-onset branch"
do i = 1, 100
mid = (lo+hi)/2
fmid = determinant(ra,a,mid)
if (flo*fmid <= 0) then
hi = mid
else
lo = mid
flo = fmid
end if
if (hi-lo < 1.0e-12_rk) exit
end do
sigma = (lo+hi)/2
end function reference_growth
subroutine reference_mode(ra,mode,sigma,z,w,dw,theta)
real(rk), intent(in) :: ra,sigma,z(:)
integer, intent(in) :: mode
real(rk), intent(out) :: w(size(z)),dw(size(z)),theta(size(z))
real(rk) :: a,c(3),normalization
integer :: i
a = 2*pi*mode
c = cross(basis(ra,a,sigma,1.0_rk,0,.false.),basis(ra,a,sigma,1.0_rk,1,.false.))
normalization = dot_product(c,basis(ra,a,sigma,0.5_rk,0,.false.))
if (abs(normalization) < tiny(1.0_rk)) error stop "singular eigenfunction normalization"
c = c/normalization
do i = 1, size(z)
w(i) = dot_product(c,basis(ra,a,sigma,z(i),0,.false.))
dw(i) = dot_product(c,basis(ra,a,sigma,z(i),1,.false.))
theta(i) = dot_product(c,basis(ra,a,sigma,z(i),0,.true.))
end do
end subroutine reference_mode
end module rb_linear_reference
Chapter 20：完整 RB 运行入口
1. 目标
把求解器组织成一个可运行程序：配置参数、生成初场、推进、输出场和诊断、保存检查点。
2. 完整代码
program run_rb
use cheb_box, only: rk
use numerical_tools, only: pi
use boussinesq_solver
use rb_linear_reference
implicit none
type(flow_solver) :: solver
type(flow_state) :: state,next
real(rk), allocatable :: u(:,:),w(:,:),vertical(:),derivative(:),thermal(:)
real(rk) :: dt,final_time,ra,critical,sigma,amplitude,d(14),old(14),flux_integral(5),begin(14),fit_t,fit_y,fit_tt,fit_ty,fit_n,growth,y
integer :: nx,nz,k,steps,n,j,unit,field_unit
character(256) :: argument,prefix,mode
mode = "super"
dt = 0.002_rk
final_time = 2
nx = 24
nz = 33
k = 12
ra = 0
prefix = "rb"
if (command_argument_count() >= 1) call get_command_argument(1,mode)
if (command_argument_count() >= 2) then
call get_command_argument(2,argument)
read(argument,*) dt
end if
if (command_argument_count() >= 3) then
call get_command_argument(3,argument)
read(argument,*) final_time
end if
if (command_argument_count() >= 4) then
call get_command_argument(4,argument)
read(argument,*) nx
end if
if (command_argument_count() >= 5) then
call get_command_argument(5,argument)
read(argument,*) nz
end if
if (command_argument_count() >= 6) then
call get_command_argument(6,argument)
read(argument,*) ra
end if
if (command_argument_count() >= 7) call get_command_argument(7,prefix)
critical = neutral_rayleigh(1)
if (ra <= 0) ra = 1.2_rk*critical
sigma = reference_growth(ra,1)
steps = nint(final_time/dt)
if (steps < 1 .or. abs(steps*dt-final_time) > 1.0e-10_rk) error stop "final time must be an integer number of steps"
call configure_solver(solver,nx,nz,k,dt,1.0_rk,ra,2.0e-10_rk,1)
solver%nonlinear = trim(mode) /= "linear"
call zero_state(solver,state)
allocate(u(nx,nz),w(nx,nz),vertical(nz),derivative(nz),thermal(nz))
call reference_mode(ra,1,sigma,solver%grid%z,vertical,derivative,thermal)
amplitude = 0.01_rk
if (trim(mode) == "linear" .or. trim(mode) == "onset") amplitude = 1.0e-5_rk
if (trim(mode) == "conduction") amplitude = 0
do j = 1, nz
u(:,j) = -amplitude*derivative(j)*sin(2*pi*solver%grid%x)/(2*pi)
w(:,j) = amplitude*vertical(j)*cos(2*pi*solver%grid%x)
state%theta(:,j) = amplitude*thermal(j)*cos(2*pi*solver%grid%x)
end do
call set_velocity(solver,state,u,w)
open(newunit=unit,file=trim(prefix)//"_history.csv",status="replace")
write(unit,'(a)') "t,E,NuV,Nub,Nut,Q,Qb,Qt,wall,temp_wall,divergence,theta_energy,viscous_dissipation,thermal_dissipation"
d = diagnostics(solver,state)
write(unit,'(*(es24.16,:,","))') d
old = d
begin = d
flux_integral = 0
fit_t=0; fit_y=0; fit_tt=0; fit_ty=0; fit_n=0
call prini(6,0)
do n = 1, steps
call advance_flow(solver,state,next)
state = next
d = diagnostics(solver,state)
write(unit,'(*(es24.16,:,","))') d
flush(unit)
flux_integral(1:3) = flux_integral(1:3)+dt*([old(4)-old(5),old(4)-old(3),old(3)-old(5)]+[d(4)-d(5),d(4)-d(3),d(3)-d(5)])/2
flux_integral(4) = flux_integral(4)+dt*(ra*(old(3)+d(3)-2)-old(13)-d(13))/2
flux_integral(5) = flux_integral(5)+dt*(old(3)+d(3)-2-old(14)-d(14))/2
old = d
if (d(1) >= final_time/4 .and. d(2) > tiny(1.0_rk)) then
y = log(sqrt(d(2)))
fit_n = fit_n+1
fit_t = fit_t+d(1)
fit_tt = fit_tt+d(1)**2
fit_y = fit_y+y
fit_ty = fit_ty+d(1)*y
end if
if (maxval(d(9:11)) > 1.0e-7_rk) error stop "RB wall/divergence budget exceeded"
end do
close(unit)
growth = 0
if (fit_n > 1) growth = (fit_ty-fit_t*fit_y/fit_n)/(fit_tt-fit_t**2/fit_n)
open(newunit=unit,file=trim(prefix)//"_summary.txt",status="replace")
write(unit,'(a,es24.16)') "critical Ra (symmetric m=1 ODE): ",critical
write(unit,'(a,es24.16)') "Ra: ",ra
write(unit,'(a,es24.16)') "reference growth: ",sigma
write(unit,'(a,es24.16)') "measured amplitude growth: ",growth
write(unit,'(a,14es24.16)') "final diagnostics: ",d
write(unit,'(a,3es24.16)') "heat window balance: ",d(6:8)-begin(6:8)-flux_integral(1:3)
write(unit,'(a,2es24.16)') "energy window balance: ",d(2)-begin(2)-flux_integral(4),d(12)-begin(12)-flux_integral(5)
write(unit,'(a,5i12)') "FGT calls/min leaves/max leaves/max level/refinements: ",solver%stats%calls,solver%stats%min_leaves,solver%stats%max_leaves,solver%stats%max_level,solver%stats%refinements
write(unit,'(a,2es24.16)') "source/output monitor maxima: ",solver%stats%source_monitor,solver%stats%output_monitor
close(unit)
call save_checkpoint(solver,state,trim(prefix)//"_checkpoint.bin")
call get_velocity(solver,state,u,w)
open(newunit=field_unit,file=trim(prefix)//"_field.csv",status="replace")
write(field_unit,'(a)') "x,z,u,w,theta,T"
do j = 1, nz
do n = 1, nx
write(field_unit,'(*(es24.16,:,","))') solver%grid%x(n),solver%grid%z(j),u(n,j),w(n,j),state%theta(n,j),1-solver%grid%z(j)+state%theta(n,j)
end do
end do
close(field_unit)
print *, "RB run finished: ",trim(prefix)
print *, "E, NuV, wall, divergence:",d(2),d(3),d(9),d(11)
print *, "amplitude growth / reference:",growth,sigma
end program run_rb
Chapter 21：密输出上的连续方程残差
1. 目标
从三个相邻数值状态构造独立时间差分，把速度、温度和压力评价到更密的空间网格，再按连续 PDE 计算残差。
2. 完整代码
module flow_residuals
use cheb_box, only: rk
use slab_grid
use boussinesq_solver
implicit none
private
public :: dense_residual
contains
subroutine dense_residual(s,before,center,after,residual,external)
type(flow_solver), intent(in) :: s
type(flow_state), intent(in) :: before,center,after
real(rk), intent(out) :: residual(5)
procedure(external_force), optional :: external
type(auxiliary_grid) :: dense
real(rk) :: u(s%grid%nx,s%grid%nz),w(s%grid%nx,s%grid%nz),ub(s%grid%nx,s%grid%nz),wb(s%grid%nx,s%grid%nz),ua(s%grid%nx,s%grid%nz),wa(s%grid%nx,s%grid%nz),p(s%grid%nx,s%grid%nz),lag
real(rk), allocatable :: ud(:,:),wd(:,:),td(:,:),pd(:,:),ut(:,:),wt(:,:),tt(:,:),rx(:,:),rz(:,:),rt(:,:),fx(:,:),fz(:,:),ft(:,:)
call make_grid(dense,2*s%grid%nx,2*s%grid%nz-1)
allocate(ud(dense%nx,dense%nz),wd(dense%nx,dense%nz),td(dense%nx,dense%nz),pd(dense%nx,dense%nz),ut(dense%nx,dense%nz),wt(dense%nx,dense%nz),tt(dense%nx,dense%nz))
allocate(rx(dense%nx,dense%nz),rz(dense%nx,dense%nz),rt(dense%nx,dense%nz),fx(dense%nx,dense%nz),fz(dense%nx,dense%nz),ft(dense%nx,dense%nz))
if (abs((center%time-before%time)-(after%time-center%time)) > 1.0e-12_rk) error stop "residual requires symmetric time samples"
lag = after%time-before%time
if (lag <= 0) error stop "invalid residual time samples"
call get_velocity(s,center,u,w)
call get_velocity(s,before,ub,wb)
call get_velocity(s,after,ua,wa)
p = flow_pressure(s,center,external)
ud = evaluate_grid(s%grid,grid_coefficients(s%grid,u),dense%x,dense%z)
wd = evaluate_grid(s%grid,grid_coefficients(s%grid,w),dense%x,dense%z)
td = evaluate_grid(s%grid,grid_coefficients(s%grid,center%theta),dense%x,dense%z)
pd = evaluate_grid(s%grid,grid_coefficients(s%grid,p),dense%x,dense%z)
ut = evaluate_grid(s%grid,grid_coefficients(s%grid,(ua-ub)/lag),dense%x,dense%z)
wt = evaluate_grid(s%grid,grid_coefficients(s%grid,(wa-wb)/lag),dense%x,dense%z)
tt = evaluate_grid(s%grid,grid_coefficients(s%grid,(after%theta-before%theta)/lag),dense%x,dense%z)
fx = 0
fz = s%ra*s%pr*td
ft = 0
if (s%thermal_feedback) ft = wd
if (present(external)) then
call external(dense,center%time,rx,rz,rt)
fx = fx+rx
fz = fz+rz
ft = ft+rt
end if
rx = ut-s%pr*laplacian(dense,ud)+dx(dense,pd)-fx
rz = wt-s%pr*laplacian(dense,wd)+dz(dense,pd)-fz
rt = tt-laplacian(dense,td)-ft
if (s%nonlinear) then
rx = rx+ud*dx(dense,ud)+wd*dz(dense,ud)
rz = rz+ud*dx(dense,wd)+wd*dz(dense,wd)
rt = rt+ud*dx(dense,td)+wd*dz(dense,td)
end if
residual(1) = max(maxval(abs(rx)),maxval(abs(rz)))
residual(2) = maxval(abs(rt))
residual(3) = maxval(abs(dx(dense,ud)+dz(dense,wd)))
residual(4) = max(maxval(abs(ud(:,[1,dense%nz]))),maxval(abs(wd(:,[1,dense%nz]))))
residual(5) = maxval(abs(td(:,[1,dense%nz])))
end subroutine dense_residual
end module flow_residuals
#应用数学#PDE 数值解