1 #define sll_transformation class(sll_c_coordinate_transformation_2d_base)
12 #include "sll_memory.h"
13 #include "sll_working_precision.h"
14 #include "sll_maxwell_solvers_macros.h"
52 sll_real64,
dimension(:, :),
pointer :: vec_norm
61 sll_real64 :: eta1_min
62 sll_real64 :: eta1_max
63 sll_real64 :: eta2_min
64 sll_real64 :: eta2_max
66 sll_real64,
dimension(:),
pointer :: massmatrix
67 sll_real64,
dimension(:, :),
pointer :: dxmatrix
68 sll_real64,
dimension(:, :),
pointer :: dymatrix
77 sll_int32 :: polarization
81 sll_real64 :: eta1_min
82 sll_real64 :: eta1_max
83 sll_real64 :: delta_eta1
84 sll_real64 :: eta2_min
85 sll_real64 :: eta2_max
86 sll_real64 :: delta_eta2
87 sll_transformation,
pointer :: tau
90 sll_real64,
pointer :: f(:, :)
91 sll_real64,
pointer :: w(:, :)
92 sll_real64,
pointer :: r(:, :)
97 sll_int32 :: flux_type
122 flux_type)
result(self)
125 sll_transformation,
pointer :: tau
126 sll_int32 :: polarization
128 sll_int32,
intent(in) :: bc_east
129 sll_int32,
intent(in) :: bc_west
130 sll_int32,
intent(in) :: bc_north
131 sll_int32,
intent(in) :: bc_south
132 sll_int32,
optional :: flux_type
136 sll_allocate(self, error)
162 sll_transformation,
pointer :: tau
163 sll_int32 :: polarization
165 sll_int32,
intent(in) :: bc_east
166 sll_int32,
intent(in) :: bc_west
167 sll_int32,
intent(in) :: bc_north
168 sll_int32,
intent(in) :: bc_south
169 sll_int32,
optional :: flux_type
173 sll_real64 :: x(degree + 1)
174 sll_real64 :: w(degree + 1)
175 sll_real64 :: dlag(degree + 1, degree + 1)
176 sll_real64 :: det, dfx, dfy
177 sll_real64 :: j_mat(2, 2)
178 sll_real64 :: inv_j(2, 2)
179 sll_real64 :: dtau_ij_mat(2, 2)
180 sll_int32 :: i, j, k, l, ii, jj, kk, ll
181 sll_real64 :: xa, xb, ya, yb
186 self%bc_south = bc_south
187 self%bc_east = bc_east
188 self%bc_north = bc_north
189 self%bc_west = bc_west
191 self%nc_eta1 = self%tau%mesh%num_cells1
192 self%nc_eta2 = self%tau%mesh%num_cells2
193 self%eta1_min = self%tau%mesh%eta1_min
194 self%eta2_min = self%tau%mesh%eta2_min
195 self%eta1_max = self%tau%mesh%eta1_max
196 self%eta2_max = self%tau%mesh%eta2_max
197 self%delta_eta1 = self%tau%mesh%delta_eta1
198 self%delta_eta2 = self%tau%mesh%delta_eta2
202 self%polarization = polarization
204 if (
present(flux_type))
then
205 self%flux_type = flux_type
210 nddl = (degree + 1)*(degree + 1)
211 ncells = self%nc_eta1*self%nc_eta2
213 sll_allocate(self%cell(self%nc_eta1, self%nc_eta2), error)
219 dtau_ij_mat(1, 1) = self%tau%mesh%delta_eta1
220 dtau_ij_mat(1, 2) = 0.0_f64
221 dtau_ij_mat(2, 1) = 0.0_f64
222 dtau_ij_mat(2, 2) = self%tau%mesh%delta_eta2
224 do j = 1, self%nc_eta2
225 do i = 1, self%nc_eta1
228 i, j, degree, self%cell(i, j))
230 allocate (self%cell(i, j)%MassMatrix(nddl))
231 self%cell(i, j)%MassMatrix(nddl) = 0.0_f64
232 allocate (self%cell(i, j)%DxMatrix(nddl, nddl))
233 self%cell(i, j)%DxMatrix = 0.0_f64
234 allocate (self%cell(i, j)%DyMatrix(nddl, nddl))
235 self%cell(i, j)%DyMatrix = 0.0_f64
237 xa = self%cell(i, j)%eta1_min; xb = self%cell(i, j)%eta1_max
238 ya = self%cell(i, j)%eta2_min; yb = self%cell(i, j)%eta2_max
240 do jj = 1, degree + 1
241 do ii = 1, degree + 1
243 j_mat = tau%jacobian_matrix(xa + x(ii)*self%delta_eta1, &
244 ya + x(jj)*self%delta_eta2)
245 j_mat(:, 1) = j_mat(:, 1)*self%delta_eta1
246 j_mat(:, 2) = j_mat(:, 2)*self%delta_eta2
248 det = (j_mat(1, 1)*j_mat(2, 2) - j_mat(1, 2)*j_mat(2, 1))
250 self%cell(i, j)%MassMatrix((ii - 1)*(degree + 1) + jj) = w(ii)*w(jj)*det
255 do jj = 1, degree + 1
256 do ii = 1, degree + 1
258 do ll = 1, degree + 1
259 do kk = 1, degree + 1
261 j_mat = tau%jacobian_matrix(xa + x(kk)*self%delta_eta1, &
262 ya + x(ll)*self%delta_eta2)
263 j_mat(:, 1) = j_mat(:, 1)*self%delta_eta1
264 j_mat(:, 2) = j_mat(:, 2)*self%delta_eta2
266 det = (j_mat(1, 1)*j_mat(2, 2) - j_mat(1, 2)*j_mat(2, 1))
268 inv_j(1, 1) = j_mat(2, 2)/det
269 inv_j(1, 2) = -j_mat(1, 2)/det
270 inv_j(2, 1) = -j_mat(2, 1)/det
271 inv_j(2, 2) = j_mat(1, 1)/det
273 k = (ii - 1)*(degree + 1) + jj
274 l = (kk - 1)*(degree + 1) + ll
279 if (jj == ll) dfx = dlag(kk, ii)
280 if (ii == kk) dfy = dlag(ll, jj)
282 self%cell(i, j)%DxMatrix(k, l) = &
283 det*w(kk)*w(ll)*(inv_j(1, 1)*dfx + inv_j(2, 1)*dfy)
285 self%cell(i, j)%DyMatrix(k, l) = &
286 det*w(kk)*w(ll)*(inv_j(1, 2)*dfx + inv_j(2, 2)*dfy)
301 sll_clear_allocate(self%w((degree + 1)*(degree + 1), 4), error)
302 sll_clear_allocate(self%r((degree + 1)*(degree + 1), 4), error)
303 sll_clear_allocate(self%f((degree + 1)*(degree + 1), 4), error)
322 sll_int32 :: left, right, node, side, bc_type, flux_type
323 sll_int32 :: i, j, k, l, ii, jj, kk
327 sll_real64 :: flux(4)
328 sll_real64 :: a(4, 4)
329 sll_real64 :: a_p(4, 4)
330 sll_real64 :: a_m(4, 4)
335 do i = 1, self%nc_eta1
336 do j = 1, self%nc_eta2
338 do jj = 1, self%degree + 1
339 do ii = 1, self%degree + 1
340 kk = (ii - 1)*(self%degree + 1) + jj
341 self%w(kk, 1) = fx%array(ii, jj, i, j)
342 self%w(kk, 2) = fy%array(ii, jj, i, j)
343 self%w(kk, 3) = fz%array(ii, jj, i, j)
344 self%w(kk, 4) = self%po%array(ii, jj, i, j)
348 self%f(:, 1) = -matmul(self%cell(i, j)%DyMatrix, self%w(:, 3)) &
349 + xi*matmul(self%cell(i, j)%DxMatrix, self%w(:, 4))
350 self%f(:, 2) = matmul(self%cell(i, j)%DxMatrix, self%w(:, 3)) &
351 + xi*matmul(self%cell(i, j)%DyMatrix, self%w(:, 4))
352 self%f(:, 3) = -matmul(self%cell(i, j)%DyMatrix, self%w(:, 1)) &
353 + matmul(self%cell(i, j)%DxMatrix, self%w(:, 2))
354 self%f(:, 4) = +matmul(self%cell(i, j)%DxMatrix, self%w(:, 1)) &
355 + matmul(self%cell(i, j)%DyMatrix, self%w(:, 2))
359 print *,
"####################################################"
360 print *,
' (i,j) ', i, j
361 print
"('Ex=',9f7.3)", self%w(:, 1)
362 print
"('Ey=',9f7.3)", self%w(:, 2)
363 print
"('Bz=',9f7.3)", self%w(:, 3)
365 print
"('dEx=',9f7.3)", self%f(:, 1)
366 print
"('dEy=',9f7.3)", self%f(:, 2)
367 print
"('dBz=',9f7.3)", self%f(:, 3)
369 print
"(a)",
'side'//
' bc '//
'left'//
' w(left) ' &
370 &//
' f(left) '//
' n1 '//
' n2 '//
' flux '
375 bc_type = self%cell(i, j)%edge(side)%bc_type
376 flux_type = self%flux_type
381 l = 1 + modulo(j - 2, self%nc_eta2)
383 k = 1 + modulo(i, self%nc_eta1)
387 l = 1 + modulo(j, self%nc_eta2)
389 k = 1 + modulo(i - 2, self%nc_eta1)
393 do jj = 1, self%degree + 1
394 do ii = 1, self%degree + 1
395 kk = (ii - 1)*(self%degree + 1) + jj
396 self%r(kk, 1) = fx%array(ii, jj, k, l)
397 self%r(kk, 2) = fy%array(ii, jj, k, l)
398 self%r(kk, 3) = fz%array(ii, jj, k, l)
399 self%r(kk, 4) = self%po%array(ii, jj, k, l)
406 do node = 1, self%degree + 1
408 left =
dof_local(side, node, self%degree)
411 n1 = self%cell(i, j)%edge(side)%vec_norm(node, 1)
412 n2 = self%cell(i, j)%edge(side)%vec_norm(node, 2)
413 r = sqrt(n1*n1 + n2*n2)
415 bc_type = self%cell(i, j)%edge(side)%bc_type
417 a(1, :) = [0.0_f64, 0.0_f64, -n2, xi*n1]
418 a(2, :) = [0.0_f64, 0.0_f64, n1, xi*n2]
419 a(3, :) = [-n2, n1, 0.0_f64, 0.0_f64]
420 a(4, :) = [xi*n1, xi*n2, 0.0_f64, 0.0_f64]
422 a_p(1, :) = [(n2*n2 + xi*n1*n1)/r, n2*n1*(xi - 1._f64)/r, -n2, xi*n1]
423 a_p(2, :) = [n2*n1*(xi - 1._f64)/r, (n1*n1 + xi*n2*n2)/r, n1, xi*n2]
424 a_p(3, :) = [-n2, n1, r, 0._f64]
425 a_p(4, :) = [n1*xi, n2*xi, 0._f64, xi*r]
427 a_m(1, :) = [-(n2*n2 + xi*n1*n1)/r, -n2*n1*(xi - 1._f64)/r, -n2, xi*n1]
428 a_m(2, :) = [-n2*n1*(xi - 1._f64)/r, -(n1*n1 + xi*n2*n2)/r, n1, xi*n2]
429 a_m(3, :) = [-n2, n1, -r, 0._f64]
430 a_m(4, :) = [n1*xi, n2*xi, 0.0_f64, -xi*r]
432 select case (bc_type)
434 case (sll_p_conductor)
436 a(1, :) = [-1.0_f64, 0.0_f64, 0.0_f64, 0.0_f64]
437 a(2, :) = [0.0_f64, -1.0_f64, 0.0_f64, 0.0_f64]
438 a(3, :) = [0.0_f64, 0.0_f64, 1.0_f64, 0.0_f64]
439 a(4, :) = [-2*xi*n1, -2*xi*n2, 0.0_f64, -1.0_f64]
441 flux = matmul(a, self%w(left, :))
442 flux = matmul(a_m, flux) + matmul(a_p, self%w(left, :))
444 self%f(left, :) = self%f(left, :) - flux
446 case (sll_p_silver_muller)
448 a(1, :) = [(n2*n2 + xi*n1*n1)/r, n2*n1*(xi - 1)/r, 0.0_f64, 0.0_f64]
449 a(2, :) = [n2*n1*(xi - 1)/r, (n1*n1 + xi*n2*n2)/r, 0.0_f64, 0.0_f64]
450 a(3, :) = [0.0_f64, 0.0_f64, r, 0.0_f64]
451 a(4, :) = [0.0_f64, 0.0_f64, 0.0_f64, xi*r]
453 self%f(left, :) = self%f(left, :) - 0.5_f64*matmul(a, self%w(left, :))
459 self%f(left, :) = self%f(left, :) &
460 - 0.5_f64*matmul(a_p, self%w(left, :)) &
461 - 0.5_f64*matmul(a_m, self%r(right, :))
464 flux = 0.5_f64*(self%w(left, :) + self%r(right, :))
466 self%f(left, :) = self%f(left, :) - matmul(a, flux)
473 print
"(3i4,5f10.4)", side, bc_type, left, &
474 self%w(left, 3), self%r(right, 3), &
475 n1, n2, -n1*self%w(left, 3)
482 print
"('fEx=',9f7.3)", self%f(:, 1)
483 print
"('fEy=',9f7.3)", self%f(:, 2)
484 print
"('fBz=',9f7.3)", self%f(:, 3)
488 do jj = 1, self%degree + 1
489 do ii = 1, self%degree + 1
490 kk = (ii - 1)*(self%degree + 1) + jj
491 dx%array(ii, jj, i, j) = self%f(kk, 1)/self%cell(i, j)%MassMatrix(kk)
492 dy%array(ii, jj, i, j) = self%f(kk, 2)/self%cell(i, j)%MassMatrix(kk)
493 dz%array(ii, jj, i, j) = self%f(kk, 3)/self%cell(i, j)%MassMatrix(kk)
505 sll_int32 :: edge, dof, degree
523 sll_int32 :: edge, dof, degree
543 sll_transformation,
pointer :: tau
546 sll_real64 :: x(d + 1)
547 sll_real64 :: w(d + 1)
548 sll_real64 :: vec_norm(d + 1, 2)
549 sll_real64 :: a, b, c1, c2
551 sll_real64 :: jac_mat(2, 2)
552 sll_real64 :: co_jac_mat(2, 2)
553 sll_real64 :: dtau_ij_mat(2, 2)
554 sll_real64 :: jac_mat_sll(2, 2)
557 sll_int32 :: bc_south
559 sll_int32 :: bc_north
564 associate(lm => tau%mesh)
568 cell%eta1_min = lm%eta1_min + real(i - 1, f64)*lm%delta_eta1
569 cell%eta2_min = lm%eta2_min + real(j - 1, f64)*lm%delta_eta2
571 cell%eta1_max = cell%eta1_min + lm%delta_eta1
572 cell%eta2_max = cell%eta2_min + lm%delta_eta2
574 dtau_ij_mat(1, 1) = lm%delta_eta1
575 dtau_ij_mat(1, 2) = 0.0_f64
576 dtau_ij_mat(2, 1) = 0.0_f64
577 dtau_ij_mat(2, 2) = lm%delta_eta2
580 sll_clear_allocate(cell%edge(side)%vec_norm(1:d + 1, 1:2), error)
581 cell%edge(side)%bc_type = sll_p_interior
584 if (j == 1) cell%edge(south)%bc_type = bc_south
585 if (i == lm%num_cells1) cell%edge(east)%bc_type = bc_east
586 if (j == lm%num_cells2) cell%edge(north)%bc_type = bc_north
587 if (i == 1) cell%edge(west)%bc_type = bc_west
600 jac_mat_sll = tau%jacobian_matrix(xk, cell%eta2_min)
601 jac_mat = matmul(jac_mat_sll, dtau_ij_mat)
602 co_jac_mat(1, 1) = jac_mat(2, 2)
603 co_jac_mat(1, 2) = -jac_mat(2, 1)
604 co_jac_mat(2, 1) = -jac_mat(1, 2)
605 co_jac_mat(2, 2) = jac_mat(1, 1)
606 length = length + sqrt(jac_mat_sll(1, 1)**2 + jac_mat_sll(2, 1)**2)*wk
607 vec_norm(k, :) = matmul(co_jac_mat, (/0._f64, -1._f64/))*wk
610 cell%edge(south)%length = length
611 cell%edge(south)%vec_norm = vec_norm/lm%delta_eta1
621 jac_mat_sll = tau%jacobian_matrix(cell%eta1_max, xk)
622 jac_mat = matmul(jac_mat_sll, dtau_ij_mat)
623 co_jac_mat(1, 1) = jac_mat(2, 2)
624 co_jac_mat(1, 2) = -jac_mat(2, 1)
625 co_jac_mat(2, 1) = -jac_mat(1, 2)
626 co_jac_mat(2, 2) = jac_mat(1, 1)
627 length = length + sqrt(jac_mat_sll(1, 2)**2 + jac_mat_sll(2, 2)**2)*wk
628 vec_norm(k, :) = matmul(co_jac_mat, (/1._f64, 0._f64/))*wk
631 cell%edge(east)%length = length
632 cell%edge(east)%vec_norm = vec_norm/lm%delta_eta2
642 jac_mat_sll = tau%jacobian_matrix(xk, cell%eta2_max)
643 jac_mat = matmul(jac_mat_sll, dtau_ij_mat)
644 co_jac_mat(1, 1) = jac_mat(2, 2)
645 co_jac_mat(1, 2) = -jac_mat(2, 1)
646 co_jac_mat(2, 1) = -jac_mat(1, 2)
647 co_jac_mat(2, 2) = jac_mat(1, 1)
648 length = length + sqrt(jac_mat_sll(1, 1)**2 + jac_mat_sll(2, 1)**2)*wk
649 vec_norm(k, :) = matmul(co_jac_mat, (/0._f64, 1._f64/))*wk
652 cell%edge(north)%length = length
653 cell%edge(north)%vec_norm = vec_norm/lm%delta_eta1
663 jac_mat_sll = tau%jacobian_matrix(cell%eta1_min, xk)
664 jac_mat = matmul(jac_mat_sll, dtau_ij_mat)
665 co_jac_mat(1, 1) = jac_mat(2, 2)
666 co_jac_mat(1, 2) = -jac_mat(2, 1)
667 co_jac_mat(2, 1) = -jac_mat(1, 2)
668 co_jac_mat(2, 2) = jac_mat(1, 1)
669 length = length + sqrt(jac_mat_sll(1, 2)**2 + jac_mat_sll(2, 2)**2)*wk
670 vec_norm(k, :) = matmul(co_jac_mat, (/-1._f64, 0._f64/))*wk
673 cell%edge(west)%length = length
674 cell%edge(west)%vec_norm = vec_norm/lm%delta_eta2
Describe different boundary conditions.
Cartesian mesh basic types.
Solve Maxwell equations on cartesian domain with Disconituous Galerkine method:
subroutine, public sll_s_dg_field_2d_init(this, degree, tau, init_function)
Gauss-Lobatto integration.
function, public sll_f_gauss_lobatto_derivative_matrix(n, y)
Construction of the derivative matrix for Gauss-Lobatto, The matrix must be already allocated of size...
real(kind=f64) function, dimension(n), public sll_f_gauss_lobatto_points(n, a, b)
Returns a 1d array of size (n) containing gauss-lobatto points in the interval [a,...
real(kind=f64) function, dimension(n), public sll_f_gauss_lobatto_weights(n, a, b)
Returns a 1d array of size (n) containing gauss-lobatto weights in the interval [a,...
type(sll_t_maxwell_2d_diga) function, pointer, public sll_f_new_maxwell_2d_diga(tau, degree, polarization, bc_south, bc_east, bc_north, bc_west, flux_type)
function dof_neighbor(edge, dof, degree)
integer(kind=i32), parameter, public sll_p_uncentered
Flux parameter.
integer(kind=i32), parameter sll_centered
Flux parameter.
function dof_local(edge, dof, degree)
subroutine compute_normals(tau, bc_south, bc_east, bc_north, bc_west, i, j, d, cell)
Compute cell normals.
subroutine, public sll_s_solve_maxwell_2d_diga(self, fx, fy, fz, dx, dy, dz)
Solve the maxwell equation.
subroutine, public sll_s_maxwell_2d_diga_init(self, tau, degree, polarization, bc_south, bc_east, bc_north, bc_west, flux_type)
Initialize Maxwell solver object using DG method.
Object to describe field data with DG numerical method.
Information about a mesh cell.
Local type with edge properties.
DG method in 2D with general coordinates.