Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_maxwell_2d_diga.F90
Go to the documentation of this file.
1 #define sll_transformation class(sll_c_coordinate_transformation_2d_base)
2 
10 
11 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12 #include "sll_memory.h"
13 #include "sll_working_precision.h"
14 #include "sll_maxwell_solvers_macros.h"
15 
17  sll_p_conductor, &
18  sll_p_interior, &
19  sll_p_silver_muller
20 
21  use sll_m_cartesian_meshes, only: &
23 
26 
27  use sll_m_dg_fields, only: &
30 
35 
36  implicit none
37 
38  public :: &
44 
45  private
46 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
47 
49  type :: edge_type
50 
51  sll_real64 :: length
52  sll_real64, dimension(:, :), pointer :: vec_norm
53  sll_int32 :: bc_type
54 
55  end type edge_type
56 
58  type :: cell_type
59 
60  sll_int32 :: i, j
61  sll_real64 :: eta1_min
62  sll_real64 :: eta1_max
63  sll_real64 :: eta2_min
64  sll_real64 :: eta2_max
65  type(edge_type) :: edge(4)
66  sll_real64, dimension(:), pointer :: massmatrix
67  sll_real64, dimension(:, :), pointer :: dxmatrix
68  sll_real64, dimension(:, :), pointer :: dymatrix
69 
70  end type cell_type
71 
74  private
75  sll_int32 :: nc_eta1
76  sll_int32 :: nc_eta2
77  sll_int32 :: polarization
78  sll_real64 :: e_0
79  sll_real64 :: mu_0
80  sll_real64 :: c
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
88  sll_int32 :: degree
89  type(cell_type), pointer :: cell(:, :)
90  sll_real64, pointer :: f(:, :)
91  sll_real64, pointer :: w(:, :)
92  sll_real64, pointer :: r(:, :)
93  sll_int32 :: bc_south
94  sll_int32 :: bc_east
95  sll_int32 :: bc_north
96  sll_int32 :: bc_west
97  sll_int32 :: flux_type
98  type(sll_t_dg_field_2d) :: po
99  sll_real64 :: xi
100 
101  contains
102 
103  procedure :: init => sll_s_maxwell_2d_diga_init
104  procedure :: solve => sll_s_solve_maxwell_2d_diga
105 
106  end type sll_t_maxwell_2d_diga
107 
109  sll_int32, parameter :: sll_centered = 20
111  sll_int32, parameter :: sll_p_uncentered = 21
112 
113 contains
114 
116  degree, &
117  polarization, &
118  bc_south, &
119  bc_east, &
120  bc_north, &
121  bc_west, &
122  flux_type) result(self)
123 
124  type(sll_t_maxwell_2d_diga), pointer :: self
125  sll_transformation, pointer :: tau
126  sll_int32 :: polarization
127  sll_int32 :: degree
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
133 
134  sll_int32 :: error
135 
136  sll_allocate(self, error)
137 
138  call sll_s_maxwell_2d_diga_init(self, &
139  tau, &
140  degree, &
141  polarization, &
142  bc_south, &
143  bc_east, &
144  bc_north, &
145  bc_west, &
146  flux_type)
147 
148  end function sll_f_new_maxwell_2d_diga
149 
151  subroutine sll_s_maxwell_2d_diga_init(self, &
152  tau, &
153  degree, &
154  polarization, &
155  bc_south, &
156  bc_east, &
157  bc_north, &
158  bc_west, &
159  flux_type)
160 
161  class(sll_t_maxwell_2d_diga) :: self
162  sll_transformation, pointer :: tau
163  sll_int32 :: polarization
164  sll_int32 :: degree
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
170 
171  sll_int32 :: nddl
172  sll_int32 :: ncells
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
182 
183  sll_int32 :: error
184 
185  self%tau => tau
186  self%bc_south = bc_south
187  self%bc_east = bc_east
188  self%bc_north = bc_north
189  self%bc_west = bc_west
190 
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
199 
200  self%xi = 0.0_f64
201  self%degree = degree
202  self%polarization = polarization
203 
204  if (present(flux_type)) then
205  self%flux_type = flux_type
206  else
207  self%flux_type = sll_centered
208  end if
209 
210  nddl = (degree + 1)*(degree + 1)
211  ncells = self%nc_eta1*self%nc_eta2
212 
213  sll_allocate(self%cell(self%nc_eta1, self%nc_eta2), error)
214 
215  x = sll_f_gauss_lobatto_points(degree + 1, 0.0_f64, 1.0_f64)
216  w = sll_f_gauss_lobatto_weights(degree + 1, 0.0_f64, 1.0_f64)
217  dlag = sll_f_gauss_lobatto_derivative_matrix(degree + 1, x)
218 
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
223 
224  do j = 1, self%nc_eta2 !Loop over cells
225  do i = 1, self%nc_eta1
226 
227  call compute_normals(tau, bc_south, bc_east, bc_north, bc_west, &
228  i, j, degree, self%cell(i, j))
229 
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
236 
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
239 
240  do jj = 1, degree + 1
241  do ii = 1, degree + 1
242 
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
247 
248  det = (j_mat(1, 1)*j_mat(2, 2) - j_mat(1, 2)*j_mat(2, 1))
249 
250  self%cell(i, j)%MassMatrix((ii - 1)*(degree + 1) + jj) = w(ii)*w(jj)*det
251 
252  end do
253  end do
254 
255  do jj = 1, degree + 1
256  do ii = 1, degree + 1
257 
258  do ll = 1, degree + 1
259  do kk = 1, degree + 1
260 
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
265 
266  det = (j_mat(1, 1)*j_mat(2, 2) - j_mat(1, 2)*j_mat(2, 1))
267 
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
272 
273  k = (ii - 1)*(degree + 1) + jj
274  l = (kk - 1)*(degree + 1) + ll
275 
276  dfx = 0.0_f64
277  dfy = 0.0_f64
278 
279  if (jj == ll) dfx = dlag(kk, ii) !Here we use transpose
280  if (ii == kk) dfy = dlag(ll, jj) !dlag (don't know why)
281 
282  self%cell(i, j)%DxMatrix(k, l) = &
283  det*w(kk)*w(ll)*(inv_j(1, 1)*dfx + inv_j(2, 1)*dfy)
284 
285  self%cell(i, j)%DyMatrix(k, l) = &
286  det*w(kk)*w(ll)*(inv_j(1, 2)*dfx + inv_j(2, 2)*dfy)
287 
288  end do
289  end do
290 
291  end do
292  end do
293 
294  end do
295  end do
296 
297  !call sll_o_display(self%cell(1,1)%MassMatrix,"f9.4")
298  !call sll_o_display(self%cell(1,1)%DxMatrix,"f9.4")
299  !call sll_o_display(self%cell(1,1)%DyMatrix,"f9.4")
300 
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)
304 
305  call sll_s_dg_field_2d_init(self%po, degree, tau)
306 
307  end subroutine sll_s_maxwell_2d_diga_init
308 
310  subroutine sll_s_solve_maxwell_2d_diga(self, fx, fy, fz, dx, dy, dz)
311 
312  class(sll_t_maxwell_2d_diga) :: self
313 
314  type(sll_t_dg_field_2d) :: fx
315  type(sll_t_dg_field_2d) :: fy
316  type(sll_t_dg_field_2d) :: fz
317 
318  type(sll_t_dg_field_2d) :: dx
319  type(sll_t_dg_field_2d) :: dy
320  type(sll_t_dg_field_2d) :: dz
321 
322  sll_int32 :: left, right, node, side, bc_type, flux_type
323  sll_int32 :: i, j, k, l, ii, jj, kk
324  sll_real64 :: n1
325  sll_real64 :: n2
326  sll_real64 :: r
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)
331  sll_real64 :: xi
332 
333  xi = self%xi
334 
335  do i = 1, self%nc_eta1
336  do j = 1, self%nc_eta2
337 
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)
345  end do
346  end do
347 
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))
356 
357 #ifdef VERBOSE
358 
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)
364 
365  print"('dEx=',9f7.3)", self%f(:, 1)
366  print"('dEy=',9f7.3)", self%f(:, 2)
367  print"('dBz=',9f7.3)", self%f(:, 3)
368 
369  print"(a)", 'side'//' bc '//'left'//' w(left) ' &
370  &//' f(left) '//' n1 '//' n2 '//' flux '
371 #endif
372 
373  do side = 1, 4 ! Loop over each side of the cell
374 
375  bc_type = self%cell(i, j)%edge(side)%bc_type
376  flux_type = self%flux_type
377  !periodic boundary conditions
378  select case (side)
379  case (south)
380  k = i
381  l = 1 + modulo(j - 2, self%nc_eta2)
382  case (east)
383  k = 1 + modulo(i, self%nc_eta1)
384  l = j
385  case (north)
386  k = i
387  l = 1 + modulo(j, self%nc_eta2)
388  case (west)
389  k = 1 + modulo(i - 2, self%nc_eta1)
390  l = j
391  end select
392 
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)
400  end do
401  end do
402 
403 #ifdef VERBOSE
404  print *, '--'
405 #endif
406  do node = 1, self%degree + 1
407 
408  left = dof_local(side, node, self%degree)
409  right = dof_neighbor(side, node, self%degree)
410 
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)
414 
415  bc_type = self%cell(i, j)%edge(side)%bc_type
416 
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]
421 
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]
426 
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]
431 
432  select case (bc_type)
433 
434  case (sll_p_conductor)
435 
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]
440 
441  flux = matmul(a, self%w(left, :))
442  flux = matmul(a_m, flux) + matmul(a_p, self%w(left, :))
443 
444  self%f(left, :) = self%f(left, :) - flux
445 
446  case (sll_p_silver_muller)
447 
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]
452 
453  self%f(left, :) = self%f(left, :) - 0.5_f64*matmul(a, self%w(left, :))
454 
455  case default
456 
457  if (self%flux_type == sll_p_uncentered) then
458 
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, :))
462  else
463 
464  flux = 0.5_f64*(self%w(left, :) + self%r(right, :))
465 
466  self%f(left, :) = self%f(left, :) - matmul(a, flux)
467 
468  end if
469 
470  end select
471 
472 #ifdef VERBOSE
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)
476 #endif
477  end do
478 
479  end do
480 
481 #ifdef VERBOSE
482  print"('fEx=',9f7.3)", self%f(:, 1)
483  print"('fEy=',9f7.3)", self%f(:, 2)
484  print"('fBz=',9f7.3)", self%f(:, 3)
485 #endif
486 
487  kk = 0
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)
494  end do
495  end do
496 
497  end do
498  end do
499 
500  end subroutine sll_s_solve_maxwell_2d_diga
501 
502  function dof_local(edge, dof, degree)
503 
504  sll_int32 :: dof_local
505  sll_int32 :: edge, dof, degree
506 
507  select case (edge)
508  case (south)
509  dof_local = (dof - 1)*(degree + 1) + 1
510  case (east)
511  dof_local = degree*(degree + 1) + dof
512  case (north)
513  dof_local = dof*(degree + 1)
514  case (west)
515  dof_local = dof
516  end select
517 
518  end function dof_local
519 
520  function dof_neighbor(edge, dof, degree)
521 
522  sll_int32 :: dof_neighbor
523  sll_int32 :: edge, dof, degree
524 
525  select case (edge)
526  case (south)
527  dof_neighbor = dof*(degree + 1)
528  case (east)
529  dof_neighbor = dof
530  case (north)
531  dof_neighbor = (dof - 1)*(degree + 1) + 1
532  case (west)
533  dof_neighbor = degree*(degree + 1) + dof
534  !dof_neighbor = (degree+1)*(degree+1)-(dof-1)*(degree+1)
535  end select
536 
537  end function dof_neighbor
538 
540  subroutine compute_normals(tau, bc_south, bc_east, bc_north, bc_west, &
541  i, j, d, cell)
542 
543  sll_transformation, pointer :: tau
544  type(cell_type) :: cell
545  sll_int32 :: i, j, d
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
550  sll_real64 :: xk, wk
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)
555  sll_real64 :: length
556  sll_int32 :: side
557  sll_int32 :: bc_south
558  sll_int32 :: bc_east
559  sll_int32 :: bc_north
560  sll_int32 :: bc_west
561  sll_int32 :: k
562  sll_int32 :: error
563 
564  associate(lm => tau%mesh)
565 
566  cell%i = i
567  cell%j = j
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
570 
571  cell%eta1_max = cell%eta1_min + lm%delta_eta1
572  cell%eta2_max = cell%eta2_min + lm%delta_eta2
573 
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
578 
579  do side = 1, 4
580  sll_clear_allocate(cell%edge(side)%vec_norm(1:d + 1, 1:2), error)
581  cell%edge(side)%bc_type = sll_p_interior
582  end do
583 
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
588 
589  x = sll_f_gauss_lobatto_points(d + 1)
590  w = sll_f_gauss_lobatto_weights(d + 1)
591 
592  length = 0._f64
593  a = cell%eta1_min
594  b = cell%eta1_max
595  c1 = 0.5_f64*(b - a)
596  c2 = 0.5_f64*(b + a)
597  do k = 1, d + 1
598  xk = c1*x(k) + c2
599  wk = c1*w(k)
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
608  end do
609 
610  cell%edge(south)%length = length
611  cell%edge(south)%vec_norm = vec_norm/lm%delta_eta1
612 
613  length = 0._f64
614  a = cell%eta2_min
615  b = cell%eta2_max
616  c1 = 0.5_f64*(b - a)
617  c2 = 0.5_f64*(b + a)
618  do k = 1, d + 1
619  xk = c1*x(k) + c2
620  wk = c1*w(k)
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
629  end do
630 
631  cell%edge(east)%length = length
632  cell%edge(east)%vec_norm = vec_norm/lm%delta_eta2
633 
634  length = 0._f64
635  a = cell%eta1_min
636  b = cell%eta1_max
637  c1 = 0.5_f64*(b - a)
638  c2 = 0.5_f64*(b + a)
639  do k = 1, d + 1
640  xk = c1*x(k) + c2
641  wk = c1*w(k)
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
650  end do
651 
652  cell%edge(north)%length = length
653  cell%edge(north)%vec_norm = vec_norm/lm%delta_eta1
654 
655  length = 0._f64
656  a = cell%eta2_min
657  b = cell%eta2_max
658  c1 = 0.5_f64*(b - a)
659  c2 = 0.5_f64*(b + a)
660  do k = 1, d + 1
661  xk = c1*x(k) + c2
662  wk = c1*w(k)
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
671  end do
672 
673  cell%edge(west)%length = length
674  cell%edge(west)%vec_norm = vec_norm/lm%delta_eta2
675 
676  end associate
677 
678  end subroutine compute_normals
679 
680 end module sll_m_maxwell_2d_diga
Cartesian mesh basic types.
Abstract class for coordinate transformations.
Solve Maxwell equations on cartesian domain with Disconituous Galerkine method:
subroutine, public sll_s_dg_field_2d_init(this, degree, tau, init_function)
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.
    Report Typos and Errors