FDTD_2D_Maxwell/gridtmz.c

/*
    2D FDTD simulator
    Copyright (C) 2019 Emilia Blåsten

    This program is free software: you can redistribute it and/or
    modify it under the terms of the GNU Affero General Public License
    as published by the Free Software Foundation, either version 3 of
    the License, or (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Affero General Public License for more details.

    You should have received a copy of the GNU Affero General Public
    License along with this program.  If not, see
    <http://www.gnu.org/licenses/>.
*/
/* gridtmz.c: Grid initialization function for a TMz simulation. 
 * Here the grid is simply homogeneous free space. */

#include "gridtmz.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "material.h"

void gridInitUpdateCoeff(Grid *g, int choice)
{
  if(g->type != tmZGrid) {
    perror("Update coefficient initialization error");
    fprintf(stderr, "Trying to initialize 2D coefficient on wrong type of grid.\n");
    exit(-1);
  }

  double imp0 = 377.0;
  int sizeX = g->sizeX;    // x size of domain
  int sizeY = g->sizeY;    // y size of domain
  double cdtds = g->cdtds; // Courant number
  materialInit(sizeX,sizeY,cdtds,imp0);

  int mm, nn;
  // set electric field update coefficients
  for(mm = 0; mm < sizeX; mm++) {
    for(nn = 0; nn < sizeY; nn++) {
      g->ezUpdateEcoeff[mm][nn] = ezCe(mm,nn,choice);
      g->ezUpdateHcoeff[mm][nn] = ezCh(mm,nn,choice);
    }
  }

  // set magnetic field update coefficients
  for(mm = 0; mm < sizeX; mm++) {
    for(nn = 0; nn < sizeY - 1; nn++) {
      g->hxUpdateEcoeff[mm][nn] = hxCe(mm,nn,choice);
      g->hxUpdateHcoeff[mm][nn] = hxCh(mm,nn,choice);
    }
  }
  for(mm = 0; mm < sizeX - 1; mm++) {
    for(nn = 0; nn < sizeY; nn++) {
      g->hyUpdateEcoeff[mm][nn] = hyCe(mm,nn,choice);
      g->hyUpdateHcoeff[mm][nn] = hyCh(mm,nn,choice);
    }
  }

  return;
}

void gridInit(Grid *g, int choice)
{
  g->maxTime = 30000;         // duration of simulation
  g->cdtds = 1.0 / sqrt(2.0); // Courant number

  gridInitUpdateCoeff(g, choice);

  return;
}


Grid *gridCreate(int sizeX, int sizeY) {
  // Memory for the struct itself
  Grid *g;
  g = (Grid *)calloc(1, sizeof(Grid));
  if(!g) {
    perror("gridCreate()");
    fprintf(stderr, "Failed to allocate memory for Grid.\n");
    exit(-1);
  }


  // Start allocating memory for the big arrays
  // tmZ grid!!! I.e. only the hx, hy and ez arrays allocated!!
  g->type = tmZGrid;
  g->sizeX = sizeX;
  g->sizeY = sizeY;

  // g->hx is actually an array storing the pointers to each y-array
  // This allows the use of g->hx[x][y] style notation
  g->hx = (double **)calloc(sizeX, sizeof(double*));
  g->hxCols = (double *)calloc(sizeX*(sizeY-1), sizeof(double));
  for(int m=0; m<sizeX; m++)
    g->hx[m] = g->hxCols + m*(sizeY-1);
  // the RHS above is a pointer pointing to m*(sizeY-1) steps
  // after where g->hxCols points to

  g->hxUpdateHcoeff = (double **)calloc(sizeX, sizeof(double*));
  g->hxUpdateHcoeffCols = (double *)calloc(sizeX*(sizeY-1), sizeof(double));
  for(int m=0; m<sizeX; m++)
    g->hxUpdateHcoeff[m] = g->hxUpdateHcoeffCols + m*(sizeY-1);

  g->hxUpdateEcoeff = (double **)calloc(sizeX, sizeof(double*));
  g->hxUpdateEcoeffCols = (double *)calloc(sizeX*(sizeY-1), sizeof(double));
  for(int m=0; m<sizeX; m++)
    g->hxUpdateEcoeff[m] = g->hxUpdateEcoeffCols + m*(sizeY-1);

  g->hy = (double **)calloc(sizeX-1, sizeof(double*));
  g->hyCols = (double *)calloc((sizeX-1)*sizeY, sizeof(double));
  for(int m=0; m<sizeX-1; m++)
    g->hy[m] = g->hyCols + m*sizeY;

  g->hyUpdateHcoeff = (double **)calloc(sizeX-1, sizeof(double*));
  g->hyUpdateHcoeffCols = (double *)calloc((sizeX-1)*sizeY, sizeof(double));
  for(int m=0; m<sizeX-1; m++)
    g->hyUpdateHcoeff[m] = g->hyUpdateHcoeffCols + m*sizeY;

  g->hyUpdateEcoeff = (double **)calloc(sizeX-1, sizeof(double*));
  g->hyUpdateEcoeffCols = (double *)calloc((sizeX-1)*sizeY, sizeof(double));
  for(int m=0; m<sizeX-1; m++)
    g->hyUpdateEcoeff[m] = g->hyUpdateEcoeffCols + m*sizeY;

  g->ez = (double **)calloc(sizeX, sizeof(double*));
  g->ezCols = (double *)calloc(sizeX*sizeY, sizeof(double));
  for(int m=0; m<sizeX; m++)
    g->ez[m] = g->ezCols + m*sizeY;

  g->ezUpdateHcoeff = (double **)calloc(sizeX, sizeof(double*));
  g->ezUpdateHcoeffCols = (double *)calloc(sizeX*sizeY, sizeof(double));
  for(int m=0; m<sizeX; m++)
    g->ezUpdateHcoeff[m] = g->ezUpdateHcoeffCols + m*sizeY;

  g->ezUpdateEcoeff = (double **)calloc(sizeX, sizeof(double*));
  g->ezUpdateEcoeffCols = (double *)calloc(sizeX*sizeY, sizeof(double));
  for(int m=0; m<sizeX; m++)
    g->ezUpdateEcoeff[m] = g->ezUpdateEcoeffCols + m*sizeY;

  if(!g->hx || !g->hxUpdateEcoeff || !g->hxUpdateHcoeff ||
     !g->hy || !g->hyUpdateEcoeff || !g->hyUpdateHcoeff ||
     !g->ez || !g->ezUpdateEcoeff || !g->ezUpdateHcoeff) {
    perror("gridCreate()");
    fprintf(stderr, "Failed to allocate memory for the arrays.\n");
    exit(-1);
  }

  return g;
}

void gridDestroy(Grid *g)
{
  free(g->hxCols);
  free(g->hx);
  free(g->hxUpdateHcoeffCols);
  free(g->hxUpdateHcoeff);
  free(g->hxUpdateEcoeffCols);
  free(g->hxUpdateEcoeff);

  free(g->hyCols);
  free(g->hy);
  free(g->hyUpdateHcoeffCols);
  free(g->hyUpdateHcoeff);
  free(g->hyUpdateEcoeffCols);
  free(g->hyUpdateEcoeff);

  free(g->ezCols);
  free(g->ez);
  free(g->ezUpdateHcoeffCols);
  free(g->ezUpdateHcoeff);
  free(g->ezUpdateEcoeffCols);
  free(g->ezUpdateEcoeff);

  free(g);

  return;
}