From 6dd784058ae8e45b306311bd303d09a32113f50e Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Sun, 3 Apr 2016 21:56:12 +0200 Subject: PML. --- src/physics.c | 128 +++++++++++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 104 insertions(+), 24 deletions(-) (limited to 'src/physics.c') diff --git a/src/physics.c b/src/physics.c index bbe2499..952e53a 100644 --- a/src/physics.c +++ b/src/physics.c @@ -16,8 +16,8 @@ void phy_sim_destroy(phy_sim *sim) } } -static float phy_calc_fake_cond(float *fake_cond, size_t *x, size_t *widths, - size_t *dims, float dt) +static void phy_calc_fc(float *fc, size_t *x, size_t *widths, + size_t *dims, float dt) { size_t i; @@ -31,31 +31,88 @@ static float phy_calc_fake_cond(float *fake_cond, size_t *x, size_t *widths, else continue; - fake_cond[i] = EPS_ZERO / (2 * dt) * pow(f, 3); + fc[i] = EPS_ZERO / (2 * dt) * pow(f, 3); } } +static void phy_calc_uec_H(float *uec, float dt, float fc0, float fc1, + float fc2, float mu) +{ + float tmp; + + tmp = 1.0f / dt + (fc1 + fc2) / (2 * EPS_ZERO) + + fc1 * fc2 * dt / sq(2.0f * EPS_ZERO); + uec[0] = 1.0f / tmp * (1.0f / dt - (fc1 + fc2) / (2 * EPS_ZERO) + - fc1 * fc2 * dt / sq(2.0f * EPS_ZERO)); + uec[1] = -1.0f / tmp * LIGHT_SPEED / mu; + uec[2] = -1.0f / tmp * LIGHT_SPEED * dt / EPS_ZERO * fc0 / mu; + uec[3] = -1.0f / tmp * dt / sq(EPS_ZERO) * fc1 * fc2; +} + +static void phy_calc_uec_E(float *uec, float dt, float fc0, float fc1, + float fc2, float eps) +{ + float tmp; + + tmp = 1.0f / dt + (fc1 + fc2) / (2 * EPS_ZERO) + + fc1 * fc2 * dt / sq(2.0f * EPS_ZERO); + uec[0] = 1.0f / tmp * (1.0f / dt - (fc1 + fc2) / (2 * EPS_ZERO) + - fc1 * fc2 * dt / sq(2.0f * EPS_ZERO)); + uec[1] = LIGHT_SPEED / tmp; + uec[2] = 1.0f / tmp * LIGHT_SPEED * dt * fc0 / EPS_ZERO; + uec[3] = 1.0f / tmp * dt / sq(EPS_ZERO) * fc1 * fc2; + + // equations above are for the D field and E = 1/eps * D + /* + uec[0] *= 1.0f / eps; + uec[1] *= 1.0f / eps; + uec[2] *= 1.0f / eps; + uec[3] *= 1.0f / eps;*/ +} + void phy_sim_reset_aux(phy_sim *sim) { - size_t x[3], i; + size_t x[3]; phy_field_info *fi = &sim->field_info; + float dt = sim->time_delta; + + memset(sim->aux, 0, fi->size1 * sizeof(phy_field_aux_point)); for (x[2] = 1; x[2] < fi->dims[2] - 1; x[2]++) for (x[1] = 1; x[1] < fi->dims[1] - 1; x[1]++) for (x[0] = 1; x[0] < fi->dims[0] - 1; x[0]++) { phy_field_aux_point *aux; + vec3_t fc = {0, 0, 0}; aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 + x[0] * fi->xstr1; - aux->eps = 1.0f / EPS_ZERO; - aux->mu = -1.0f / MU_ZERO; + aux->eps = EPS_ZERO; + aux->mu = MU_ZERO; v3_set(aux->int_curl_E, 0, 0, 0); v3_set(aux->int_curl_H, 0, 0, 0); - phy_calc_fake_cond(aux->fake_cond, x, sim->pml_widths, - fi->dims, sim->time_delta); + phy_calc_fc(fc, x, sim->pml_widths, fi->dims, dt); + /* + for (size_t i = 0; i < 3; i++) { + aux->uec_H[i][0] = 1; + aux->uec_H[i][1] = -dt / MU_ZERO; + aux->uec_H[i][2] = 0; + aux->uec_H[i][3] = 0; + aux->uec_E[i][0] = 1; + aux->uec_E[i][1] = +dt / EPS_ZERO; + aux->uec_E[i][2] = 0; + aux->uec_E[i][3] = 0; + }*/ + + phy_calc_uec_H(aux->uec_H[0], dt, fc[0], fc[1], fc[2], 1); + phy_calc_uec_H(aux->uec_H[1], dt, fc[1], fc[0], fc[2], 1); + phy_calc_uec_H(aux->uec_H[2], dt, fc[2], fc[0], fc[1], 1); + + phy_calc_uec_E(aux->uec_E[0], dt, fc[0], fc[1], fc[2], 1); + phy_calc_uec_E(aux->uec_E[1], dt, fc[1], fc[0], fc[2], 1); + phy_calc_uec_E(aux->uec_E[2], dt, fc[2], fc[0], fc[1], 1); } } @@ -96,9 +153,15 @@ int phy_sim_create(phy_sim *sim) printf("phy_sim_create: Courant number is %f\n", 3 * LIGHT_SPEED * sim->time_delta / fi->spacing); + #if 1 sim->pml_widths[0] = fi->dims[0] / 5; sim->pml_widths[1] = fi->dims[1] / 5; sim->pml_widths[2] = fi->dims[2] / 5; + #else + sim->pml_widths[0] = 0; + sim->pml_widths[1] = 0; + sim->pml_widths[2] = 0; + #endif fi->xstr = 3; fi->ystr = fi->xstr * fi->dims[0]; @@ -207,25 +270,30 @@ void phy_sim_update_em_H(phy_sim *sim) { phy_field_info *fi = &sim->field_info; size_t x, y, z; - float *Hn, *Hp, dt; + float *Hn, *Hp; Hn = sim->fields[2].H; Hp = sim->fields[0].H; - dt = sim->time_delta; for (z = 1; z < fi->dims[2] - 1; z++) for (y = 1; y < fi->dims[1] - 1; y++) for (x = 1; x < fi->dims[0] - 1; x++) { - size_t offs; + size_t offs, i; phy_field_aux_point *aux; 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 * aux->mu * aux->curl_E[0]; - Hn[offs+1] = Hp[offs+1] + dt * aux->mu * aux->curl_E[1]; - Hn[offs+2] = Hp[offs+2] + dt * aux->mu * aux->curl_E[2]; + for (i = 0; i < 3; i++) + Hn[offs+i] = aux->uec_H[i][0] * Hp[offs+i] + + aux->uec_H[i][1] * aux->curl_E[i] + + aux->uec_H[i][2] * aux->int_curl_E[i] + + aux->uec_H[i][3] * aux->int_H[i]; + + //Hn[offs+0] = Hp[offs+0] + dt * aux->mu * aux->curl_E[0]; + //Hn[offs+1] = Hp[offs+1] + dt * aux->mu * aux->curl_E[1]; + //Hn[offs+2] = Hp[offs+2] + dt * aux->mu * aux->curl_E[2]; } } @@ -233,25 +301,30 @@ 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; + float *En, *Ep; En = sim->fields[2].E; Ep = sim->fields[0].E; - dt = sim->time_delta; for (z = 1; z < fi->dims[2] - 1; z++) for (y = 1; y < fi->dims[1] - 1; y++) for (x = 1; x < fi->dims[0] - 1; x++) { - size_t offs; + size_t offs, i; phy_field_aux_point *aux; offs = x * fi->xstr + y * fi->ystr + z * fi->zstr; aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + z * fi->zstr1; - En[offs+0] = Ep[offs+0] + dt * aux->eps * aux->curl_H[0]; - En[offs+1] = Ep[offs+1] + dt * aux->eps * aux->curl_H[1]; - En[offs+2] = Ep[offs+2] + dt * aux->eps * aux->curl_H[2]; + for (i = 0; i < 3; i++) + En[offs+i] = aux->uec_E[i][0] * Ep[offs+i] + + aux->uec_E[i][1] * aux->curl_H[i] + + aux->uec_E[i][2] * aux->int_curl_H[i] + + aux->uec_E[i][3] * aux->int_E[i]; + + //En[offs+0] = Ep[offs+0] + dt * aux->eps * aux->curl_H[0]; + //En[offs+1] = Ep[offs+1] + dt * aux->eps * aux->curl_H[1]; + //En[offs+2] = Ep[offs+2] + dt * aux->eps * aux->curl_H[2]; } } @@ -342,8 +415,9 @@ void phy_sim_debug(phy_sim *sim, size_t *coords) 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("x = [%zu %zu %zu], offs1=%zu, aux=%p\n", coords[0], coords[1], + coords[2], offs1, aux); + printf("eps = %f, mu = %f\n", aux->eps, aux->mu); printf("CE = [%E %E %E]\n", aux->curl_E[0], aux->curl_E[1], aux->curl_E[2]); printf("CH = [%E %E %E]\n", aux->curl_H[0], aux->curl_H[1], @@ -352,8 +426,14 @@ void phy_sim_debug(phy_sim *sim, size_t *coords) 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("FC = [%E %E %E]\n", aux->fake_cond[0], aux->fake_cond[1], - aux->fake_cond[2]); + printf("UEC_E = [[%E %E %E %E] [%E %E %E %E] [%E %E %E %E]]\n", + aux->uec_E[0][0], aux->uec_E[0][1], aux->uec_E[0][2], aux->uec_E[0][3], + aux->uec_E[1][0], aux->uec_E[1][1], aux->uec_E[1][2], aux->uec_E[1][3], + aux->uec_E[2][0], aux->uec_E[2][1], aux->uec_E[2][2], aux->uec_E[2][3]); + printf("UEC_H = [[%E %E %E %E] [%E %E %E %E] [%E %E %E %E]]\n", + aux->uec_H[0][0], aux->uec_H[0][1], aux->uec_H[0][2], aux->uec_H[0][3], + aux->uec_H[1][0], aux->uec_H[1][1], aux->uec_H[1][2], aux->uec_H[1][3], + aux->uec_H[2][0], aux->uec_H[2][1], aux->uec_H[2][2], aux->uec_H[2][3]); printf("=======================\n"); } -- cgit