Solving a coupled differential equations system using time splitting

91 Views Asked by At
/******************************************************************************

                            Online C Compiler.
                Code, Compile, Run and Debug C program online.
Write your code in this editor and press "Run" button to compile and execute it.

*******************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

#define PI  3.141592

void read_input(double *D, double *L, int *nx, double *t_F);

double main(void) {
  /******************************/
  /* Declarations of parameters */
  /******************************/
  /* Number of grid points */                  
  int nx;
  /* Length of domain */
  double L;
  /* Equation coefficients */
  double D;
  /* Length of time to run simulation. */
  double t_F;

  /* Read in from file; */
  read_input(&D, &L, &nx, &t_F);

  /* Grid spacing */
  double dx = L/nx;
  double invdx2 = 1.0/(dx*dx);      
  /* Time step */
  double dt = 0.25/invdx2; // changed to 0.25/dx^2 to satisfy the stability condition

  /************************************************/
  /* Solution Storage at Current / Next time step */
  /************************************************/
  double *uc, *un, *vc, *vn;
  /* Time splitting solutions */
  double *uts1, *uts2, *vts1, *vts2;
  /* Derivative used in finite difference */
  double deriv;

  /* Allocate memory according to size of nx */
  uc = malloc(nx  * sizeof(double));
  un = malloc(nx * sizeof(double));
  vc = malloc(nx * sizeof(double));
  vn = malloc(nx * sizeof(double));
  uts1 = malloc(nx * sizeof(double));
  uts2 = malloc(nx * sizeof(double));
  vts1 = malloc(nx * sizeof(double));
  vts2 = malloc(nx * sizeof(double));

  /* Check the allocation pointers */
  if (uc==NULL||un==NULL||vc==NULL||vn==NULL||uts1==NULL||
  uts2==NULL||vts1==NULL||vts2==NULL) {
    printf("Memory allocation failed\n");
    return 1;
  }
  
  int k;
  double x;
  /* Current time */
  double ctime;

  /* Initialise arrays */
  for(k = 0; k < nx; k++) {
    x = k*dx;
    uc[k]  = 1.0 + sin(2.0*PI*x/L);
    vc[k]  = 0.0;
    /* Set other arrays to 0 */
    uts1[k] = 0; uts2[k] = 0;
    vts1[k] = 0; vts2[k] = 0;
  }

  /* Loop over timesteps */ 
  while (ctime < t_F){
  
    /* Rotation factors for time-splitting scheme. */
    double cfac = cos(dt); //changed from 2*dt to dt
    double sfac = sin(dt);

    /* First substep for diffusion equation, A_1 */ 
    for (k = 0; k < nx; k++) {
      x = k*dx;
      /* Diffusion at half time step. */
      deriv = (uc[k-1] + uc[k+1] - 2*uc[k])*invdx2 ;
      uts1[k] = uc[k] + (D * deriv + vc[k])* 0.5*dt; //
      deriv = (vc[k-1] + vc[k+1] - 2*vc[k])*invdx2;
      vts1[k] = vc[k] + (D * deriv - uc[k]) * 0.5*dt;
    }

    /* Second substep for decay/growth terms, A_2 */
    for (k = 0; k < nx; k++) {
      x = k*dx;
      /* Apply rotation matrix to u and v, */
      uts2[k] = cfac*uts1[k] + sfac*vts1[k];
      vts2[k] = -sfac*uts1[k] + cfac*vts1[k];
    }

    /* Third substep for diffusion terms, A_1 */
    for (k = 0; k < nx; k++) {
      x = k*dx;
      deriv = (uts2[k-1] + uts2[k+1] - 2*uts2[k])*invdx2;
      un[k] = uts2[k] + (D * deriv + vts2[k]) * 0.5*dt;
      deriv = (vts2[k-1] + vts2[k+1] - 2*vts2[k])*invdx2;
      vn[k] = vts2[k] + (D * deriv - uts2[k]) * 0.5*dt;
    }
       
    /* Copy next values at timestep to u, v arrays. */
    memcpy(uc,un, sizeof(double) * nx);
    memcpy(vc,vn, sizeof(double) * nx);

    /* Increment time. */   
    ctime += dt;
    for (k = 0; k < nx; k++ ) {
      x = k*dx;
      printf("%g %g %g %g\n",ctime,x,uc[k],vc[k]);
    }
    
  }
  /* Free allocated memory */
  free(uc); free(un);
  free(vc); free(vn);
  free(uts1); free(uts2);
  free(vts1); free(vts2);

  return 0;
}

// The lines below don't contain any bugs! Don't modify them
void read_input(double *D, double *L, int *nx, double *t_F) {
   FILE *infile;
   if(!(infile=fopen("input.txt","r"))) {
       printf("Error opening file\n");
       exit(1);
   }
   if(4!=fscanf(infile,"%lf %lf %d %lf",D,L,nx,t_F)) {
       printf("Error reading parameters from file\n");
       exit(1);
   }
   fclose(infile);
}

So this is the code. It is meant to solve the following differential equations:

du/dt - Dd^2u/dx^2 - v = 0

dv/dt - Dd^2v/dx^2 + u = 0

It splits the equations into two parts. The second x derivative part(A1) and the decay part which contains u and v(A2) . It uses two half steps(0.5dt) for A1 and 1 full step of dt for A2. I know how to do time splitting but i dont know whether i have done it correctly here.

This is for an assignment and i have fixed all the errors and i am just trying to make the code work as intended. I have never had to solve something similar to this so i am definitely very stuck right now. The solution converges but i think its wrong. Any ideas why? Am not looking for someone to outright tell me what am doing wrong, just guide me in the right direction if you know what i mean.

PS: When i compile the code with gcc i get a warning about double main(void). Why might that be?

0

There are 0 best solutions below