Skip to content

Instantly share code, notes, and snippets.

@buhanec
Created November 27, 2016 23:23
Show Gist options
  • Save buhanec/89f1c773c5a51718b23ee3d77cf0ecf5 to your computer and use it in GitHub Desktop.
Save buhanec/89f1c773c5a51718b23ee3d77cf0ecf5 to your computer and use it in GitHub Desktop.
#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
#define oi printf
#define OI fprintf
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 isclose(double a, double b)
{
return fabs(a - b) <= (FLOAT_CMP_ATOL + FLOAT_CMP_RTOL * fabs(b));
}
double X(double T)
{
return X_S + (X_P - X_S) * (1 + tanh(T - T_C / T_w)) / 2;
}
double* diffs(double *a, double l_padding, double r_padding, int length, int offset)
{
double *r;
r = malloc(length * sizeof(double));
int i = 0;
for (; i < offset; i++) {
r[i] = l_padding;
}
for (; (offset + i) < length - 1; i++) {
r[i + offset] = a[i] - a[i + 1];
}
for (; i < length; i++) {
r[i] = r_padding;
}
return r;
}
void swap_mem();
int main(int argc, const char * argv[])
{
FILE *input = fopen( "input.txt", "r" );
if (input == NULL) {
OI(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);
if (isclose(x_w, 0)) {
OI(stderr, "Warning: x_w = 0 will result in division by zero");
}
if (isclose(T_w, 0)) {
OI(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));
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;
}
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]);
T_[x_n] = T[x_n - 1] + dX / dx * dT / dx + X(T[x_n - 1]) * dT2 / pow(dx, 2);
if (!(t_n % SIM_GRANULARITY || x_n % SIM_GRANULARITY)) {
oi("%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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment