Skip to content

Commit ad33264

Browse files
lmdiazangulolmdiazangulo
authored
Merge pull request #117 from OpenSEMBA/feature/sgbc
Feature/sgbc
2 parents 8c4ca57 + ca2ed1f commit ad33264

File tree

7 files changed

+143745
-12
lines changed

7 files changed

+143745
-12
lines changed

‎src_main_pub/healer.F90‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -222,7 +222,7 @@ SUBROUTINE CreateVolumeMM (layoutnumber, Mtag, tags, numertag, MMiEx, MMiEy, MMi
222222
DO i = punto%XI, puntoPlus1%XE
223223
medio = MMiHx (i, j, k)
224224
! IF (medio /= 0) THEN !ojo esto estaba antes de 031016 y daba maxima prioridad al medio 0 PEC. Ahora puedo tener medios con mas prioridad!!! !?!? cambio agresivo 031016!!!
225-
IF (med(indicemedio)%Priority >= med(medio)%Priority) THEN
225+
IF (med(indicemedio)%Priority > med(medio)%Priority) THEN
226226
MMiHx (i, j, k) = indicemedio
227227
Mtag(i,j,k)=64*numertag
228228
tags%face%x(i,j,k) = 64*numertag
@@ -238,7 +238,7 @@ SUBROUTINE CreateVolumeMM (layoutnumber, Mtag, tags, numertag, MMiEx, MMiEy, MMi
238238
DO i = punto%XI, punto%XE
239239
medio = MMiHy (i, j, k)
240240
! IF (medio /= 0) THEN !ojo esto estaba antes de 031016 y daba maxima prioridad al medio 0 PEC. Ahora puedo tener medios con mas prioridad!!! !?!? cambio agresivo 031016!!!
241-
IF (med(indicemedio)%Priority >= med(medio)%Priority) THEN
241+
IF (med(indicemedio)%Priority > med(medio)%Priority) THEN
242242
MMiHy (i, j, k) = indicemedio
243243
Mtag(i,j,k)=64*numertag
244244
tags%face%y(i,j,k) = 64*numertag
@@ -254,7 +254,7 @@ SUBROUTINE CreateVolumeMM (layoutnumber, Mtag, tags, numertag, MMiEx, MMiEy, MMi
254254
DO i = punto%XI, punto%XE
255255
medio = MMiHz (i, j, k)
256256
! IF (medio /= 0) THEN !ojo esto estaba antes de 031016 y daba maxima prioridad al medio 0 PEC. Ahora puedo tener medios con mas prioridad!!! !?!? cambio agresivo 031016!!!
257-
IF (med(indicemedio)%Priority >= med(medio)%Priority) THEN
257+
IF (med(indicemedio)%Priority > med(medio)%Priority) THEN
258258
MMiHz (i, j, k) = indicemedio
259259
Mtag(i,j,k)=64*numertag
260260
tags%face%z(i,j,k) = 64*numertag
@@ -361,12 +361,12 @@ SUBROUTINE CreateSurfaceMM (layoutnumber, Mtag, tags, numertag, MMiEx, MMiEy, MM
361361
CALL AddToShared (iEz, i, j, k, indicemedio, medio, Eshared)
362362
END IF
363363
END DO
364-
END DO
364+
END DO
365365
DO j = punto%YI, punto%YE
366366
DO k = punto%ZI, punto%ZE
367367
medio = MMiHx (i, j, k)
368368
! IF (medio /= 0) THEN !ojo esto estaba antes de 031016 y daba maxima prioridad al medio 0 PEC. Ahora puedo tener medios con mas prioridad!!! !?!? cambio agresivo 031016!!!
369-
IF (med(indicemedio)%Priority >= med(medio)%Priority) then
369+
IF (med(indicemedio)%Priority > med(medio)%Priority) then
370370
MMiHx (i, j, k) = indicemedio;
371371
Mtag(i,j,k)=64*numertag ! if (.true..or.(Mtag(i,j,k)==0).or.(int(Mtag(i,j,k)/64) == numertag)) Mtag(i,j,k) = IBSET(64*numertag,3);
372372
tags%face%x(i,j,k) = 64*numertag
@@ -408,7 +408,7 @@ SUBROUTINE CreateSurfaceMM (layoutnumber, Mtag, tags, numertag, MMiEx, MMiEy, MM
408408
DO k = punto%ZI, punto%ZE
409409
medio = MMiHy (i, j, k)
410410
! IF (medio /= 0) THEN !ojo esto estaba antes de 031016 y daba maxima prioridad al medio 0 PEC. Ahora puedo tener medios con mas prioridad!!! !?!? cambio agresivo 031016!!!
411-
IF (med(indicemedio)%Priority >= med(medio)%Priority) then
411+
IF (med(indicemedio)%Priority > med(medio)%Priority) then
412412
MMiHy (i, j, k) = indicemedio;
413413
Mtag(i,j,k)=64*numertag ! if (.true..or.(Mtag(i,j,k)==0).or.(int(Mtag(i,j,k)/64) == numertag)) Mtag(i,j,k) = IBSET(64*numertag,4);;
414414
tags%face%y(i,j,k) = 64*numertag
@@ -449,7 +449,7 @@ SUBROUTINE CreateSurfaceMM (layoutnumber, Mtag, tags, numertag, MMiEx, MMiEy, MM
449449
DO j = punto%YI, punto%YE
450450
medio = MMiHz (i, j, k)
451451
! IF (medio /= 0) THEN !ojo esto estaba antes de 031016 y daba maxima prioridad al medio 0 PEC. Ahora puedo tener medios con mas prioridad!!! !?!? cambio agresivo 031016!!!
452-
IF (med(indicemedio)%Priority >= med(medio)%Priority) then
452+
IF (med(indicemedio)%Priority > med(medio)%Priority) then
453453
MMiHz (i, j, k) = indicemedio;
454454
Mtag(i,j,k)=64*numertag ! if (.true..or.(Mtag(i,j,k)==0).or.(int(Mtag(i,j,k)/64) == numertag)) Mtag(i,j,k) = IBSET(64*numertag,5);
455455
tags%face%z(i,j,k) = 64*numertag

‎test/pyWrapper/test_full_system.py‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -268,6 +268,75 @@ def test_sgbc_structured_resistance(tmp_path):
268268
assert np.allclose(i.array[-101:-1], np.ones(100)*i.array[-100], rtol=1e-3)
269269
assert np.allclose(-1/i.array[-101:-1], np.ones(100)*(50+45), rtol=0.05)
270270

271+
272+
def test_pec_overlapping_sgbcs(tmp_path):
273+
""" Test that PEC surfaces overlapping SGBC surfaces prioritize PEC.
274+
"""
275+
fn = CASES_FOLDER + 'sgbcOverlapping/sgbcOverlapping.fdtd.json'
276+
277+
# Runs case without overlap.
278+
solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path)
279+
solver.run()
280+
p = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0])
281+
t = p['time'].to_numpy()
282+
iSGBC = p['current'].to_numpy()
283+
284+
# Adds current SGBC elements as PEC. Now both are defined over same surface.
285+
sgbcElementIds = solver["materialAssociations"][1]["elementIds"]
286+
solver['materialAssociations'][0]["elementIds"].extend(sgbcElementIds)
287+
solver.cleanUp()
288+
solver.run()
289+
iPEC = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0])['current'].to_numpy()
290+
291+
292+
# For debugging only.
293+
# plt.figure()
294+
# plt.plot(t, iSGBC,'.-', label='SGBC case')
295+
# plt.plot(t, iPEC,'.-', label='PEC overlapping')
296+
# plt.grid(which='both')
297+
# plt.legend()
298+
# plt.show()
299+
300+
301+
# Checks values are different due to PEC prioritization.
302+
assert np.all(np.greater(np.abs(iPEC[1000:]), np.abs(iSGBC[1000:])))
303+
304+
def test_sgbc_overlapping_sgbc(tmp_path):
305+
""" Test that SGBC surfaces overlapping SGBC surfaces prioritize first in MatAss.
306+
"""
307+
fn = CASES_FOLDER + 'sgbcOverlapping/sgbcOverlapping.fdtd.json'
308+
309+
# Runs case without overlap.
310+
solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path)
311+
# Changes materialId in first SGBC in MatAss to material with larger conductivity.
312+
solver['materialAssociations'][1]["materialId"] = 6
313+
solver.cleanUp()
314+
solver.run()
315+
p = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0])
316+
317+
t = p['time'].to_numpy()
318+
iSGBC_top = p['current'].to_numpy()
319+
320+
# Changes materialId in second SGBC in MatAss to material with larger conductivity.
321+
solver['materialAssociations'][1]["materialId"] = 2
322+
solver['materialAssociations'][2]["materialId"] = 6
323+
solver.cleanUp()
324+
solver.run()
325+
iSGBC_bottom = Probe(solver.getSolvedProbeFilenames("Bulk probe")[0])['current'].to_numpy()
326+
327+
328+
# For debugging only.
329+
# plt.figure()
330+
# plt.plot(t, iSGBC_top,'.-', label='SGBC sigma = 40 S/m, top')
331+
# plt.plot(t, iSGBC_bottom,'.-', label='SGBC sigma = 20 S/m, bottom')
332+
# plt.grid(which='both')
333+
# plt.legend()
334+
# plt.show()
335+
336+
337+
# Checks values are different due to prioritization of first written.
338+
assert np.all(np.greater(np.abs(iSGBC_top[1000:]), np.abs(iSGBC_bottom[1000:])))
339+
271340
def test_dielectric_transmission(tmp_path):
272341
_FIELD_TOLERANCE = 0.05
273342

0 commit comments

Comments
 (0)