//--------------------------------------------------------------------------
// untmee_2dheat_crank.c
// Solves the unsteady heat equations in a living room of rectangular domain
// with prescribed temperature on walls, windows and door of 20C and on 
// on fireplace of 80C. 
// The solution method is Crank_Nicoloson scheme, and aSOR method for solving 
// the linear system of equations is 
// Successive-Over-Relaxiation (SOR).
// Author:         Zhi-Gang Feng, http://www.mee.unt.edu/feng/
// Last Modified:  03/16/2008
//--------------------------------------------------------------------------
#include <stdio.h>
#include <math.h>
#define NXMAX 102
#define NYMAX 102
#define max_iter 2000
main()
{
  float lx,ly,lb,lfire,u_wall,u_fire,dx,dy,x[NXMAX],y[NYMAX];
  float u[NXMAX][NYMAX],unew_iter[NXMAX][NYMAX],uiter[NXMAX][NYMAX];
  float err,err_tol,omega,beta2,dt,t,tmax,alpha,r,tmp=0,rhs;
  int mx,my,nx,ny,nb,ne,i,j,iter;
  FILE *fp;
  
  // physical dimension and parameters
  lx=6.; ly=4.; lb=2; lfire=1;
  u_wall=20; u_fire=80;
  t=0; tmax=20; alpha=0.072; 

  //--------------------------------------------------------------------------
  // numerical domain
  mx=60; my=40; nx=mx+1; ny=my+1;
  dx=lx/mx; dy=ly/my; 
  beta2=(dx*dx)/(dy*dy);
  dt=0.1;
  // diffusion number; stability condition requires r<=0.25 in the case of dx=dy.
  r=alpha*dt/(dx*dx);  

  if(r>0.25){
    printf("r=%.2f is over the stability limit (r=0.25), you may have to use a smaller time step.\n",r);
  }
  else{
      printf("r=%.2f\n",r);
  }

  for(i=0;i<nx;i++){
    x[i]=i*dx;
  }

  for(j=0;j<ny;j++){
    y[j]=j*dy;
  }
  //--------------------------------------------------------------------------

  // allocate u and initize it to be u_wall
  for(i=0;i<nx;i++){
    for(j=0;j<ny;j++){
      u[i][j]=u_wall;
    }
  }

  //--------------------------------------------------------------------------
  // boundary conditions
  //left and right b.c.
  for(j=0;j<ny;j++){
    u[0][j]=u_wall;  
    u[mx][j]=u_wall; 
  }
  // fireplace
  nb=(int)(lb/dx); ne=(int)((lb+lfire)/dx);

  for(i=0;i<nb;i++){
    u[i][0]=u_wall; 
  }

  for(i=nb;i<ne;i++){
    u[i][0]=u_fire; 
  }
  
  for(i=ne;i<nx;i++){
    u[i][0]=u_wall; 
  }

  // top boundary
  for(i=0;i<nx;i++){
    u[i][my]=u_wall; 
  }
  //--------------------------------------------------------------------------
  
  err_tol=0.001; err=2*err_tol;
  iter=0; omega=1.8;
  //--------------------------------------------------------------------------
  // main loop w.r.t. time

  t=0;

  for(i=0;i<nx;i++){
    for(j=0;j<ny;j++){
      uiter[i][j]=u[i][j];
      unew_iter[i][j]=u[i][j];
    }
  }
  while (t<tmax+dt){
    
    // solving implicit equations
        


    err=2*err_tol; iter=0;
    while ((err>err_tol) && (iter<max_iter)){
      for(i=1;i<mx;i++){
	for(j=1;j<my;j++){

	  rhs=(2./r-2*(1+beta2))*u[i][j]+u[i-1][j]+u[i+1][j]+beta2*(u[i][j-1]+u[i][j+1]);
	  unew_iter[i][j]=(unew_iter[i-1][j]+uiter[i+1][j]+
			   beta2*(unew_iter[i][j-1]+uiter[i][j+1])+rhs)/(2/r+2+2*beta2);
	  unew_iter[i][j]=(1-omega)*uiter[i][j]+omega*unew_iter[i][j];
	}
      }
      
      err=0;
      for(i=0;i<nx;i++){
	for(j=0;j<ny;j++){
	  tmp=fabs(unew_iter[i][j]-uiter[i][j]);
	  if(tmp>err){
	    err=tmp;
	  }
	  uiter[i][j]=unew_iter[i][j];
	}
      }

      iter=iter+1;
      //printf("error=%f iter=%d\n",err,iter);
    }

    // move to the next time step
    for(i=1;i<mx;i++){
      for(j=1;j<my;j++){
	u[i][j]=uiter[i][j];
      }
    }

    if((int)(t/dt)%20==0) {
      printf("time=%.2f, temperature at center=%6.3f\n",t,u[(nx-1)/2][(ny-1)/2]);
    }
    t=t+dt;
  }
  
  //--------------------------------------------------------------------------
  // output the results at the final time into a data file
  fp=fopen("heat.dat","w");
  for(i=0;i<nx;i++){
    for(j=0;j<ny;j++){
      fprintf(fp,"%f  %f  %f\n",x[i],y[j],u[i][j]);
    }
    fprintf(fp,"\n");
  }
  fclose(fp);
}