From f1e955e7e9719aa82d1b32dac2bce62f26f75df7 Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Sat, 2 Apr 2016 21:10:03 +0200 Subject: Auxiliary fields (WIP). --- src/physics.c | 227 +++++++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 161 insertions(+), 66 deletions(-) (limited to 'src/physics.c') diff --git a/src/physics.c b/src/physics.c index 7926478..a56ad97 100644 --- a/src/physics.c +++ b/src/physics.c @@ -16,6 +16,27 @@ void phy_sim_destroy(phy_sim *sim) } } +void phy_sim_reset_aux(phy_sim *sim) +{ + size_t x, y, z; + phy_field_info *fi = &sim->field_info; + + for (z = 0; z < fi->depth; z++) + for (y = 0; y < fi->height; y++) + for (x = 0; x < fi->width; x++) { + phy_field_aux_point *point; + + point = sim->aux + z * fi->zstr1 + y * fi->ystr1 + + x * fi->xstr1; + + point->eps = 1.0f / EPS_ZERO; + point->mu = -1.0f / MU_ZERO; + + v3_set(point->int_curl_E, 0, 0, 0); + v3_set(point->int_curl_H, 0, 0, 0); + } +} + void phy_sim_reset(phy_sim *sim) { size_t i; @@ -28,9 +49,7 @@ void phy_sim_reset(phy_sim *sim) memset(sim->fields[i].H, 0, fi->size * sizeof(float)); } - memset(sim->field_eps, 0, fi->size1 * sizeof(float)); - memset(sim->field_eps, 0, fi->size1 * sizeof(float)); - phy_sim_compute_const_fields(sim); + phy_sim_reset_aux(sim); sim->time = 0; sim->running = false; @@ -49,8 +68,11 @@ int phy_sim_create(phy_sim *sim) fi->width = 100; fi->height = 100; fi->depth = 100; - fi->spacing = 1.5e10f / max3(fi->width, fi->height, fi->depth); - sim->time_delta = 0.2; + fi->spacing = 10e+9f / max3(fi->width, fi->height, fi->depth); + sim->time_delta = 0.1; + + printf("phy_sim_create: Courant number is %f\n", + 3 * LIGHT_SPEED * sim->time_delta / fi->spacing); fi->xstr = 3; fi->ystr = fi->xstr * fi->width; @@ -71,9 +93,8 @@ int phy_sim_create(phy_sim *sim) } } - sim->field_eps = malloc(fi->size1 * sizeof(float)); - sim->field_mu = malloc(fi->size1 * sizeof(float)); - if (!sim->field_eps || !sim->field_mu) { + sim->aux = malloc(fi->size1 * sizeof(phy_field_aux_point)); + if (!sim->aux) { con_printf("phy_sim_create: out of memory\n"); goto error; } @@ -92,83 +113,119 @@ error: return 1; } -void phy_sim_compute_const_fields(phy_sim *sim) +void phy_sim_update_aux_H(phy_sim *sim) { - size_t x, y, z; phy_field_info *fi = &sim->field_info; + size_t x, y, z; - for (z = 0; z < fi->depth; z++) - for (y = 0; y < fi->height; y++) - for (x = 0; x < fi->width; x++) { - float *eps, *mu; + for (z = 1; z < fi->depth - 1; z++) + for (y = 1; y < fi->height - 1; y++) + for (x = 1; x < fi->width - 1; x++) { + size_t offs, offs_px, offs_py, offs_pz; + phy_field_aux_point *aux; + float *E = sim->fields[1].E, *H = sim->fields[1].H; + + offs = x * fi->xstr + y * fi->ystr + z * fi->zstr; + offs_px = offs + fi->xstr; + offs_py = offs + fi->ystr; + offs_pz = offs + fi->zstr; - eps = sim->field_eps + z * fi->zstr1 + y * fi->ystr1 + - x * fi->xstr1; - mu = sim->field_mu + z * fi->zstr1 + y * fi->ystr1 + - x * fi->xstr1; + aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + + z * fi->zstr1; - *eps = 1.129409066837282e+11; - *mu = -7.957747154594767e+5; + aux->curl_E[0] = (E[offs_py+2] - E[offs+2] - + E[offs_pz+1] + E[offs+1]) / fi->spacing; + aux->curl_E[1] = (E[offs_pz+0] - E[offs+0] - + E[offs_px+2] + E[offs+2]) / fi->spacing; + aux->curl_E[2] = (E[offs_px+1] - E[offs+1] - + E[offs_py+0] + E[offs+0]) / fi->spacing; + + v3_acc(aux->int_H, H + offs); + v3_acc(aux->int_curl_E, aux->curl_E); } } -void phy_sim_step_E(phy_field_info *fi, float *mul_const, float dt, - float *En, float *Ep, float *H) +void phy_sim_update_aux_E(phy_sim *sim) { + phy_field_info *fi = &sim->field_info; size_t x, y, z; - float dx = fi->spacing; for (z = 1; z < fi->depth - 1; z++) for (y = 1; y < fi->height - 1; y++) for (x = 1; x < fi->width - 1; x++) { size_t offs, offs_nx, offs_ny, offs_nz; - size_t offs1; - float dExdt, dEydt, dEzdt; + phy_field_aux_point *aux; + float *E = sim->fields[1].E, *H = sim->fields[1].H; offs = x * fi->xstr + y * fi->ystr + z * fi->zstr; offs_nx = offs - fi->xstr; offs_ny = offs - fi->ystr; offs_nz = offs - fi->zstr; - offs1 = x * fi->xstr1 + y * fi->ystr1 + z * fi->zstr1; + aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + + z * fi->zstr1; - dExdt = (H[offs+2] - H[offs_ny+2] - H[offs+1] + H[offs_nz+1]); - dEydt = (H[offs+0] - H[offs_nz+0] - H[offs+2] + H[offs_nx+2]); - dEzdt = (H[offs+1] - H[offs_nx+1] - H[offs+0] + H[offs_ny+0]); + aux->curl_H[0] = (H[offs+2] - H[offs_ny+2] - + H[offs+1] + H[offs_nz+1]) / fi->spacing; + aux->curl_H[1] = (H[offs+0] - H[offs_nz+0] - + H[offs+2] + H[offs_nx+2]) / fi->spacing; + aux->curl_H[2] = (H[offs+1] - H[offs_nx+1] - + H[offs+0] + H[offs_ny+0]) / fi->spacing; - En[offs+0] = Ep[offs+0] + dt * mul_const[offs1] * dExdt / dx; - En[offs+1] = Ep[offs+1] + dt * mul_const[offs1] * dEydt / dx; - En[offs+2] = Ep[offs+2] + dt * mul_const[offs1] * dEzdt / dx; + v3_acc(aux->int_E, E + offs); + v3_acc(aux->int_curl_H, aux->curl_H); } } -void phy_sim_step_H(phy_field_info *fi, float *mul_const, float dt, - float *Hn, float *Hp, float *E) +void phy_sim_update_em_H(phy_sim *sim) { + phy_field_info *fi = &sim->field_info; size_t x, y, z; - float dx = fi->spacing; + float *Hn, *Hp, dt; + + Hn = sim->fields[2].H; + Hp = sim->fields[0].H; + dt = sim->time_delta; for (z = 1; z < fi->depth - 1; z++) for (y = 1; y < fi->height - 1; y++) for (x = 1; x < fi->width - 1; x++) { - size_t offs, offs_px, offs_py, offs_pz; - size_t offs1; - float dHxdt, dHydt, dHzdt; + size_t offs; + phy_field_aux_point *aux; offs = x * fi->xstr + y * fi->ystr + z * fi->zstr; - offs_px = offs + fi->xstr; - offs_py = offs + fi->ystr; - offs_pz = offs + fi->zstr; + aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + + z * fi->zstr1; + + Hn[offs+0] = Hp[offs+0] + dt * aux->eps * aux->curl_E[0]; + Hn[offs+1] = Hp[offs+1] + dt * aux->eps * aux->curl_E[1]; + Hn[offs+2] = Hp[offs+2] + dt * aux->eps * aux->curl_E[2]; + } +} + +void phy_sim_update_em_E(phy_sim *sim) +{ + phy_field_info *fi = &sim->field_info; + size_t x, y, z; + float *En, *Ep, dt; + + En = sim->fields[2].E; + Ep = sim->fields[0].E; + dt = sim->time_delta; - offs1 = x * fi->xstr1 + y * fi->ystr1 + z * fi->zstr1; + for (z = 1; z < fi->depth - 1; z++) + for (y = 1; y < fi->height - 1; y++) + for (x = 1; x < fi->width - 1; x++) { + size_t offs; + phy_field_aux_point *aux; - dHxdt = (E[offs_py+2] - E[offs+2] - E[offs_pz+1] + E[offs+1]); - dHydt = (E[offs_pz+0] - E[offs+0] - E[offs_px+2] + E[offs+2]); - dHzdt = (E[offs_px+1] - E[offs+1] - E[offs_py+0] + E[offs+0]); + offs = x * fi->xstr + y * fi->ystr + z * fi->zstr; + aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + + z * fi->zstr1; - Hn[offs+0] = Hp[offs+0] + dt * mul_const[offs1] * dHxdt / dx; - Hn[offs+1] = Hp[offs+1] + dt * mul_const[offs1] * dHydt / dx; - Hn[offs+2] = Hp[offs+2] + dt * mul_const[offs1] * dHzdt / dx; + En[offs+0] = Ep[offs+0] + dt * aux->mu * aux->curl_H[0]; + En[offs+1] = Ep[offs+1] + dt * aux->mu * aux->curl_H[1]; + En[offs+2] = Ep[offs+2] + dt * aux->mu * aux->curl_H[2]; } } @@ -178,13 +235,16 @@ void phy_sim_add_sources_E(phy_sim *sim) float *E; size_t x, y, z; - x = fi->width / 2; - y = fi->height / 2; - z = fi->depth / 2; + x = 0; + for (y = 1; y < fi->height - 1; y++) + for (z = 1; z < fi->depth - 1; z++) { + E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + + x * fi->xstr; - E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + x * fi->xstr; - E[0] += sin(sim->time * 0.5) * - exp(-0.1f * pow(sim->time - 3.0f, 2)) * 50; + E[0] += sin(sim->time) * 100; + E[1] += 0; + E[2] += 0; + } } void phy_sim_add_sources_H(phy_sim *sim) @@ -193,27 +253,31 @@ void phy_sim_add_sources_H(phy_sim *sim) float *H; size_t x, y, z; - x = fi->width / 2; - y = fi->height / 2; - z = fi->depth / 2; + x = 0; + for (y = 1; y < fi->height - 1; y++) + for (z = 1; z < fi->depth - 1; z++) { + H = sim->fields[2].H + z * fi->zstr + y * fi->ystr + + x * fi->xstr; - H = sim->fields[2].H + z * fi->zstr + y * fi->ystr + x * fi->xstr; - H[1] -= exp(-0.1f * pow(sim->time - 5.0f - sim->time_delta / 2, 2)) * 5; + H[0] += 0; + H[1] += 0; + H[2] += sin(sim->time + sim->time_delta / 2) / IMPEDANCE_ZERO * 100; + } } + void phy_sim_step(phy_sim *sim) { - phy_field_info *fi = &sim->field_info; phy_field_em *fields = sim->fields, tmp; int64_t start_time; start_time = get_time(); - phy_sim_step_E(fi, sim->field_mu, sim->time_delta, - fields[2].E, fields[0].E, fields[1].H); - phy_sim_add_sources_E(sim); - phy_sim_step_H(fi, sim->field_eps, sim->time_delta, - fields[2].H, fields[0].H, fields[1].E); + phy_sim_update_aux_H(sim); + phy_sim_update_em_H(sim); phy_sim_add_sources_H(sim); + phy_sim_update_aux_E(sim); + phy_sim_update_em_E(sim); + phy_sim_add_sources_E(sim); SDL_mutexP(sim->rotate_lock); @@ -230,6 +294,31 @@ void phy_sim_step(phy_sim *sim) { SDL_mutexV(sim->rotate_lock); } +void phy_sim_debug(phy_sim *sim, size_t *coords) +{ + phy_field_info *fi = &sim->field_info; + size_t offs1; + phy_field_aux_point *aux; + + offs1 = coords[2] * fi->zstr1 + coords[1] * fi->ystr1 + + coords[0] * fi->xstr1; + + aux = sim->aux + offs1; + + printf("===[ phy_sim_debug ]===\n"); + printf("x = [%zu %zu %zu], offs1=%zu\n", coords[0], coords[1], + coords[2], offs1); + printf("CE = [%E %E %E]\n", aux->curl_E[0], aux->curl_E[1], + aux->curl_E[2]); + printf("CH = [%E %E Ef]\n", aux->curl_H[0], aux->curl_H[1], + aux->curl_H[2]); + printf("ICE = [%E %E %E]\n", aux->int_curl_E[0], aux->int_curl_E[1], + aux->int_curl_E[2]); + printf("ICH = [%E %E %E]\n", aux->int_curl_H[0], aux->int_curl_H[1], + aux->int_curl_H[2]); + printf("=======================\n"); +} + void phy_run_command(phy_sim *sim, int number, void *data) { switch (number) { @@ -249,12 +338,17 @@ void phy_run_command(phy_sim *sim, int number, void *data) case PHY_CMD_RESET: phy_sim_reset(sim); break; + + case PHY_CMD_DEBUG: + phy_sim_debug(sim, data); + free(data); + break; } } int phy_thread(phy_sim *sim) { - int rv, cmd_num; + int cmd_num; void *cmd_data; while (1) { @@ -286,5 +380,6 @@ cont_main_loop: out: printf("phy_thread: quitting\n"); + return 0; } -- cgit