Compute fake conductivity (for PML).

1 files changed, 89 insertions, 51 deletions

diff --git a/src/physics.c b/src/physics.c
index a56ad97..bbe2499 100644
--- a/src/physics.c
+++ b/src/physics.c
@@ -1,4 +1,4 @@
-#include "common.h"
+#include "common.h"                                 
 #include "itc.h"
 #include "physics.h"
 
@@ -16,24 +16,46 @@ 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)
+{
+	size_t i;
+
+	for (i = 0; i < 3; i++) {
+		float f;
+
+		if (x[i] <= widths[i])
+			f = 1.0f - (float)(x[i] - 1) / widths[i];
+		else if (x[i] >= dims[i] - widths[i] - 1)
+			f = (float)(x[i] + widths[i] - dims[i] + 2) / widths[i];
+		else
+			continue;
+
+		fake_cond[i] = EPS_ZERO / (2 * dt) * pow(f, 3);
+	}
+}
+
 void phy_sim_reset_aux(phy_sim *sim)
 {
-	size_t x, y, z;
+	size_t x[3], i;
 	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;
+	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;
+
+		aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 +
+			x[0] * fi->xstr1;
 
-		point = sim->aux + z * fi->zstr1 + y * fi->ystr1 +
-			x * fi->xstr1;
+		aux->eps = 1.0f / EPS_ZERO;
+		aux->mu = -1.0f / MU_ZERO;
 
-		point->eps = 1.0f / EPS_ZERO;
-		point->mu = -1.0f / MU_ZERO;
+		v3_set(aux->int_curl_E, 0, 0, 0);
+		v3_set(aux->int_curl_H, 0, 0, 0);
 
-		v3_set(point->int_curl_E, 0, 0, 0);
-		v3_set(point->int_curl_H, 0, 0, 0);
+		phy_calc_fake_cond(aux->fake_cond, x, sim->pml_widths,
+		                   fi->dims, sim->time_delta);
 	}
 }
 
@@ -65,24 +87,28 @@ int phy_sim_create(phy_sim *sim)
 
 	fi = &sim->field_info;
 
-	fi->width = 100;
-	fi->height = 100;
-	fi->depth = 100;
-	fi->spacing = 10e+9f / max3(fi->width, fi->height, fi->depth);
+	fi->dims[0] = 50;
+	fi->dims[1] = 50;
+	fi->dims[2] = 50;
+	fi->spacing = 10e+9f / max3(fi->dims[0], fi->dims[1], fi->dims[2]);
 	sim->time_delta = 0.1;
 
 	printf("phy_sim_create: Courant number is %f\n",
 	       3 * LIGHT_SPEED * sim->time_delta / fi->spacing);
 
+	sim->pml_widths[0] = fi->dims[0] / 5;
+	sim->pml_widths[1] = fi->dims[1] / 5;
+	sim->pml_widths[2] = fi->dims[2] / 5;
+
 	fi->xstr = 3;
-	fi->ystr = fi->xstr * fi->width;
-	fi->zstr = fi->ystr * fi->height;
-	fi->size = fi->zstr * fi->depth;
+	fi->ystr = fi->xstr * fi->dims[0];
+	fi->zstr = fi->ystr * fi->dims[1];
+	fi->size = fi->zstr * fi->dims[2];
 
 	fi->xstr1 = 1;
-	fi->ystr1 = fi->xstr1 * fi->width;
-	fi->zstr1 = fi->ystr1 * fi->height;
-	fi->size1 = fi->zstr1 * fi->depth;
+	fi->ystr1 = fi->xstr1 * fi->dims[0];
+	fi->zstr1 = fi->ystr1 * fi->dims[1];
+	fi->size1 = fi->zstr1 * fi->dims[2];
 
 	for (i = 0; i < 3; i++) { 
 		sim->fields[i].E = malloc(fi->size * sizeof(float));
@@ -118,9 +144,9 @@ void phy_sim_update_aux_H(phy_sim *sim)
 	phy_field_info *fi = &sim->field_info;
 	size_t x, y, z;
 
-	for (z = 1; z < fi->depth - 1; z++)
-	for (y = 1; y < fi->height - 1; y++)
-	for (x = 1; x < fi->width - 1; x++) {
+	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, offs_px, offs_py, offs_pz;
 		phy_field_aux_point *aux;
 		float *E = sim->fields[1].E, *H = sim->fields[1].H;
@@ -150,9 +176,9 @@ void phy_sim_update_aux_E(phy_sim *sim)
 	phy_field_info *fi = &sim->field_info;
 	size_t x, y, z;
 
-	for (z = 1; z < fi->depth - 1; z++)
-	for (y = 1; y < fi->height - 1; y++)
-	for (x = 1; x < fi->width - 1; x++) {
+	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, offs_nx, offs_ny, offs_nz;
 		phy_field_aux_point *aux;
 		float *E = sim->fields[1].E, *H = sim->fields[1].H;
@@ -187,19 +213,19 @@ void phy_sim_update_em_H(phy_sim *sim)
 	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++) {
+	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;
 		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;
+		      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];
+		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];
 	}
 }
 
@@ -213,9 +239,9 @@ void phy_sim_update_em_E(phy_sim *sim)
 	Ep = sim->fields[0].E;
 	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++) {
+	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;
 		phy_field_aux_point *aux;
 
@@ -223,9 +249,9 @@ void phy_sim_update_em_E(phy_sim *sim)
 		aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 +
 		         z * fi->zstr1;
 
-		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];
+		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];
 	}
 }
 
@@ -235,34 +261,44 @@ void phy_sim_add_sources_E(phy_sim *sim)
 	float *E;
 	size_t x, y, z;
 
-	x = 0;
-	for (y = 1; y < fi->height - 1; y++)
-	for (z = 1; z < fi->depth - 1; z++) {
+
+	x = fi->dims[0]/2;
+	y = fi->dims[1]/2;
+	z = fi->dims[2]/2;
+	E = sim->fields[2].E + z * fi->zstr + y * fi->ystr +
+		     x * fi->xstr;
+
+	E[0] += exp(-40.0f * pow(sim->time - 3, 2)) * 50;
+
+/*
+	x = 10;
+	for (y = 1; y < fi->dims[1] - 1; y++)
+	for (z = 1; z < fi->dims[2] - 1; z++) {
 		E = sim->fields[2].E + z * fi->zstr + y * fi->ystr +
 		     x * fi->xstr;
 
-		E[0] += sin(sim->time) * 100;
+		E[0] += sin(sim->time * 3) * 1;
 		E[1] += 0;
 		E[2] += 0;
-	}
+	}*/
 }
 
 void phy_sim_add_sources_H(phy_sim *sim)
 {
-	phy_field_info *fi = &sim->field_info;
+	/*phy_field_info *fi = &sim->field_info;
 	float *H;
 	size_t x, y, z;
 
 	x = 0;
-	for (y = 1; y < fi->height - 1; y++)
-	for (z = 1; z < fi->depth - 1; z++) {
+	for (y = 1; y < fi->dims[1] - 1; y++)
+	for (z = 1; z < fi->dims[2] - 1; z++) {
 		H = sim->fields[2].H + z * fi->zstr + y * fi->ystr +
 		     x * fi->xstr;
 
 		H[0] += 0;
 		H[1] += 0;
 		H[2] += sin(sim->time + sim->time_delta / 2)  / IMPEDANCE_ZERO * 100;
-	}
+	}*/
 }
 
 
@@ -310,12 +346,14 @@ void phy_sim_debug(phy_sim *sim, size_t *coords)
 	       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],
+	printf("CH = [%E %E %E]\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("FC = [%E %E %E]\n", aux->fake_cond[0], aux->fake_cond[1],
+	       aux->fake_cond[2]);
 	printf("=======================\n");
 }