summaryrefslogtreecommitdiff
path: root/src/physics.c
diff options
context:
space:
mode:
authorPaweł Redman <pawel.redman@gmail.com>2016-04-03 21:56:12 +0200
committerPaweł Redman <pawel.redman@gmail.com>2016-04-03 21:56:12 +0200
commit6dd784058ae8e45b306311bd303d09a32113f50e (patch)
treee730668b9b1bd212898bb8f36c670e516fec438c /src/physics.c
parent3ef209bad0cd7c1ae83a65702633e43f5a91c97f (diff)
PML.
Diffstat (limited to 'src/physics.c')
-rw-r--r--src/physics.c128
1 files changed, 104 insertions, 24 deletions
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");
}