Eigen库的共轭梯度法(Conjugate Gradient)代码分析

ConjugateGradient类

Eigen使用ConjugateGradient类进行共轭梯度法计算。一个代码示例为(https://eigen.tuxfamily.org/dox/classEigen_1_1ConjugateGradient.html ):

int n = 10000;
VectorXd x(n), b(n);
SparseMatrix<double> A(n,n);
// fill A and b
ConjugateGradient<SparseMatrix<double>, Lower|Upper> cg;
cg.compu[......]

继续阅读

数值密度SPH(Number Density SPH)的一些推导

本文提到的方法见于:Tartakovsky, Alexandre M., and Paul Meakin. “A smoothed particle hydrodynamics model for miscible flow in three-dimensional fractures and the two-dimensional Rayleigh–Taylor instability.” Journal of Computational Physics 207.2 (2005): 610-624.

普通SPH形式

密度: $\rho_i=\sum_{j}m_jW_{ij}.$

其中$[……]

继续阅读

关于SPH通用核函数的一些推导

值的表示

对于一维情况,现假设我们有一个定义在非负实数域上的标量权值函数$\phi_1(x)$,截断至1,即满足:$\int{\phi_1(x)}\mathrm{d}x=1$,且$\phi_1(x)=0$对$\forall x>1$成立。

我们希望基于它设计一个关于$\vec{r}$的对称核函数$W_1(\vec{r})$,截断至$h$,那么这个核函数的形式应为

$$W_1(\vec{r})=\alpha_1 \phi_1\left(\frac{|\vec{r}|}{h}\right).$$

它的积分应为1,因此

$$\int W_1(x)\mathrm{d}x=h\alpha_1\i[……]

继续阅读

Robert Saye:VIIM论文算例简介

3.10收敛性测试和数值验证

首先是每算例都做三组:

  • 固定[latex]\epsilon[/latex],改变网格尺寸h,检查[latex]h\to 0[/latex]的收敛情况。
  • 固定比例关系[latex]\epsilon=\alpha h[/latex],其中[latex]\alpha[/latex]为一常数,检查[latex]h\to 0[/latex]的收敛情况。
  • 交换极限,先算一个[latex]\epsilon\to 0^+[/latex]的内极限,然后检查[latex]h\to 0[/latex]的收敛情况。

所有测试都是在规则正方形/立方体网格上完成的。

3.[……]

继续阅读

PARIS模拟器算例介绍

文章链接:https://doi.org/10.1016/j.cpc.2021.107849

算例部分在第四章。

4.1基础测试

4.1.1泊肃叶流动

设置一个8*8*2网格,模拟二维正方形水槽内的泊肃叶流动。参数[latex]||\nabla p||=\mu=\rho=1[/latex].在水槽左端右端设置压强边界条件。模拟至流场以1e-3精度稳定,约在1000步后,时刻1达成。使用显式viscous diffusion.上下边界处使用粘性边界条件u=0.

4.1.2二维斯托克斯流

和泊肃叶流动相同,但计算区域中心有一直径为0.5的圆盘。不进行advect[……]

继续阅读

关于MAC网格不可压流体求解中投影矩阵的笔记

现假设我们有一流体格i,和与之相邻的格子j。我们研究投影矩阵A中,代表(i,j)格子对的元素(下简记为A[i,j])和i格子所对应的对角元(下简记为A[i,i])的数值情况。注意,i和j都是格子,本身坐标是一个向量,不要与矩阵坐标相混淆。

现列举j的几种情况。

情况一:j亦是流体格。此时,j格压强也是未知数,A[i,j]=-1,而A[i,i]中有一项1,如下图。(或在porous flow中,面积分数为a,那就是-a和a)

情况二:j亦是流体格,但i和j中间有一Neumann边界条件(虽然这很奇怪,我怀疑它会不会出现)。此时格子i的方程中,i-j相关联的项直接被取消,[……]

继续阅读

Robert Bridson书中关于自由表面求解的笔记

Fluid Simulation for Computer Graphics第二版,8.2节。

在求解自由表面流体时,简单地把计算域划分为液体格和气体格,会导致artifact.如果能更精确地追踪表面,可以修改一下在水-气界面附近计算压强差的方式。

根据Gibou et al. 2002的ghost fluid方法,可以把更新流体速度节点[latex]u_{i+1/2,j,k}[/latex]的方程写作

假设(i,j,k)在流体内部,而(i+1,j,k)在空气中,即[latex]\phi_{i,j,k}\leq 0<\phi_{i+1,j,k}[/latex].认为[……]

继续阅读

Ubuntu 20.04 Server软路由配置折腾记录

记录一下最近的一次折腾,使用一台安装了Ubuntu 20.04 Server的工控机搭建软路由,并成功配置家庭NAS的内网访问。计算机网络没好好学,全靠本能debug,感谢lcy同志的全程技术科普(

本文主要参考了https://blog.lcy.im/2017/09/14/%E4%BB%8Eubuntu-server-%E5%BC%80%E5%A7%8B%E9%85%8D%E7%BD%AE%E8%87%AA%E5%8A%A8%E5%88%86%E6%B5%81%E7%9A%84%E8%BD%AF%E8%B7%AF%E7%94%B1/

网络拓扑设计

瞎画的,是这意思就行([……]

继续阅读

Real-time exploration of regular volume data by adaptive reconstruction of isosurfaces读书笔记

本文讨论使用Octree(八叉树)重建表面的问题。

3 Continuous isosurfaces

讨论连续性问题。

如果在Octree上储存数据,那么在不同大小的格子交界的地方,数据可能会不连续。

本文采用的方法:在大小格子交界处,让大格子在小格子上sample数值:

也就是用方块节点插值出来三角节点的值,把原来的扔掉。这样一来,小格子上的数值就是线性的,那么marching cubes就会在大小格子上得到同样的intersection point,也就是说数值是连续的。

但是,虽然数值连续,但算出来的表面不一定连续。如图3:

本文的方法,[……]

继续阅读

关于马拉戈尼效应(Marangoni Effect)

本文基于维基百科内容写作:

https://en.wikipedia.org/wiki/Marangoni_effect

https://en.wikipedia.org/wiki/Marangoni_number

示例视频(需科学上网):

https://upload.wikimedia.org/wikipedia/commons/3/38/Marangoni_effect_experimental_demonstration.ogv

马拉戈尼效应由表面张力的梯度产生。在简化情况下,流体被马拉戈尼效应驱动而流动的特征速度[latex]u\sim\Delta\ga[……]

继续阅读