forrestguang的博客

这里写自定义目录标题多重网格算法介绍问题描述——五点差分法求解二维泊松方程五点差分法Gauss迭代算法限制算子介绍提升算子二重网格算法多重网格算法多重网格cuda代码编写串行代码mg.c两重网格cuda并行代码jacobi迭代的cuda编程device_jacobiMakefilecuda_mg.cucuda_mg.hmain.c多重网格算法介绍问题描述——五点差分法求解二维泊松方程{−Δu=f, in Ω,u=g, on ∂Ω,\left\{\begin{al

含参优化控制问题描述{min⁡(y(x,μ),u(x,μ))∈Y×UQ(y(x,μ),u(x,μ))+β(x,μ)j(u(x,μ)), s.t. u(x,μ)∈Uad(μ)={v:ua(μ)≤v≤ub(μ)}⊂U,\left\{\begin{aligned} &\min _{(y(\mathbf{x}, \boldsymbol{\mu}), u(\mathbf{x}, \boldsymbol{\mu})) \in Y \times U} Q(y(\mathbf{x},

问题描述二维泊松方程{−Δu=f, in Ω,u=g, on ∂Ω,\left\{\begin{aligned} -\Delta u &=f, \text { in } \Omega, \\ u &= g ,\text { on }\partial \Omega,\\ \end{aligned}\right.{−Δuu​=f, in Ω,=g, on ∂Ω,​

dolfin安装dolfin是我们发现的一个专门求解PDE约束优化问题的软件，或者说是一个python的库，参考网页添加链接描述，这个软件暂时缺乏用户手册，不过网页上有不少demo，本人就是根据网页上的demo来学习使用dolfin的，具体的原理也请参考网页，本人不做过多解释。这个软件的使用需要结合fenics库，参考fenics官网添加链接描述,fenics有用户手册，自己可以到官网下载。下面重点讲解dolfin以及相关的库的安装方法，本人都是在linux服务器安装的首先还是先建一个虚拟环境，本人建

伴随方法介绍引入算例基本形式min⁡u∈Uad,yJ(y,u)=12∥y−yd∥L2(Ω)2+α2∥u∥L2(Ω)2,Ω=[0,1]d. s.t. {−Δy=u,in  Ω,y=y0,on  ∂Ω.\begin{aligned} &\min_{u \in U_{ad},y} J(y,u) = \frac{1}{2}\|y - y_d\|^{2}_{L^{2}(\Omega)} + \frac{\alpha}{2}\|u\|^2_{L^{2}(\Omega)} ,\qu

基本形式{min⁡y,uJ(y,u), s.t. F(y,u)=0,u∈Uad={ua≤u≤ub}.\left\{\begin{array}{l} \min_{y,u} J(y,u),\\ \text { s.t. }\mathbf{F}(y,u) = 0 ,\\ u \in U_{ad} = \{ u_a\leq u \leq u_b\}. \end{array}\right.⎩⎨⎧​miny,u​J(y,u), s.t. F(y,

min⁡y,uJ(y,u)=∫Ω(y−yΩ)2dxdy s.t. −Δy=u+eΩ,x∈Ω,y=0,x∈∂Ω.−1≤u≤1.\begin{array}{cl} \min _{y,u} & J(y,u) = \int_{\Omega} (y - y_{\Omega})^2 dx dy \\ \text { s.t. } & - \Delta y = u + e_{\Omega},\mathbf{x} \in \Omega, \\{ }&y = 0, \math

1：Error: CUDA driver version is insufficient for CUDA runtime version，这个问题一般是任务在登录节点运行出现的，如果放到计算节点（slurm提交）就会消失。2：__global__开头的函数没有被GPU调用，这种情况可能是slurm文件有问题，比如说没有指定分区gpu（许多人直接拿openmp的slurm用，但是openmp指定分区cpu）。3：所有代码都正确，但是结果有问题，经过考察发现这种问题一般是设置的blockdim=(NX,N

稀疏矩阵主要指的是存在大量零元素的矩阵，常见的稀疏矩阵比如说社交媒体的邻接矩阵，微分方程数值解钟差分法形成的矩阵和有限元方法搭建的刚度矩阵。对于稀疏矩阵，如果对于矩阵每个元素都分配内存存储，将会造成大量的内存浪费，同时做矩阵运算的时候由于读取过程中反复读取零元素将会增大内存访问的时间。因此，为了降低存储空间同时提高访问效率，一种高效的稀疏矩阵存储方法势在必行，下面主要介绍两种稀疏矩阵乘法，我们以一个5×55 \times 55×5的稀疏矩阵存储来举例子。

⇒⇒CPU端定义变量⇒⇒将数据加载到GPU⇒⇒在GPU端执行计算⇒⇒将计算结果拷贝回CPU.CPU端的执行需要一个.c文件来编写代码.GPU端的执行需要一个.cu文件来编写代码.由于C++不能直接调用cuda kernel，还需要编写一个.h文件辅助.即对于一个一般的计算任务，CUDA编程至少应该包含.c文件，.cu文件和.h三个文件。CPU和GPU。

直接使用anaconda，下面红色位置使用命令：pip install mpi4py即可，第一次写代码使用导入MPI库，第一次运行代码可能会报错ImportError: DLL load failed点击链接添加链接描述，安装msmpisetup.exe和MS-MPISDK。下面重点介绍MPI的python使用comm = MPI.COMM_WORLD是全局通信器comm.send(msg, dest=1)发送数据，comm.send(data,dest,tag)这是句法。s = comm

MPI,OpenMP的差别run.slurmMakefile九点差分法的矩阵形式omp_orphaning.comp_end.comp_simd_end.c编译这个omp_simd_end.c的时候，需要先module load gcc，这步是为了使用最新版本的gcc，然后需要在Makefile里面写入相关的向量化命令申明函数也需要利用#pragma omp declare simd //申明向量化函数#-----------------------------------

注意以下安装过程都是默认在docker环境下的服务器进行的，用户需要被服务器管理员添加锦docker组1:百度搜索firedrake官网，进去，然后Ctrl + f，搜索docker，找到下面这个部分2：点击进入上面提到的docker image，复制右边的命令：docker pull firedrakeproject/firedrake，把这段命令复制到服务器终端，然后服务器会开始下载跟firedrake相关的库，大概有两三个G.3：下载完了以后，键入命令：docker run -it – p

MPI的特点：每个进程都是独立的处理器，有自己私有的内存空间，不同进程之间只能通过消息发送接收来传递信息OpenMP的程序有一系列线程（thread)控制，其中一个进程可以开启多个线程，每个线程都有自己私有的变量，并且线程之间有共享内存，可以通过共享内存交互信息。我们以最简单的代码来对OpenMP做一个基本的介绍omp_hello.c运行代码首先需要设置环境变量，这里假设使用4个线程，-lm连接数学库，不过对于这个代码来说不需要使用。下面我们来看OpenMP代码的基本组成。1：加载头文件#incl

给定一个泊松方程如下,{−Δu=f in Ω,u=u0 on ∂Ω.\left\{\begin{aligned}-\Delta u &= f &&\text{ in } \Omega, \\ u &=u_0 &&\text { on } \partial \Omega.\end{aligned}\right.{−Δuu​=f=u0​​​ in Ω, on ∂Ω.​其中Ω=[0,1]2\Omega = [0,1]^2Ω=[0,1]2，对于该泊松方程，往往可以采取差分法求出计算区域离散点的近似解，具体做

#include <iostream>#include <cmath>#include <ctime>using namespace std;#define PI M_PIconst int M = 11;const int N = 11;const double area[2][2] = {{0,1},{0,1}};double UU(double x,double y,int order[2]){ if(order[0] == 0 &a

代码主要包括三个部分第一：定义PDE约束第二：定义目标函数第三：编写优化过程管道流问题初始形状参考链接的geo文件。参考链接添加链接描述不知道为什么原链接代码无法收敛，会报错，因此本人做了一些修改本人把u∇uu\nabla uu∇u去掉了import firedrake as fdimport fireshape as fsimport ROLimport fireshape.zoo as fszfrom fireshape import PdeConstraintfrom

fireshape的使用安装以及基本算例参考链接添加链接描述，这里本人将针对这个介绍做一些中文的翻译以及在代码上做一些修改和解释fireshape的安装参考前一篇博客添加链接描述特殊说明，后面可能会用到fireshape.zoo库，这个库使用简单的source firedrake/activate/bin 无法调用，原因未知，这里的处理方式是，写一个pre.sh文件，文件内容是source firedrake/bin/activatecd fireshape-mastersource activa

bumpkflow.py还是这个区域，现在考虑NS方程，其中γ=1/50,f0=0,f1=0,g0=2.5∗(1−y2),g1=0\gamma=1/50,f_0=0,f_1=0,g_0=2.5*(1-y^2),g_1=0γ=1/50,f0​=0,f1​=0,g0​=2.5∗(1−y2),g1​=0

界面问题−∇⋅(σ∇u)+u=5,Ω-\nabla \cdot(\sigma \nabla u) + u = 5,\Omega−∇⋅(σ∇u)+u=5,Ωu∣∂Ω1=0u|_{\partial \Omega_1} = 0u∣∂Ω1​​=0[σ∇u⋅n]=3,in∂Ω2[\sigma\nabla u \cdot n] = 3,in \partial \Omega_2[σ∇u⋅n]=3,in∂Ω2​∂Ω1\partial \Omega_1∂Ω1​为矩形边界∂Ω2\partial \Omega_2∂Ω2

ADMMimport matplotlib.pyplot as pltimport numpy as npimport timefrom matplotlib import cmdef Cholesky(matrix): w = matrix.shape[0] G = np.zeros((w,w))#实际上只用一半的空间就可以完成矩阵分解 for i in range(w): G[i,i] = (matrix[i,i] - np.dot(G[i,:i]...

问题引入给出两个优化问题：问题1min⁡y,uJ(y,u)=12∫Ω(y−yΩ)2+u2dxdy+∫∂ΩeΓyds. s.t. −Δy+y=u+eΩ,x∈Ω,∂y∂n=0,x∈∂Ω.0≤u≤1.\begin{array}{cl} \min _{y,u} & J(y,u) = \frac{1}{2}\int_{\Omega} (y - y_{\Omega})^2 + u^2 dx dy + \int_{\partial \Omega} e_{\Gamma}y

区域分解算法以下都以泊松方程来举例问题引入给定二维空间的泊松方程：−Δu=f,x∈Ω-\Delta u = f,x \in \Omega−Δu=f,x∈Ωu=g,x∈∂Ωu = g,x \in \partial \Omegau=g,x∈∂Ω我们把区域分成若干部分分别求解，然后相邻区域做一个通信，这样就可以得到一个并行算法。下面我们先把区域分成两部分做一个介绍：考虑区域Ω1,Ω2\Omega_1,\Omega_2Ω1​,Ω2​其中红色的是Ω1\Omega_1Ω1​，蓝色的是Ω2\Omega

一维完全气体欧拉方程的引入在这里我们将利用LF格式，MacCormack格式和Roe格式对一维完全气体欧拉方程进行数值求解。{(ρρuE)t+(ρuρu2+pu(E+p))x=0,p=(γ−1)(E−12ρu2),γ=1.4.\left\{\begin{array}{c}\left(\begin{array}{c}\rho \\\rho u \\E\end{array}\right)_{t}+\left(\begin{array}{c}\rho u \\\rho u^{2}+p \\u(

L型区域上的有限元方法之前的传统数值方法基本上都是针对规则的矩形区域，按道理有限元方法的应用就是为不规则区域量身打造，为此本人尝试使用有限元方法来求解L型区域上的泊松方程，使用的同样是Galerkin原理。混合边界首先我们先尝试求解混合边界上的泊松方程，参考以下图像：其中蓝色代表Dirichlet边界，橙色代表Neumann边界，有限元的基础知识本人已经在上一篇博客介绍有限元因此我就不花篇幅介绍基函数以及试探检验函数空间的选取了，本人使用的仍然是矩形有限元基函数。下面本人重点介绍这样的L型区域在

重点强调：以下内容都是本人的上机作业，仅仅为了方便同僚学习交流，希望不会有抄袭我博客以及代码当成自己原创作业的现象发生问题提出{−Δu=f,x∈Ω,u=u0(x),x∈∂Ω0,∂u∂n=g−βu,x∈∂Ω1.\left \{ \begin{array}{r l} - \Delta u = f,x \in \Omega,\\ u = u_{0}(x),x \in \partial \Omega_{0},\\ \frac{\partial u}{\partial n} = g

差分法求解矩形区域椭圆方程Δu=f,x∈Ω=[−1,1]∗[−1,1]\Delta u = f, x \in \Omega = [-1,1]*[-1,1]Δu=f,x∈Ω=[−1,1]∗[−1,1]u=g,x∈∂Ωu = g, x \in \partial \Omegau=g,x∈∂Ω在x轴和y轴上分别取分划长度为dx,dy，令M=int(2/dx)+1,N=int(2/dx)+1M = int(2/dx) + 1,N = int(2/dx) + 1M=int(2/dx)+1,N=int(2/dx)+

三个线性对流方程{ut+ux=0,x∈Ω,t>0,u(x,0)=u0(x),x∈Ω.\left \{ \begin{array}{r l} u_t + u_x = 0,x \in \Omega,t > 0, \\ u(x,0) = u_0(x),x \in \Omega. \end{array} \right.{ut​+ux​=0,x∈Ω,t>0,u(x,0)=u0​(x),x∈Ω.​∙\bullet∙ Intermittentu0(x)={1,x

G-S迭代def GS(A,b,X,N): X_old = X.copy() n = X.shape[0] for k in range(N): X[0,0] = (b[0,0] - A[0:1,1:n]@X_old[1:n,0])/A[0,0] for i in range(1,n): X[i,0] = (b[i,0] - A[i:i + 1,0:i]@X[0:i] - A[i: i + 1,i + 1:n]@X_old[i

稳定性和后验误差问题提出ut+aux=buxx+cu+f,x∈Ω,t>0,u_t + au_x = bu_{xx} + cu + f,x \in \Omega,t > 0,ut​+aux​=buxx​+cu+f,x∈Ω,t>0,u(x,0)=u0(x),x∈Ω,u(x,0) = u_0(x),x \in \Omega,u(x,0)=u0​(x),x∈Ω,u(α,t)=uα(t),t>0,u(\alpha,t) = u_{\alpha}(t),t > 0,u(α,t)=u

#采用五点差分格式，首先给出精确解和右端项的定义。'''def UU(X, order): temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])**2 + 0.5 if order[0]==0 and order[1]==0: return torch.log(temp) if order[0]==1 and order[1]==0: return temp**(-1) * (20*(X[:,0]+X[:,1

利用Galerkin方法def UU(X, order): temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])**2 + 0.5 if order[0]==0 and order[1]==0: return torch.log(temp) if order[0]==1 and order[1]==0: return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X

使用迭代法求解1：## 向前欧拉格式（注意，这种方法不稳定，不可取）```pythonimport numpy as npimport matplotlib.pyplot as pltimport torchimport timeimport torch.nn as nndef UU(prob,X,t,ox,ot):#精确解，将空间坐标点排成（M*N，2）的形式计算 if prob == 1: if ox == [0,0] and ot == 0:

差分法求解矩形区域椭圆方程Δu=f,x∈Ω=[−1,1]∗[−1,1]\Delta u = f, x \in \Omega = [-1,1]*[-1,1]Δu=f,x∈Ω=[−1,1]∗[−1,1]u=g,x∈∂Ωu = g, x \in \partial \Omegau=g,x∈∂Ω在x轴和y轴上分别取分划长度为dx,dy，令M=int(2/dx)+1,N=int(2/dx)+1M = int(2/dx) + 1,N = int(2/dx) + 1M=int(2/dx)+1,N=int(2/dx)+

需要添加人工的Dirichlet边界考虑的问题仍然是矩形区域的泊松方程。如果整个区域都采用Neumann边界条件，即∂ΩN=∂Ω\partial\Omega_{N}=\partial\Omega∂ΩN​=∂Ω，则这个问题没有定解，因为对于某个精确解uuu而言，针对任意的常数C，u+cu+cu+c都是解为此，我们常常需要添加一个人工提供的Dirichlet边界条件，比如说在这里，我们选取{−1}×[−1,1]\{-1\}\times[-1,1]{−1}×[−1,1]为Dirichlet边界，也就是说矩