I am trying to solve the linear finite element equation $M\ddot{u}+Ku=F(t)$, where $M$ is the mass matrix ,$K$ the stiffness matrix and $F(t)$ the external load vector, parallel in time using XBraid package.Using Newmark's method for time integration we take three equations:
- $(K+\frac{M}{dt^2β})u_{n+1}=F_{n+1}+\frac{M}{dt^2β}(u_n+dt\dot{u}_n+(0.5-β)dt^2\ddot{u}_n)$
- $\ddot{u}_{n+1}=\frac{1}{βdt^2}(u_{n+1}-u_n-dt\dot{u}_n-dt^2(0.5-β)\ddot{u}_n)$
- $\dot{u}_{n+1}=\dot{u}_{n}+dt((1-γ)\ddot{u}_{n}+γ\ddot{u}_{n+1})$
We arrange these three unknowns in a block vector $w_{n+1}=[u_{n+1},\ddot{u}_{n+1},\dot{u}_{n+1}]$ and with some further effort we take something like $$ Aw=f $$ where $w=[w_0\;w_1 \dots w_{n+1}]^T$, $f=[f_0\;f_1 \dots f_{n+1}]^T$ and
$$A=\begin{bmatrix} I \\ Φ && I \\ & &Φ&& \ddots \\ &&&&\ddots \\ &&&&&Φ &&&I \\ \end{bmatrix} $$
with $Φ$ be the time step matrix which take the solution from one time to another.
I tried to implement this idea in Xbraid,however I noticed some strange results. For example I noticed that for very small time-points(50 time points) the residual provided by Xbraid grows at each iteration and at some point decreases rapidly and the method converges giving me the same results as the serial problem. For larger time-steps (over 1000 for example) the residual grows rapidly (takes very large values) and after some point it start to decrease very slowly reaching a point where it stucks and no furhter degradation takes place, leading to non-convergence. Of course this behaviour is not normal in both cases, but I don't know what can cause this behaviour. I provide the basic core of Xbraid implementation in case someone familiar with the software can suggest a solution.
/*Time-step Function. Newmark Method is used to take the solution from t->t+1*/
void Newmark_BetaMethod(double t,double dt,double gamma,double beta,int NN,int NEy,int *dbc,double
**K_global,double **M_global,double **Matrix_global,double *Fext,double *RHS,Solver *solver,double
*tempDispl,double *tempVeloc,double *tempAccel,double *displ,double *veloc,double *accel){
/*create Matrix_global*/
for(int i=0;i<2*NN;i++){
for(int j=0;j<2*NN;j++){
Matrix_global[i][j]=K_global[i][j]+(M_global[i][j])/(beta*dt*dt);
}
}
/*hold values at time step t*/
for(int i=0;i<2*NN;i++){
tempDispl[i]=displ[i];
tempVeloc[i]=veloc[i];
tempAccel[i]=accel[i];
}
/*calculate RHS*/
for(int i=0;i<2*NN;i++){
RHS[i]=0.0;
}
for(int i=0;i<2*NN;i++){
for(int j=0;j<2*NN;j++){
RHS[i]+=(M_global[i][j]*(displ[j]+dt*veloc[j]+(0.5-beta)*(dt*dt)*accel[j]))/(beta*dt*dt);
}
RHS[i]+=t*Fext[i];
}
/*zero the rows where homogeneous dirichlet_bc refer to*/
BoundaryCondForRHS(NN,NEy,dbc,RHS);
/*Find displacement at t+1*/
solver->ConjugateGradient(Matrix_global,RHS,displ);
/*Update acceleration*/
for(int i=0;i<2*NN;i++){
accel[i]= (displ[i]-tempDispl[i]-dt*tempVeloc[i]-(dt*dt)*(0.5-beta)*tempAccel[i])/(beta*dt*dt);
}
/*Update velocity*/
for(int i=0;i<2*NN;i++){
veloc[i]= tempVeloc[i]+dt*((1-gamma)*tempAccel[i]+gamma*accel[i]);
}
}
//include all time-dependent variables->create a block vector
class BraidVector{
public:
double *displ; //a displacements vector
double *veloc; //velocities vector
double *accel; //acceleration
BraidVector(double *u,double *v,double *acc);
virtual ~BraidVector() {};
};
BraidVector::BraidVector(double *u,double *v,double *acc){
displ=u;
veloc=v;
accel=acc;
}
//include all time-independent information here
class MyBraidApp:public BraidApp
{
protected:
public:
MyBraidApp(MPI_Comm comm_t_, int rank_,double **Matrix_global_,double **M_global_,double **K_global_,double *Fext_,double *RHS_,double gamma_,double beta_,double *tempD,double *tempV,double *tempA,int NN,int NEy_,int *dbc,Solver *solver_ ,double tstart_ =0.0,double tstop_=1.0,int ntime_=100);
int rank;
double **Matrix_global; //This matrix pointer holds the sum of K+M/beta*dt*dt
double **M_global; //global mass matrix
double **K_global; //global stiffness matrix
double *Fext; //global total load vector
double *RHS; //temporary vector that holds t*Fext+(M/beta*dt*dt)*{displ+dt*velocity+(0.5-beta)*dt*dt*acceleration}
double beta; //Newmark parameter
double gamma; //Newmark parameter
double *tempDispl; //temporary vector holds old values of displacements
double *tempVeloc;
double *tempAccel;
int nodes; //total number of nodes in spatial mesh
int NEy; //number of elemenets in y->helps to find the size of boundary nodes
int *displ_bc; //holds the dirichlet boundary nodes where u=0;
Solver *solver; //pointer to a solver object (needed for solving tthe linear system at each step)
virtual ~MyBraidApp();
virtual int Init(double t,braid_Vector *w_ptr);
virtual int Step(braid_Vector w_,braid_Vector wstop, braid_Vector fstop,BraidStepStatus &pstatus);
virtual int Clone(braid_Vector w_,braid_Vector *s_ptr);
virtual int Sum(double alpha,braid_Vector x,double beta,braid_Vector y);
virtual int BufSize(int *size_ptr,BraidBufferStatus &status);
virtual int BufPack(braid_Vector w_,void *buffer,BraidBufferStatus &status);
virtual int BufUnpack(void *buffer,braid_Vector *w_ptr,BraidBufferStatus &status);
virtual int SpatialNorm(braid_Vector w_,double *norm_ptr);
virtual int Access(braid_Vector w_,BraidAccessStatus &astatus);
virtual int Free(braid_Vector w_);
// Not needed in this example
virtual int Residual(braid_Vector w_,
braid_Vector r_,
BraidStepStatus &pstatus) { return 0; }
// Not needed in this example
virtual int Coarsen(braid_Vector fw_,
braid_Vector *cw_ptr,
BraidCoarsenRefStatus &status) { return 0; }
// Not needed in this example
virtual int Refine(braid_Vector cw_,
braid_Vector *fw_ptr,
BraidCoarsenRefStatus &status) { return 0; }
};
/*Construct the braid app*/
MyBraidApp::MyBraidApp(MPI_Comm comm_t_, int rank_,double **Matrix_global_,double **M_global_,double **K_global_,double *Fext_,double *RHS_,double gamma_,double beta_,double *tempD,double *tempV,double *tempA,int NN,int NEy_,int *dbc,Solver *solver_ ,double tstart_ ,double tstop_,int ntime_)
:BraidApp(comm_t_, tstart_, tstop_, ntime_)
{
rank=rank_;
Matrix_global=Matrix_global_;
M_global=M_global_;
K_global=K_global_;
Fext=Fext_;
RHS=RHS_;
beta=beta_;
gamma=gamma_;
nodes=NN;
tempDispl=tempD;
tempAccel=tempA;
tempVeloc=tempV;
NEy=NEy_;
displ_bc=dbc;
solver=solver_;
}
MyBraidApp::~MyBraidApp(){
delete []tempDispl;
delete []tempVeloc;
delete []tempAccel;
delete []Fext;
delete []RHS;
delete []displ_bc;
int size=2*nodes;
for(int i=0;i<size;i++){
delete M_global[i];
delete Matrix_global[i];
delete K_global[i];
}
delete [] M_global;
delete [] K_global;
delete [] Matrix_global;
}
/*Initialize the block vector w*/
int MyBraidApp::Init(double t,braid_Vector *w_ptr){
int size=2*nodes;
double *displ=new double[size];
double *veloc=new double[size];
double *accel=new double[size];
BraidVector *w=new BraidVector(displ,veloc,accel);
/*w->displ=new double[size];
w->veloc=new double[size];
w->accel=new double[size];
w->size=size;*/
//Initial conditions at t=0
if(t==tstart){
for(int i=0;i<size;i++){
w->displ[i]=0.0;
w->veloc[i]=0.0;
w->accel[i]=0.0;
}
}
else{
//for t>0 use random values as initial guess
for(int i=0;i<size;i++){
w->displ[i]=1e-2;
w->veloc[i]=1e-2;
w->accel[i]=1e-2;
}
}
*w_ptr=(braid_Vector) w;
return 0;
}
/*Takes a time step using Newmark method*/
int MyBraidApp::Step(braid_Vector w_,braid_Vector wstop, braid_Vector fstop,BraidStepStatus &pstatus){
BraidVector *w=(BraidVector*)w_;
double t1,t2;
// Get time step information->Return tstart and tstop
pstatus.GetTstartTstop(&t1, &t2);
double dt=t2-t1; //calculate dt
/*Call Newmark method to advance in time*/
Newmark_BetaMethod(t2,dt,gamma,beta,nodes,NEy,displ_bc,K_global,M_global,Matrix_global,Fext,RHS,solver,tempDispl,tempVeloc,tempAccel,w->displ,w->veloc,w->accel);
return 0;
}
/*Create a new clone vector s of block vector w*/
int MyBraidApp::Clone(braid_Vector w_,braid_Vector *s_ptr){
int size=2*nodes;
BraidVector *w = (BraidVector*) w_;
double *displ=new double[size];
double *veloc=new double[size];
double *accel=new double[size];
BraidVector *s=new BraidVector(displ,veloc,accel);
/*s->displ=new double[size];
s->veloc=new double[size];
s->accel=new double[size];
s->size=size;*/
for(int i=0;i<size;i++){
s->displ[i]=w->displ[i];
s->veloc[i]=w->veloc[i];
s->accel[i]=w->accel[i];
}
*s_ptr=(braid_Vector) s;
return 0;
}
int MyBraidApp::Sum(double alpha,braid_Vector x,double beta,braid_Vector y){
int size=2*nodes;
BraidVector *vec1 = (BraidVector*) x;
BraidVector *vec2 = (BraidVector*) y;
for(int i=0;i<size;i++){
vec2->displ[i]=alpha*(vec1->displ[i])+beta*(vec2->displ[i]);
vec2->veloc[i]=alpha*(vec1->veloc[i])+beta*(vec2->veloc[i]);
vec2->accel[i]=alpha*(vec1->accel[i])+beta*(vec2->accel[i]);
}
return 0;
}
/*Create the buffer's size that holds the data of braid vector*/
int MyBraidApp::BufSize(int *size_ptr,BraidBufferStatus &status){
int buffsize=3*(2*nodes); //buffsize holds the size of the block vector w =size of displ+ size of veloc +size of accel
*size_ptr=buffsize*sizeof(double);
return 0;
}
/*pack displacement,velocity and acceleration vector into a buffer->This buffer is send via proccesses*/
int MyBraidApp::BufPack(braid_Vector w_,void *buffer,BraidBufferStatus &status){
BraidVector *w=(BraidVector*)w_;
double *dbuffer=(double *)buffer;
int size=2*nodes;
for(int i=0;i<size;i++){
dbuffer[i]=w->displ[i];
dbuffer[size+i]=w->veloc[i];
dbuffer[2*size+i]=w->accel[i];
}
status.SetSize((3*size)*sizeof(double));
return 0;
}
/*Create a new braid vector,unpack the data, and initialize it*/
int MyBraidApp::BufUnpack(void *buffer,braid_Vector *w_ptr,BraidBufferStatus &status){
int size=2*nodes;
double *dbuffer = (double *) buffer;
double *displ=new double[size];
double *veloc=new double[size];
double *accel=new double[size];
BraidVector *w=new BraidVector(displ,veloc,accel);
/*w->displ=new double[size];
w->veloc=new double[size];
w->accel=new double[size];*/
for(int i=0;i<size;i++){
w->displ[i]=dbuffer[i];
w->veloc[i]=dbuffer[size+i];
w->accel[i]=dbuffer[2*size+i];
}
*w_ptr = (braid_Vector) w;
return 0;
}
/*Euclidean norm is used for convergece*/
int MyBraidApp::SpatialNorm(braid_Vector w_,double *norm_ptr){
int size=2*nodes;
BraidVector *w=(BraidVector*)w_;
double norm=0.0;
for(int i=0;i<size;i++){
norm+=(w->displ)[i]*(w->displ)[i]+(w->veloc)[i]*(w->veloc)[i]+(w->accel)[i]*(w->accel)[i];
}
*norm_ptr=sqrt(norm);
return 0;
}
/*Clean up*/
int MyBraidApp::Free(braid_Vector w_){
BraidVector *w=(BraidVector*)w_;
delete []w->displ;
delete []w->veloc;
delete []w->accel;
delete w;
return 0;
}
int MyBraidApp::Access(braid_Vector w_,BraidAccessStatus &astatus){
char filename[255];
FILE *file;
BraidVector *w = (BraidVector*) w_;
int size=2*nodes;
int done,iter,index,level;
double t;
astatus.GetTILD(&t, &iter, &level, &done);
astatus.GetTIndex(&index);
// Print information to file
sprintf(filename, "%s.%07d.%05d", "Plane stress", index, rank);
file = fopen(filename, "w");
//fprintf(file, "%d\n", Npoints);
fprintf(file, "%.14e\n", tstart );
fprintf(file, "%.14e\n", tstop );
//fprintf(file, "%d\n", size );
fprintf(file, "%s\n" , "*Displacements*");
for (int i=0;i<size;i++){
fprintf(file, "%.14e\n", w->displ[i]);
}
fprintf(file, "%s\n" , "*Velocities*");
for (int i=0;i<size;i++){
fprintf(file, "%.14e\n", w->veloc[i]);
}
fprintf(file, "%s\n" , "*Accelerations*");
for (int i=0;i<size;i++){
fprintf(file, "%.14e\n", w->accel[i]);
}
fflush(file);
fclose(file);
return 0;
}
I need to mention that I ran Xbraid sequentially as mentioned in the manual and the results were exactly the same as my serial code.
