Created
November 27, 2016 23:23
-
-
Save buhanec/89f1c773c5a51718b23ee3d77cf0ecf5 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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