I need help configuring my code to simulate a droplet impact on a smooth surface.

[WushiDashi]

Hello everyone, I am trying to run this code on my Basilisk software to simulate a droplet impact on a smooth surface all while capturing the temperature field. only that I keep hitting road blocks, below is my code: "/*
Lit1.c - reduced-memory test variant
3D octree droplet impact (water) onto a heated smooth surface

• Reduced memory footprint for local testing:
  • MAXLEVEL = 8 (was 10)
  • coarser adapt threshold for f
  • less frequent Tecplot output / previews
• Static contact angle, temperature-dependent properties (as before)
• Compile with your existing qcc command.
  */

#include "grid/octree.h"
#include "navier-stokes/centered.h"
#include "vof.h"

#include "two-phase.h"
#include "tension.h"
#include "diffusion.h"
#include "contact.h"
#include "view.h"

#include "tecplot.h"

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

#define MAXLEVEL 8 /* reduced for memory /
#define BASELEVEL 5
#define OUTPUT_DT 0.02 / less frequent Tecplot output /
#define PREVIEW_DT 0.02
#define END_TIME 0.150
#define ADAPT_FREQ 12
#define ADAPT_THRESH_F 5e-3 / coarser refinement for f to reduce cells */
#define ADAPT_THRESH_U 0.05
#define ADAPT_THRESH_T 1.0

/* Case parameters */
double droplet_radius = 0.0015; // 1.5 mm
double droplet_diameter = 2.0 * 0.0015;
double initial_center_z = 0.0085; // 8.5 mm
double ambient_T = 293.15; // 20 C (K)
double contact_angle_static_deg = 80.0; // static contact angle (deg)
double we_target = 5.65;

/* Fields */
scalar T[]; // temperature [K]
scalar kappa[]; // curvature for output
vector height_func[]; // vector height function for contact-angle
struct OutputTec tecplot_data;

/* Air properties */
const double rho_air = 1.204;
const double mu_air = 1.82e-5;

/* Bulk liquid properties (updated periodically) */
double liquid_rho = 998.0;
double liquid_mu = 1.002e-3;
double liquid_sigma = 0.072;

/* Temperature-dependent viscosity table (°C -> Pa.s) */
static double mu_T_C[] = { 20.0, 40.0, 80.0, 120.0 };
static double mu_vals[] = { 1.002e-3, 6.53e-4, 3.55e-4, 2.00e-4 };

double interp_mu_C(double TC) {
int n = 4;
if (TC <= mu_T_C[0]) return mu_vals[0];
if (TC >= mu_T_C[n-1]) return mu_vals[n-1];
for (int i = 0; i < n-1; i++) {
if (TC >= mu_T_C[i] && TC <= mu_T_C[i+1]) {
double t = (TC - mu_T_C[i])/(mu_T_C[i+1]-mu_T_C[i]);
return mu_vals[i] * (1 - t) + mu_vals[i+1] * t;
}
}
return mu_vals[n-1];
}

double rho_water_TK(double TK) {
double TC = TK - 273.15;
return 998.2 - 0.3 * (TC - 20.0);
}

double sigma_water_TK(double TK) {
double TC = TK - 273.15;
double sigma20 = 0.072;
double dsig_dT = -1.5e-4;
double s = sigma20 + dsig_dT * (TC - 20.0);
if (s < 0.01) s = 0.01;
return s;
}

/* Update global bulk liquid properties using volume-weighted average temperature */
void update_bulk_liquid_properties() {
double sumT = 0.0, vol = 0.0;
foreach(reduction(+:sumT) reduction(+:vol)) {
if (f[] > 0.5) {
sumT += T[] * f[] * dv();
vol += f[] * dv();
}
}
if (vol > 0) {
double Tavg = sumT / vol;
double TC = Tavg - 273.15;
liquid_mu = interp_mu_C(TC);
liquid_rho = rho_water_TK(Tavg);
liquid_sigma = sigma_water_TK(Tavg);
}
rho1 = liquid_rho;
rho2 = rho_air;
mu1 = liquid_mu;
mu2 = mu_air;
f.sigma = liquid_sigma;
}

/* Estimate impact velocity from We */
double estimate_impact_velocity() {
double rho_ref = 998.0;
double sigma_ref = 0.072;
double D = droplet_diameter;
return sqrt(we_target * sigma_ref / (rho_ref * D));
}

/* Estimate contact-line radial speed near wall (z <= 2.5local_dx) /
double estimate_contact_line_speed(double local_dx) {
double cx = L0 * 0.5, cy = L0 * 0.5;
double sumu = 0.0, count = 0.0;
foreach(reduction(+:sumu) reduction(+:count)) {
double zc = z;
if (zc <= 2.5local_dx) {
if (f[] > 1e-4 && f[] < 1.0 - 1e-4) {
double rx = x - cx, ry = y - cy;
double r = sqrt(rxrx + ry*ry) + 1e-15;
double ur = (rx * u.x[] + ry * u.y[]) / r;
sumu += ur;
count += 1.0;
}
}
}
if (count > 0.0) return sumu / count;
return 0.0;
}

/* Main: set domain and N appropriate for octree /
int main() {
fprintf(stderr, "Lit1.c (small) : 3D octree droplet impact, MAXLEVEL=%d\n", MAXLEVEL);
L0 = 0.012; / 12 mm domain /
N = 1 << MAXLEVEL; / ensure N is a power-of-two */
run();
return 0;
}

/* init event: create droplet and initial impact velocity /
event init (i = 0) {
init_grid (1 << BASELEVEL); / use power-of-two base cells /
double cx = L0 * 0.5, cy = L0 * 0.5, cz = initial_center_z;
double Uimp = estimate_impact_velocity();
foreach() {
double dx = x - cx, dy = y - cy, dz = z - cz;
double r2 = dxdx + dydy + dzdz;
if (r2 <= sq(droplet_radius)) {
f[] = 1.0;
T[] = ambient_T;
u.x[] = 0.0;
u.y[] = 0.0;
u.z[] = -Uimp;
} else {
f[] = 0.0;
T[] = ambient_T;
u.x[] = u.y[] = u.z[] = 0.0;
}
}
boundary({f, T, u});
update_bulk_liquid_properties();
f.height = height_func;
double theta = contact_angle_static_deg * M_PI / 180.0;
/* use the same pattern as the example: set tangential components */
height_func.t[bottom] = contact_angle(theta);
height_func.r[bottom] = contact_angle(theta);
}

/* Thermal advection-diffusion (explicit) /
event thermal (i++) {
advection ({T}, u, dt);
const double k_water = 0.6;
const double cp_water = 4182.0;
const double k_air = 0.025;
const double cp_air = 1005.0;
foreach() {
if (f[] > 0.5) {
double alpha = k_water / (liquid_rho * cp_water);
T[] += dt * alpha * (T[1,0,0] + T[-1,0,0] + T[0,1,0] + T[0,-1,0] + T[0,0,1] + T[0,0,-1] - 6.0T[])/sq(Delta);
} else {
double alpha = k_air / (rho_air * cp_air);
T[] += dt * alpha * (T[1,0,0] + T[-1,0,0] + T[0,1,0] + T[0,-1,0] + T[0,0,1] + T[0,0,-1] - 6.0*T[])/sq(Delta);
}
}
boundary({T});
}

/* Periodically update bulk properties */
event properties_update (i += 5) {
update_bulk_liquid_properties();
}

/* Contact logging for later analysis */
event contact_logging (i += 10) {
double local_dx = L0 / (1 << MAXLEVEL);
double ucl = estimate_contact_line_speed(local_dx);
static FILE *fp = NULL;
if (!fp) fp = fopen("contact_log.txt","w");
if (fp) fprintf(fp, "%g %g %g %g\n", t, ucl, liquid_mu, liquid_sigma);
fflush(fp);
}

/* Update surface tension */
event update_sigma (i += 10) {
f.sigma = liquid_sigma;
}

/* Curvature compute */
event curvature_output (i += 20) {
curvature(f, kappa);
}

/* Adaptive mesh refinement */
event adapt (i += ADAPT_FREQ) {
scalar fe = f;
adapt_wavelet({fe, u.x, u.y, u.z, T}, (double[]){ADAPT_THRESH_F, ADAPT_THRESH_U, ADAPT_THRESH_U, ADAPT_THRESH_U, ADAPT_THRESH_T}, MAXLEVEL, BASELEVEL);
}

/* Preview image */
event preview (t += PREVIEW_DT) {
static int frame = 0;
clear();
view(fov = 20, tx = -0.5, ty = -0.1, bg = {1,1,1}, width = 800, height = 600);
draw_vof("f", color = "T", min = ambient_T, max = ambient_T + 130.0);
char title[128];
sprintf(title, "t=%.4f s mu=%.2e Pa.s sigma=%.3e N/m", t, liquid_mu, liquid_sigma);
draw_string(title, size = 20, lw = 2);
char filename[128];
sprintf(filename, "preview_%04d.ppm", frame++);
save(filename);
}

/* Tecplot output */
event tecplot_out (t += OUTPUT_DT) {
curvature(f, kappa);
scalar *list = NULL;
list = list_append(list, f);
list = list_append(list, kappa);
list = list_append(list, u.x);
list = list_append(list, u.y);
list = list_append(list, u.z);
list = list_append(list, T);
tecplot_data.tec_cc = list;
tecplot_data.num_scalars = 6;
sprintf(tecplot_data.varname, "f kappa ux uy uz T");
output_tec(tecplot_data, i, t);
free(list);
fprintf(stderr, "[TEC] step=%d t=%.6f mu=%.2e sigma=%.3e\n", i, t, liquid_mu, liquid_sigma);
}

/* Monitoring /
event monitoring (i += 50) {
stats sT = statsf(T);
stats sU = statsf(u.z);
double max_r = 0.0;
foreach(reduction(max:max_r)) {
if (f[] > 0.5) {
double rx = x - L00.5, ry = y - L00.5;
double r = sqrt(rxrx + ry*ry);
if (r > max_r) max_r = r;
}
}
double spread_factor = (2.0 * max_r) / (2.0 * droplet_radius);
double droplet_vol = 0.0;
foreach(reduction(+:droplet_vol)) {
if (f[] > 0.5) droplet_vol += f[] * dv();
}
double initial_vol = (4.0/3.0) * M_PI * pow(droplet_radius, 3.0);
double vol_ratio = droplet_vol / initial_vol;
fprintf(stderr, "i=%d t=%.4f spread=%.3f vol_ratio=%.4f T_max=%.1f u_zmax=%.3f\n",
i, t, spread_factor, vol_ratio, sT.max, sU.max);
}

/* Time step control (CFL + capillary) */
event timestep (i++) {
double umax = 0.0;
foreach(reduction(max:umax)) {
double vm = sqrt(sq(u.x[]) + sq(u.y[]) + sq(u.z[]));
if (vm > umax) umax = vm;
}
double local_dx = L0 / (1 << MAXLEVEL);
double dt_cfl = umax > 1e-15 ? 0.2 * local_dx / umax : 1e-4;
double dt_cap = 1e-4;
if (liquid_sigma > 0.0) {
double rho_local = max(liquid_rho, rho_air);
dt_cap = 0.25 * sqrt(rho_local * pow(local_dx, 3.0) / (2.0 * M_PI * liquid_sigma + 1e-20));
}
DT = dt_cfl < dt_cap ? dt_cfl : dt_cap;
if (DT < 1e-8) DT = 1e-8;
}

/* End */
event end (t = END_TIME) {
fprintf(stderr, "=== Simulation ended at t=%.6f s ===\n", t);
curvature(f, kappa);
scalar *list = NULL;
list = list_append(list, f);
list = list_append(list, kappa);
list = list_append(list, u.x);
list = list_append(list, u.y);
list = list_append(list, u.z);
list = list_append(list, T);
tecplot_data.tec_cc = list;
tecplot_data.num_scalars = 6;
sprintf(tecplot_data.varname, "f kappa ux uy uz T");
output_tec(tecplot_data, -1, t);
free(list);
dump ("final_state_small");
fprintf(stderr, "Final state saved and tecplot exported.\n");
}", the issue is that the simulation stops after 1 minute and it does not produce results
Kindly if anyone could help Ill be most appreciative