@疏影横斜水清浅 OpenFOAM不可压求解器的密度是被约掉的，我这里的空气密度取的是1。当时的写法不太严谨，没检查量纲。

算是解决了吧，其实大可不必费力去猜，先写个Info看看这玩意执行的时间点合不合适，然后直接访问fvMesh对象，拿个引用出来就完事了

codeCorrect
#{
volVectorField &U = mesh().lookupObjectRef<volVectorField>("U");
#}

有一项工作需要在每个时间步对物理场进行修正，之前采用的方法是参考fvOptions/corrections的源码编译自己的库(OpenFOAM-v2206)，但是这种方法比较麻烦。现在我注意到coded source的参考文档有这么一句话：

The coded option is available to all primitive field types, i.e.

scalarCodedSource: scalar equations
...
It provides hooks to:

include sources: codeAddSup
constrain values before the equation is solved: codeSetValue
apply corrections after the equation has been solved: codeCorrect

看上去似乎能在求解控制方程的前后使用codeSetValue或codeCorrect调整场量，从源码看似乎传入了一个场的引用：

template<class Type>
void Foam::fv::CodedSource<Type>::correct
(
    GeometricField<Type, fvPatchField, volMesh>& field
)
{
    DebugInfo
        << "fv::CodedSource<" << pTraits<Type>::typeName
        << ">::correct for source " << name_ << endl;

    updateLibrary(name_);
    redirectOption().correct(field);
}

但是最后并没有找到用法...在codeCorrect代码段里试了试变量名称和field都显示变量未声明，不知哪位大佬有正确的使用方法？还是我理解有误，这个函数并不能实现修改场量的功能？

@李东岳 可以的，李老师！不过做这部分内容的时候我对OF的认知还比较浅，很多细节可能还比较糙..

我写过类似的东西，网格可以根据当前的场值进行变动，用于研究风场作用下的积雪沉积/侵蚀。代码和算例可以参考[https://github.com/fightingxiaoxiao/driftScalarDyFoam]，如果有用的话，也可以引下我的论文[https://doi.org/10.3389/feart.2022.822140]。

大致的思路是这样的：

1. 在dynamicMeshDict中定义速度相关的动网格：

dynamicFvMesh dynamicMotionSolverFvMesh;

motionSolverLibs (fvMotionSolvers);

motionSolver velocityLaplacian;   //displacementLaplacian;
diffusivity  quadratic inverseDistance 3(bottom.snow roof top);

2. 定义pointMotionU边界条件为codedFixedValue

    bottom.snow
    {
        type            codedFixedValue;
        value          uniform (0 0 0);

        name         erosionDeposition;
    }

3. 最后在codeDict嵌入代码，这里比较大的一个坑就是动网格使用节点值进行位移，但是实际计算拿到的都是面心值和体心值，需要做一个插值。另外，deltaH是我自己求解器里定义的场，需要替换成你自己的内容。

erosionDeposition
{
    code
    #{
        
        label patchIndex = patch().index();
        const volScalarField& deltaH = this->db().objectRegistry::lookupObject<volScalarField>("deltaH");

        primitivePatchInterpolation facePointInterp(deltaH.mesh().boundaryMesh()[patchIndex]);   //初始化插值类

        const scalarField& deltaHp = deltaH.boundaryField()[patchIndex];
        auto deltaHpp = facePointInterp.faceToPointInterpolate(deltaHp);                      //面心向节点插值
        Info << "Successfully interpolate "<< deltaHpp().size() << " points." << endl;
        Info << "Local points number = " << patch().localPoints().size() << endl;
        pointField pVec(deltaHpp().size());

        const scalar lowerLimit = 0.;

        forAll(pVec, i)
        {
            pVec[i][0] = 0;
            pVec[i][1] = 0;
            pVec[i][2] = deltaHpp()[i];

            if(patch().localPoints()[i][2] + pVec[i][2] * this->db().time().deltaTValue() < lowerLimit)
            {    
                pVec[i][2] = (lowerLimit - patch().localPoints()[i][2])/this->db().time().deltaTValue();
            }
        }

        (*this) == pVec;
    #};

    codeInclude
    #{
        #include "fvCFD.H"
        #include "primitivePatchInterpolation.H"
    #};

    codeOptions 
    #{
        -I$(LIB_SRC)/finiteVolume/lnInclude \
        -I$(LIB_SRC)/meshTools/lnInclude 
    #};

4. 最终效果如下，一个方块周边的时变积雪
   

@星星星星晴 我基本参照的是这个教程https://github.com/BrushXue/OpenFOAM-AppleM1

@李东岳 是的，现在消费级平台已经开始普及DDR5内存了。苹果这个芯片里的CPU和GPU是共享内存的，所以理论内存带宽很大，但是似乎8个CPU核心无法饱和那么大的带宽。

入手了一台MacBook，在macOS上原生编译的OpenFOAM-v2206居然和Ubuntu ARM虚拟机上的OpenFOAM-10半斤八两。

CPU型号：Apple M1 Max 8P+2E (MacBook Pro 16 2021)
系统：macOS Monterey 12.3 真系统
内存：64GB LPDDR5 6400MHz
OpenFOAM版本：OpenFOAM-v2206

8 166.85
4 277.05
2 467.17
1 798.8

CPU型号：Apple M1 Max 8P+2E (MacBook Pro 16 2021)
系统：Ubuntu ARM 20.04 虚拟机
内存：32GB LPDDR5 6400MHz
OpenFOAM版本：OpenFOAM-10

8 154
4 289.836
2 453.656
1 803.381

如果颗粒刚度比较大的话，可以用重叠网格（overset mesh），也就是背景网格+颗粒周边网格，两套网格重合的地方会进行插值处理。这样可以回避网格重画的难点，不过这东西算二维还行，三维想算快算好比较折磨。

Cyclic似乎只适合几何拓扑一致的网格边界吧，你这种情况用cyclicAMI边界可能会好一些。