Petsc求解非线性方程，SNES对象的介绍-1

最新推荐文章于 2023-12-20 08:00:00 发布

Galerkin码农选手

最新推荐文章于 2023-12-20 08:00:00 发布

阅读量978

点赞数 1

分类专栏：高性能计算

本文链接：https://blog.csdn.net/forrestguang/article/details/113820245

版权

高性能计算专栏收录该内容

47 篇文章 32 订阅 ¥239.90 ¥99.00

订阅专栏

非线性方程的求解

引入非线性方程：
$F(X)=\left[\begin{array}{c}2\alpha(x_0-1)+4x_0(x_0^2+x_1^2-0.25)\\2\alpha(x_1-1)+4x_1(x_0^2+x_1^2-0.25) \end{array}\right]$
我们采用基本牛顿法来求解这个优化问题的零点，基本牛顿法： $x_{k+1}=x_k+d,J_Fd=-F(X)$ ， $J_F$ 是Jacobi矩阵。

new.c

下面我们先使用基本的操作来求解这个问题，首先我们需要定义这样一个向量函数GG，需要自己手动求出Jacobi矩阵定义为HH，然后我们需要定一个训练函数Train来更新X。Petsc类似于C语言和python的结合体，无法返回一个数组作为函数的返回值，所以定义GG的时候，我们需要向量U来储存F在向量X的取值，还需要接受 $\alpha$ 作为参数，参考下面代码：
这里注意向量X和U的赋值方式，首先*ax必须是可读的const，否则后面编译会报错，这里还引入了VecGetArray和VecGetArrayRead，这个Get和Restore是成对出现的。

PetscErrorCode GG(Vec U,Vec X,PetscReal b){
    const PetscReal *ax;
    //PetscReal *ax;
    PetscReal *au;
    //PetscInt i;

    VecGetArray(U,&au);
    VecGetArrayRead(X,&ax);//read the Vec X
    //au[0] = PetscExpReal(b*ax[0])/b - ax[1];
    //au[1] = ax[0]*ax[0] + ax[1]*ax[1] - 1;
    au[0] = 2*b*(ax[0] - 1) + 4*ax[0]*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25);
    au[1] = 2*b*(ax[1] - 1) + 4*ax[1]*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25);
    VecRestoreArrayRead(X,&ax);
    VecRestoreArray(U,&au);

    return 0;
}

我们需要事先设定迭代次数N和跳出循环的条件，具体的原理各位参考本人优化博客的代码，下面直接放代码和makefile：
makeffile

include ${PETSC_DIR}/lib/petsc/conf/variables
include ${PETSC_DIR}/lib/petsc/conf/rules
CFLAGS += -pedantic -std=c99

ecjac: ecjac.o
	-${CLINKER} -o ecjac ecjac.o ${PETSC_LIB}
	${RM} ecjac.o


newton: newton.o
	-${CLINKER} -o newton newton.o ${PETSC_LIB}
	${RM} newton.o
new: new.o
	-${CLINKER} -o new new.o ${PETSC_LIB}
	${RM} new.o

# testing
runnewton_1:
	-@../testit.sh newton "-snes_fd" 1 1
runnewton_2:
	-@../testit.sh newton "-snes_mf" 1 2
runecjac_1:
	-@../testit.sh ecjac "" 1 1

runecjac_2:
	-@../testit.sh ecjac "-snes_mf_operator" 1 2



#--------------------------

runnew_1:
	-@../testit.sh new "" 1 1


test_newton: runnewton_1 runnewton_2

test_ecjac: runecjac_1 runecjac_2

test_new: runnew_1 

test: test_newton test_new test_ecjac

# etc

.PHONY: distclean runexpcircle_1 runexpcircle_2 runecjac_1 runecjac_2 runreaction_1 runreaction_2 runreaction_3 runreaction_4 runreaction_5 runreaction_6 test test_expcircle test_ecjac test_reaction

distclean:
	@rm -f *~ expcircle ecjac reaction *tmp

下面是Petsc代码：
new.c

static char help[] = "solve nonlinear system in newton method.\n"
"the dimension is 2.\n";

#include <petsc.h>


extern PetscErrorCode GG(Vec,Vec,PetscReal);
extern PetscErrorCode HH(Mat,Vec,PetscReal);
extern PetscErrorCode Train(Vec,KSP,PetscInt,PetscReal,PetscReal);

int main(int argc,char **argv){
    PetscErrorCode ierr;
    Vec X;
    PetscInt N = 20;
    PetscReal b = 0.5,X0[2] = {1.0,-1.0},eps = 1.0e-10;
    PetscInt ind[2] = {0,1};
    KSP ksp;

    PetscInitialize(&argc,&argv,NULL,help);
    /***如果没有ierr，会警告
    PetscOptionsBegin(PETSC_COMM_WORLD,"new_","options for new",""); 
    PetscOptionsInt("-N","iteration number of optimal problem","new.c",N,&N,NULL);
    PetscOptionsEnd();
    ***/ 
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"new_","options for new"," test ");CHKERRQ(ierr);
    ierr = PetscOptionsInt("-n","iteration number of optimal problem","new.c",N,&N,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnd();CHKERRQ(ierr);
    
    VecCreate(PETSC_COMM_WORLD,&X);
    VecSetSizes(X,PETSC_DECIDE,2);
    VecSetFromOptions(X);
    
    VecSetValues(X,2,ind,X0,INSERT_VALUES);
    VecAssemblyBegin(X);
    VecAssemblyEnd(X);

    KSPCreate(PETSC_COMM_WORLD,&ksp);
    KSPSetFromOptions(ksp);
    Train(X, ksp, N, b,eps);
    VecView(X,PETSC_VIEWER_STDOUT_WORLD);
    
    return PetscFinalize();
}

PetscErrorCode GG(Vec U,Vec X,PetscReal b){
    const PetscReal *ax;
    //PetscReal *ax;
    PetscReal *au;
    //PetscInt i;

    VecGetArray(U,&au);
    VecGetArrayRead(X,&ax);//read the Vec X
    //au[0] = PetscExpReal(b*ax[0])/b - ax[1];
    //au[1] = ax[0]*ax[0] + ax[1]*ax[1] - 1;
    au[0] = 2*b*(ax[0] - 1) + 4*ax[0]*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25);
    au[1] = 2*b*(ax[1] - 1) + 4*ax[1]*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25);
    VecRestoreArrayRead(X,&ax);
    VecRestoreArray(U,&au);

    return 0;
}
PetscErrorCode HH(Mat A,Vec X,PetscReal b){
    const PetscReal *ax;
    PetscReal v[4];
    PetscInt row[2] = {0,1},col[2] = {0,1};

    VecGetArrayRead(X,&ax);
    //v[0] = PetscExpReal(b*ax[0]);v[1] = -1.0;
    //v[2] = 2*ax[0];v[3] = 2*ax[1];
    v[0] = 2*b + 4*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25) + 8*ax[0]*ax[0];v[1] = 8*ax[0]*ax[1];
    v[2] = 8*ax[0]*ax[1];v[3] = 2*b + 4*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25) + 8*ax[1]*ax[1];
    MatSetValues(A,2,row,2,col,v,INSERT_VALUES);
    VecRestoreArrayRead(X,&ax);

    return 0;
}
PetscErrorCode Train(Vec X,KSP ksp,PetscInt N,PetscReal b,PetscReal eps){
    PetscInt i;
    PetscReal res;
    Vec dir,grad,fun;
    Mat A;

    MatCreate(PETSC_COMM_WORLD,&A);
    MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);
    MatSetFromOptions(A);
    MatSetUp(A);
    VecDuplicate(X,&dir);
    VecDuplicate(X,&grad);
    VecDuplicate(X,&fun);

    for(i = 0; i < N; i++){
	GG(fun,X,b);
	VecNorm(fun,NORM_2,&res);
	GG(grad,X,b);
	HH(A,X,b);
	MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    	MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
	KSPSetOperators(ksp,A,A);
    	KSPSolve(ksp,grad,dir);
	if(i%2 == 0){
	    PetscPrintf(PETSC_COMM_WORLD,
			"the %d iteration of residual is %.4f.\n",i,res);
	}	
	if(res > eps){
	    VecAXPY(X,-1.0,dir);
	    
	}
	else{
	    break;
	}
    }
    PetscPrintf(PETSC_COMM_WORLD,
	"the end iteration is %d, residual is %.4f.\n",i + 1,res);
	
    return 0;
}

运行代码
make new
./new -new_n 10
在这里插入图片描述

个人觉得这个代码很基础，没什么必要解释，用到的函数前面几篇博客都介绍过。下面重点介绍SNES对象：

SNES对象

SNES对象的创造搭建配置和前面介绍的Mat几乎一样，这里本人先把代码放上去再开始解释：
newton.c

static char help[] = "solve nonlinear system in newton method.\n"
"the dimension is 2.\n"
"run with snes_mf snes_fd.\n";

#include <petsc.h>

extern PetscErrorCode FormFunction(SNES, Vec,Vec,void*);
int main(int argc,char **argv){
    
    SNES snes;
    Vec x,r;

    PetscInitialize(&argc,&argv,NULL,help);

    VecCreate(PETSC_COMM_WORLD,&x);
    VecSetSizes(x,PETSC_DECIDE,2);
    VecSetFromOptions(x);
    VecSet(x,1.0);
    VecDuplicate(x,&r);

    SNESCreate(PETSC_COMM_WORLD,&snes);
    SNESSetFunction(snes,r,FormFunction,NULL);
    SNESSetFromOptions(snes);
    SNESSolve(snes,NULL,x);
    VecView(x,PETSC_VIEWER_STDOUT_WORLD);

    SNESDestroy(&snes);VecDestroy(&x);VecDestroy(&r);
    
    return PetscFinalize();
} 
PetscErrorCode FormFunction(SNES snes,Vec x,Vec F,void *ctx){
    const PetscReal *ax;
    PetscReal b = 0.5,*au;

    VecGetArrayRead(x,&ax);
    VecGetArray(F,&au);
    au[0] = 2*b*(ax[0] - 1) + 4*ax[0]*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25);
    au[1] = 2*b*(ax[1] - 1) + 4*ax[1]*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25);
    VecRestoreArrayRead(x,&ax);
    VecRestoreArray(F,&au);

    return 0;
}

相信惟一的难点就是这个FormFunction函数来，FormFunction函数接受第一个Vec参数x为输入，第二个Vec参数F来储存函数值。惭愧，最后这个*ctx的作用本人也还没有看懂。
在这里插入图片描述

运行newton.c
make newton
./newton -snes_fd -snes_monitor
上面命令可以看到迭代过程惨差的变化。在使用牛顿迭代的过程中我们需要使用Jacobi矩阵，但是上面这个代码newton.c我们只搭建定义函数F，为此需要通过有限差分来近似求解Jacobi矩阵，而选项-snes_fd（-snes_mf）会利用 $J_{i,j}=\frac{F(x_i+hu_j)-F(x_i)}{h}$ 指导Petsc寻找步长h来近似求解导数。在求解这个问题中使用的默认求解器是-snes_type newtonls(这是一个牛顿求解器).下面我们介绍一下这个过程SNES都做了什么：
1：通过在SNESSetFunction()设置的回调函数(FromFunction)来计算 $F(x_k)$ 。
2：Jacobi矩阵J通过下面的某一种方法来计算得到：
a:使用者直接提供，然后通过SNESSetJacobian()中的回调函数FormJacobian来计算（ecjac.c就是这个么做的)。
b:通过调用 $N^2$ 次FormFunction计算 $F(x_i+hu_j)$ 来组装搭建J，
c:通过一些特殊的向量 $v$ ，调用 $O (N)$ 次FormFunction计算 $F(x_i+hv_j)$ 来组装搭建J。
d:不组装J，通过Krylov方法迭代求解线性方程组，有限差分法计算 $J y$ ，也就是方程组左端项。
3：使用KSP，PC求解线性方程组。
4：求步长 $\lambda_k$ 。
5：更新迭代点 $x_{k+1}$ 。
6：检查收敛性条件。
我们的代码newton.c使用的是2b。在第一步和第二步的b，每次牛顿迭代需要调用N+1次FormFunction。N如果很大就很麻烦。
在这里插入图片描述 ecjac.c
下面我们利用结构体和FormJacobian来计算矩阵，
利用下面所示结构体定义b。

typedef struct{
    PetscReal b;
}AppCtx;

在我们定义的函数FormJacobian(SNES snes,Vec x,Mat J,Mat P,void *ctx)中，J就是Jacobian矩阵，P是预处理器。

PetscErrorCode FormJacobian(SNES snes,Vec x,Mat J,Mat P,void *ctx){
    AppCtx *user = (AppCtx*)ctx;
    const PetscReal *ax,b = user->b;
    PetscReal v[4];
    PetscInt row[2] = {0,1},col[2] = {0,1};

    VecGetArrayRead(x,&ax);
    //v[0] = PetscExpReal(b*ax[0]);v[1] = -1.0;
    //v[2] = 2*ax[0];v[3] = 2*ax[1];
    v[0] = 2*b + 4*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25) + 8*ax[0]*ax[0];v[1] = 8*ax[0]*ax[1];
    v[2] = 8*ax[0]*ax[1];v[3] = 2*b + 4*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25) + 8*ax[1]*ax[1];
    MatSetValues(P,2,row,2,col,v,INSERT_VALUES);
    MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
    VecRestoreArrayRead(x,&ax);
    if(J != P){
	MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
        MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    }

    return 0;

}

下面直接展示代码和运行结果：

static char help[] = "solve nonlinear system by snes object."
"printf the Jacobi matrix";
#include <petsc.h>
typedef struct{
    PetscReal b;
}AppCtx;
extern PetscErrorCode FormJacobian(SNES,Vec, Mat,Mat,void*);
extern PetscErrorCode FormFunction(SNES, Vec,Vec,void*);
int main(int argc,char **argv){
    SNES snes;
    Vec x,r;
    Mat J;
    AppCtx user;

    PetscInitialize(&argc,&argv,NULL,help);
    user.b = 0.5;

    VecCreate(PETSC_COMM_WORLD,&x);
    VecSetSizes(x,PETSC_DECIDE,2);
    VecSetFromOptions(x);
    VecSet(x,1.0);
    VecDuplicate(x,&r);

    MatCreate(PETSC_COMM_WORLD,&J);
    MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
    MatSetFromOptions(J);
    MatSetUp(J);

    SNESCreate(PETSC_COMM_WORLD,&snes);
    SNESSetFunction(snes,r,FormFunction,&user);
    SNESSetJacobian(snes,J,J,FormJacobian,&user);
    SNESSetFromOptions(snes);

    SNESSolve(snes,NULL,x);
    VecView(x,PETSC_VIEWER_STDOUT_WORLD);

    SNESDestroy(&snes);VecDestroy(&x);VecDestroy(&r);MatDestroy(&J);
    
    return PetscFinalize();
} 
PetscErrorCode FormJacobian(SNES snes,Vec x,Mat J,Mat P,void *ctx){
    AppCtx *user = (AppCtx*)ctx;
    const PetscReal *ax,b = user->b;
    PetscReal v[4];
    PetscInt row[2] = {0,1},col[2] = {0,1};

    VecGetArrayRead(x,&ax);
    //v[0] = PetscExpReal(b*ax[0]);v[1] = -1.0;
    //v[2] = 2*ax[0];v[3] = 2*ax[1];
    v[0] = 2*b + 4*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25) + 8*ax[0]*ax[0];v[1] = 8*ax[0]*ax[1];
    v[2] = 8*ax[0]*ax[1];v[3] = 2*b + 4*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25) + 8*ax[1]*ax[1];
    MatSetValues(P,2,row,2,col,v,INSERT_VALUES);
    MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
    VecRestoreArrayRead(x,&ax);
    if(J != P){
	MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
        MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    }

    return 0;

}
PetscErrorCode FormFunction(SNES snes,Vec x,Vec F,void *ctx){
    AppCtx *user = (AppCtx*)ctx;
    const PetscReal *ax,b = user->b;
    PetscReal *au;

    VecGetArrayRead(x,&ax);
    VecGetArray(F,&au);
    au[0] = 2*b*(ax[0] - 1) + 4*ax[0]*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25);
    au[1] = 2*b*(ax[1] - 1) + 4*ax[1]*(ax[0]*ax[0] + ax[1]*ax[1] - 0.25);
    VecRestoreArrayRead(x,&ax);
    VecRestoreArray(F,&au);

    return 0;
}

在这里插入图片描述
最后那个-snes_test_jacobian是在比较差分法求解的Jacobi矩阵和真实的矩阵的差别。

Galerkin码农选手

关注

1
点赞
踩
1

收藏

觉得还不错? 一键收藏
打赏
0
评论
Petsc求解非线性方程，SNES对象的介绍-1

非线性方程的求解引入非线性方程：F(X)=[2α(x0−1)+4x0(x02+x12−0.25)2α(x1−1)+4x1(x02+x12−0.25)]F(X)=\left[\begin{array}{c}2\alpha(x_0-1)+4x_0(x_0^2+x_1^2-0.25)\\2\alpha(x_1-1)+4x_1(x_0^2+x_1^2-0.25) \end{array}\right]F(X)=[2α(x0−1)+4x0(x02+x12−0.25)2α(x1−1)+4x1(x02+x
复制链接

扫一扫