Category Archives: Numeric

SAS PROC OPTMODEL Examples

July 29, 2020NumericZhouYao

Integer programming

proc optmodel;
	num n=544, t1=26842838.69, t2=-9536251.37;
	var x{1..n} binary;
	num v1{1..n}, v2{1..n};
	read data dt into [_N_] v1='第一维度'n v2='第二维度'n; 	

	min k=sum{i in 1..n} x[i];

	con a: t1-1.5<=sum{i in 1..n} v1[i]*x[i]<=t1+1.5;
	con b: t2-1.5<=sum{i in 1..n} v2[i]*x[i]<=t2+1.5;

	solve;
quit;

Continuous

Example 1

Solution: $[0,\cdots,0]$ .

proc optmodel;
number D init 1000;
var x{1..D} >=-32 <=32;
min z=sum{ i in 1..D} (x[i]**2-10*cos(2*constant('pi')*x[i])+10);
solve with lso/ maxgen=2000;
create data answer from [i] x z;
quit;

solve with lso: lso is based on GA.
DON’T forget parentheses after sum{} as sum{} x+y is equal to (sum{} x)+y which is very different from sum{}(x+y).
In create data answer from [i] x z; [i] is used to index x and z.

Example 2

Shifted function:

proc optmodel;
number D init 1000;
number xt{i in 1..D} init i/5;
var x{i in 1..D} init i/10;
min z=sum{ i in 1..D} ((x[i]-xt[i])**2-10*cos(2*constant('pi')*(x[i]-xt[i]))+10);
solve with lso/ maxgen=2000 primalin;
create data answer from [i] x z xt;
quit;

Notes:

primalin is used to tell that initial value of x is used as an agent.

Example 3

Use result of LSO as starting value of NLP.

proc optmodel;
number D init 1000;
number xt{i in 1..D} init i/5;
var x{i in 1..D} init i/10;
min z=sum{ i in 1..D} ((x[i]-xt[i])**2-10*cos(2*constant('pi')*(x[i]-xt[i]))+10);
solve with lso/ maxgen=2000 primalin;
create data answer from [i] x z xt;
quit;

proc optmodel;
number D init 1000;
number xt{i in 1..D} init i/5;
var x{1..D};
read data answer into [i] x;
min z=sum{ i in 1..D} ((x[i]-xt[i])**2-10*cos(2*constant('pi')*(x[i]-xt[i]))+10);
solve with nlp;
create data anser_nlp from [i] x z;
quit;

Matlab常用函数

July 24, 2020Numeric, ProductivityZhouYao

methods(obj), properties(obj)列出成员函数和成员变量。class(obj)查看变量属于哪个类。

OpenFOAM求解器

April 7, 2020NumericZhouYao

	压缩	不可压缩
稳态	rhoSimpleFoam	simpleFoam overSimpleFoam
瞬态	rhoPimpleFoam overRhoPimpleDyMFoam	pisoFoam pimpleFoam overPimpleDyMFoam

单相流

	压缩	不可压缩
双相	compressibleInterFoam compressibleInterDyMFoam	interFoam overInterDyMFoam
多相	multiphaseEulerFoam	interMixingFoam(3) multiphaseInterFoam

多相流

名字中带DyM的表示动网格，over的表示嵌套网格。

OpenFOAM动网格手册(尚未完成)

April 7, 2020NumericZhouYao

OpenFOAM中自带了两种类型的动网格，一是传统的网格，通过边界移动或变形来控制网格点内部的运动，通过局部重划防止网格质量的降低。这种网格的特点是当运动非常复杂时，非常容易出现计算错误。

二是嵌套（overset）网格，这是将每个部件的网格组合在一起，只需要控制单个部件的网格的运动，这样不涉及单元的形状变化，显然要容易得多。

传统网格按照不同的维度，又可以进行分类。比如按运动是显式还是隐式的进行分类。显式运动时，物体的运动是可以提前确定的；隐式运动，比如球浮在水上，需要通过流体的状态来求解物体运动。

也按边界是否变形来分类。刚体运动，比如一个球在正方形区域内运动，外边界为正方开，内边界为圆，两种边界不发生变形。边界变形的运动，比如手捏装有水的气体，多相流也可以视为边界变形的隐式运动，或自由边界运动。

按拓扑是否改变来分类。拓扑就是网格点间的连接关系，新增网格点，把一个单元划分为多个单元均为拓扑改变。拓扑不改变的情形通常有两种：一是不同区域的网格发生相对移动，但内部的网格点间不相对移动。比如将一个正方形网格旋转一个角度。又比如在一个正方形区域内有一个圆形区域，这个圆形区域自转。二是内部网格点发生相对移动，但移动幅度很小，可以忽略不计。比如正方形内的圆形区域绕直线小幅来回运动，此时网格单元会拉伸，但不需要将单元分割或合并。

在OpenFOAM中完整实现网格运动，一是需要指定动网格类型，二是要指定运动求解器，需要指定motionSolverLibs和motionSolver，三是指定流体力学求解器。通常动网格类型中带MotionSolver的表示支持运动求解器，流体力学求解器名字中带DyM的表示支持动网格，带Overset的表示支持嵌套网格。

传统网格

传统网格共有9种。

嵌套网格

目前只有一种：dynamicOversetFvMesh.

运动求解器

运动求解器用于控制网格运动，主要是显式与隐式。最重要的是lib是fvMotionSolvers。

显式指定运动

指定不同时刻的偏移，中间时刻用插值得到：displacementLaplacian。指定不同时刻的速度，中间时刻用插值得到：velocityComponentLaplacian。

指定偏移需要在0/constant中添加pointDisplacement文件。有三种形式：一是固定值，二是通过列表给不同时刻指定不同值，三是在这个文件中直接插入C++代码来计算每一个时刻的值。显然第三种形式最灵活。

pointDisplacement文件设置

它的格式与设定压强、温度等初场类似。实际上也可以通过类似形式来给不同时刻指定不同的边界值。

固定值：

patch
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

时变值：

patch
{
    type            uniformFixedValue;
    uniformValue    table
    (
        (0   (0 0 0))
        (5   (10 0 0))
    );
}

以上表示5时刻偏移为(10,0,0)，注意这里的偏移是累积偏移，指相对于初始时刻的偏移。

用C++代码计算：

//标量声。
patch
{
    type            codedFixedValue;
    value           uniform 0;

    // Name of generated boundary condition
    redirectType    rampedFixedValue;

    code
    #{
        const scalar t = this->db().time().value();
        operator==(min(10, 0.1*t));
    #};
}
//向量场

标量场code块中指定C++代码，此时表示偏移为 $min(10,0.1t)$ 。

隐式运动

sixDoFRigidBodyMotion。

OpenFOAM源码解读（一）：动网格调用流程

April 5, 2020NumericZhouYao

首先明确，OpenFOAM中涉及动网格与基本网格的代码是分散在两种文件夹中的，动网格是派生的：

基本网格数据结构如场等，在src/OpenFOAM中，比如里面的fields文件夹；有限体积法涉及的网格数据结构在src/finiteVolume/fvMesh中。
动网格基本的motionSolver在src/dynamicMesh/motionSolvers中，涉及拓扑改变的polyTopoChange也在src/dynamicMesh中，派生的motionSolver在src/fvMotionSolver中
动网格大致分为两类，一类是拓扑不改变，适用于微小形变，位于src/dynamicFvMesh，另一类是拓扑改变，位于src/topoChangerFvMesh中。

以rhoCentralDyMSolver为例（文件位于applications/solvers/compressible/rhoCentralFoam中），求解器构造动网格是通过代码：
#include "createDynamicFvMesh.H"
调用是通过
mesh.update()

createDynamicFvMesh.H文件调用dynamicFvMesh类构造动网格，这是一个基类，基类通过case文件的constant/dynamicMeshDict取得应该调用哪个派生类，然后通过run time selection table调用子类的构造函数，所以调用子类构造函数的代码是在基类dynamicFvMesh中，输出类似于：
"Selecting dynamicFvMesh "
的消息的代码也在这个基类中。

有些动网格需要motionSolver，通常在网格名中包含了MotionSolver，比如dynamicMotionSolverTopoFvMesh(位于src/topoChangerFvMesh中)，在这个类中会构造motionSolver:
Foam::dynamicMotionSolverTopoFvMesh::dynamicMotionSolverTopoFvMesh ( const IOobject& io ) : topoChangerFvMesh(io), motionPtr_(motionSolver::New(*this)) {}

motionSolver::New(*this)就是构造运动求解器。与dynamicFvMesh类似，motionSolver也会调用派生类的函数，依据的还是dynamicMeshDict和run time selection table.使用motionSolver的代码也是在dynamicMotionSolverTopoFvMesh，比如motionSolver.updateMesh()。

因此假如我们需要了解每种solver可以使用哪些参数，可以去派生类的构造函数中找，比如src/fvMotionSolver/fvMotionSolvers中找。如果需要编写自己的motionSolver，可以参考这些文件。

dynamicMotionSolverTopoFvMesh

考查一下继承过程：

dynamicMotionSolverTopoFvMesh
topoChangerFvMesh
dynamicFvMesh
fvMesh(位于src/finiteVolume/fvMesh)
polyMesh

文首提到的mesh.update()函数中会使用的两个最重要的类就是polyTopoChanger（位于src/dynamicMesh/polyTopoChanger，类成员在topoChangerFvMesh中声明）和motionSolver。polyTopoChanger会调用polyTopoChange，后者负责refinement等等，且polyTopoChange.C文件达到90KB，值得挖掘。

而motionSolver使用的一个主要函数是updateMesh()，这是一个虚函数，完整的定义在不同的类中，参考下文。

displacementLaplacianFvMotionSolver

以displacementLaplacianFvMotionSolver考查一下它关于motionSolver的继承过程，从高到低如下：

displacementLaplacianFvMotionSolver，读取配置信息主要在这个类的构造函数中。
displacementMotionSolver,fvMotionSolver。在displacementMotionSolver类中会访问0/pointDisplacement文件，具体负责读取的是pointVectorField，定义为typedef GeometricField<vector,pointPatchField,pointMesh> pointVectorField。GeometricField位于src/OpenFOAM/fields中。
points0motionSolver。属于基类，位于src\dynamicMesh\motionSolvers\displacement\points0。
motionSolver。

注意fvMotionSolver虽然名字里面带了MotionSolver，但它并非继承自motionSolver，而是基类，内部包含了一个fvMesh成员变量。

对于这个motionSolver，updateMesh()的完整定义在points0motionSolver中。