main.c

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

#define FLOAT_CMP_RTOL 1e-05
#define FLOAT_CMP_ATOL 1e-08
#define SIM_GRANULARITY 5

double x_L, x_R, x_T, x_w;
int N;
double t_f, t_D;
double X_S, X_P;
double T_C, T_w, T_0, T_1;
double P;

int is_close(double a, double b);

double X(double T);

void swap_mem();

int main(int argc, const char * argv[]) {
    // Read inputs
    FILE *input = fopen( "input.txt", "r" );
    if (input == NULL) {
        fprintf(stderr, "Can't open input file \n");
        return 1;
    }
    fscanf(input, "%lf", &x_L);
    fscanf(input, "%lf", &x_R);
    fscanf(input, "%d",  &N);
    fscanf(input, "%lf", &t_f);
    fscanf(input, "%lf", &t_D);
    fscanf(input, "%lf", &X_S);
    fscanf(input, "%lf", &X_P);
    fscanf(input, "%lf", &T_C);
    fscanf(input, "%lf", &T_w);
    fscanf(input, "%lf", &T_0);
    fscanf(input, "%lf", &T_1) ;
    fscanf(input, "%lf", &x_T);
    fscanf(input, "%lf", &x_w);
    fscanf(input, "%lf", &P);

    // Input warnings
    if (is_close(x_w, 0)) {
        fprintf(stderr, "Warning: x_w = 0 will result in division by zero");
    }
    if (is_close(T_w, 0)) {
        fprintf(stderr, "Warning: T_w = 0 will result in division by zero");
    }

    // Simulation steps
    int x_steps = N * SIM_GRANULARITY;
    int t_steps = (int) round(t_f / t_D) * SIM_GRANULARITY;

    // Simulation deltas
    double dx = (x_R - x_L) / x_steps;
    double dt = t_f / t_steps;

    // Temporal slices
    double *T, *T_;
    T = malloc(x_steps * sizeof(double));
    T_ = malloc(x_steps * sizeof(double));

    // Initialise first slice
    for (int x_n = 0; x_n < x_steps; x_n++) {
        double x = x_L + x_n * dx;
        T[x_n] = (1 + tanh(x - x_T / x_w)) / 2 * (T_1 - T_0) + T_0;
    }

    // Perform simulation
    for (int t_n = 0; t_n < t_steps; t_n++) {
        double t = t_n * dt;
        T_[0] = dt / pow(dx, 2) * (X(T[1]) + X(T[0])) * (T[1] - T[0]) + T[0];
        double _dT = T[1] - T[0];
        for (int x_n = 1; x_n < x_steps; x_n++) {
            double x = x_L + x_n * dx;
            double dT = T[x_n] - T[x_n - 1];
            double dT2 = dT - _dT;
            double dX = X(T[x_n]) - X(T[x_n - 1]);

            double old = T[x_n - 1],
                   B = dX * dT,
                   chi_old = X(T[x_n - 1]),
                   dees = dt / pow(dx, 2),
                   change = (B + chi_old * dT2) * dees;
            double new = old + change;
            T_[x_n] = new;

            if (!(t_n % SIM_GRANULARITY || x_n % SIM_GRANULARITY)) {
                printf("%lf %lf %lf\n", t, x, T[x_n]);
            }

            _dT = dT;
            swap_mem(&T, &T_);
        }
    }
    return 0;
}

void swap_mem(long *a, long *b) {
    *a ^= *b;
    *b ^= *a;
    *a ^= *b;
}

double X(double T) {
    return X_S + (X_P - X_S) * (1 + tanh(T - T_C / T_w)) / 2;
}

int is_close(double a, double b) {
    return fabs(a - b) <= (FLOAT_CMP_ATOL + FLOAT_CMP_RTOL * fabs(b));
}