commit 27a869b4d3e2f7cfa02be9a858d4a5b3a3874a9b
parent 1d8e369516e373fbc0ed0bebd2dcbaad34d08e86
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 6 Apr 2020 11:54:00 +0200
Fix indentation and use consistent comment style
Diffstat:
| M | fluid.c | | | 68 | +++++++++++++++++++++++++++++++++++--------------------------------- |
| M | simulation.c | | | 21 | +++++++++++---------- |
2 files changed, 46 insertions(+), 43 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -10,9 +10,8 @@ hydrostatic_fluid_pressure_distribution()
{
int i;
for (i=0; i<sim.nz; ++i)
- sim.p_f_ghost[i+1] = sim.p_f_top +
- sim.phi[i]*sim.rho_f*sim.G*
- (sim.L_z - sim.z[i]);
+ sim.p_f_ghost[i+1] = sim.p_f_top
+ + sim.phi[i]*sim.rho_f*sim.G*(sim.L_z - sim.z[i]);
}
/* Determines the largest time step for the current simulation state. Note
@@ -52,36 +51,36 @@ set_largest_fluid_timestep(const double safety)
static double
sine_wave(const double time,
- const double amplitude,
- const double frequency,
+ const double ampl,
+ const double freq,
const double phase)
{
- return amplitude*sin(2.0*PI*frequency*time + phase);
+ return ampl*sin(2.0*PI*freq*time + phase);
}
static double
triangular_pulse(const double time,
- const double peak_amplitude,
- const double frequency,
+ const double peak_ampl,
+ const double freq,
const double peak_time)
{
- if (peak_time - 1.0/(2.0*frequency) < time && time <= peak_time)
- return peak_amplitude*2.0*frequency*(time - peak_time) + peak_amplitude;
- else if (peak_time < time && time < peak_time + 1.0/(2.0*frequency))
- return peak_amplitude*2.0*frequency*(peak_time - time) + peak_amplitude;
+ if (peak_time - 1.0/(2.0*freq) < time && time <= peak_time)
+ return peak_ampl*2.0*freq*(time - peak_time) + peak_ampl;
+ else if (peak_time < time && time < peak_time + 1.0/(2.0*freq))
+ return peak_ampl*2.0*freq*(peak_time - time) + peak_ampl;
else
return 0.0;
}
static double
square_pulse(const double time,
- const double peak_amplitude,
- const double frequency,
+ const double peak_ampl,
+ const double freq,
const double peak_time)
{
- if (peak_time - 1.0/(2.0*frequency) < time &&
- time < peak_time + 1.0/(2.0*frequency))
- return peak_amplitude;
+ if (peak_time - 1.0/(2.0*freq) < time &&
+ time < peak_time + 1.0/(2.0*freq))
+ return peak_ampl;
else
return 0.0;
}
@@ -90,7 +89,7 @@ static void
set_fluid_bcs(const double p_f_top)
{
set_bc_dirichlet(sim.p_f_ghost, sim.nz, +1, p_f_top);
- sim.p_f_ghost[sim.nz] = p_f_top; /* Include top node in BC */
+ sim.p_f_ghost[sim.nz] = p_f_top; /* Include top node in BC */
set_bc_neumann(sim.p_f_ghost,
sim.nz,
-1,
@@ -161,21 +160,24 @@ darcy_solver_1d(const int max_iter,
theta = 1.0;
if (isnan(sim.p_f_mod_pulse_time))
- p_f_top = sim.p_f_top + sine_wave(sim.t,
- sim.p_f_mod_ampl,
- sim.p_f_mod_freq,
- sim.p_f_mod_phase);
+ p_f_top = sim.p_f_top
+ + sine_wave(sim.t,
+ sim.p_f_mod_ampl,
+ sim.p_f_mod_freq,
+ sim.p_f_mod_phase);
else
if (sim.p_f_mod_pulse_shape == 1)
- p_f_top = sim.p_f_top + square_pulse(sim.t,
- sim.p_f_mod_ampl,
- sim.p_f_mod_freq,
- sim.p_f_mod_pulse_time);
+ p_f_top = sim.p_f_top
+ + square_pulse(sim.t,
+ sim.p_f_mod_ampl,
+ sim.p_f_mod_freq,
+ sim.p_f_mod_pulse_time);
else
- p_f_top = sim.p_f_top + triangular_pulse(sim.t,
- sim.p_f_mod_ampl,
- sim.p_f_mod_freq,
- sim.p_f_mod_pulse_time);
+ p_f_top = sim.p_f_top
+ + triangular_pulse(sim.t,
+ sim.p_f_mod_ampl,
+ sim.p_f_mod_freq,
+ sim.p_f_mod_pulse_time);
/* set fluid BCs (1 of 2) */
set_fluid_bcs(p_f_top);
@@ -230,14 +232,14 @@ darcy_solver_1d(const int max_iter,
#endif
p_f_ghost_new[i+1] = p_f_ghost_old[i+1]
- + epsilon*dp_f_dt_impl[i]*sim.dt;
+ + epsilon*dp_f_dt_impl[i]*sim.dt;
if (epsilon < 1.0)
p_f_ghost_new[i+1] += (1.0 - epsilon)
- *dp_f_dt_expl[i]*sim.dt;
+ *dp_f_dt_expl[i]*sim.dt;
p_f_ghost_new[i+1] = p_f_ghost_old[i+1]*(1.0 - theta)
- + p_f_ghost_new[i+1]*theta;
+ + p_f_ghost_new[i+1]*theta;
r_norm[i] = fabs((p_f_ghost_new[i+1] - sim.p_f_ghost[i+1])
/(sim.p_f_ghost[i+1] + 1e-16));
diff --git a/simulation.c b/simulation.c
@@ -290,12 +290,12 @@ compute_critical_state_friction()
if (sim.fluid)
for (i=0; i<sim.nz; ++i)
sim.mu_c[i] = sim.mu_wall/
- (sim.sigma_n_eff[i]/(sim.P_wall - sim.p_f_top));
+ (sim.sigma_n_eff[i]/(sim.P_wall - sim.p_f_top));
else
for (i=0; i<sim.nz; ++i)
sim.mu_c[i] = sim.mu_wall/
- (1.0 + (1.0 - sim.phi[i])*sim.rho_s*sim.G*
- (sim.L_z - sim.z[i])/sim.P_wall);
+ (1.0 + (1.0 - sim.phi[i])*sim.rho_s*sim.G*
+ (sim.L_z - sim.z[i])/sim.P_wall);
}
void
@@ -364,7 +364,7 @@ compute_shear_strain_rate_plastic()
int i;
for (i=0; i<sim.nz; ++i)
sim.gamma_dot_p[i] = shear_strain_rate_plastic(sim.g_ghost[i+1],
- sim.mu[i]);
+ sim.mu[i]);
}
void
@@ -372,8 +372,8 @@ compute_shear_velocity()
{
int i;
- // TODO: implement iterative solver
- // Dirichlet BC at bottom
+ /* TODO: implement iterative solver for v_x from gamma_dot_p field */
+ /* Dirichlet BC at bottom */
sim.v_x[0] = sim.v_x_bot;
for (i=1; i<sim.nz; ++i)
@@ -428,7 +428,7 @@ local_fluidity(const double p,
if (mu - C/p <= mu_s)
return 0.0;
else
- return sqrt(p/rho_s*d*d) * ((mu - C/p) - mu_s)/(b*mu);
+ return sqrt(p/rho_s*d*d)*((mu - C/p) - mu_s)/(b*mu);
}
void
@@ -497,10 +497,11 @@ poisson_solver_1d_cell_update(int i,
double coorp_term;
coorp_term = dz*dz/(2.0*pow(xi[i], 2.0));
- g_out[i+1] = 1.0/(1.0 + coorp_term)*(coorp_term*
- g_local[i] + g_in[i+2]/2.0 + g_in[i]/2.0);
+ g_out[i+1] = 1.0/(1.0 + coorp_term)
+ *(coorp_term*g_local[i] + g_in[i+2]/2.0 + g_in[i]/2.0);
- r_norm[i] = pow(g_out[i+1] - g_in[i+1], 2.0) / (pow(g_out[i+1], 2.0) + 1e-16);
+ r_norm[i] = pow(g_out[i+1] - g_in[i+1], 2.0)
+ /(pow(g_out[i+1], 2.0) + 1e-16);
#ifdef DEBUG
printf("-- %d --------------\n", i);
