10 REAL(KIND=8),
PARAMETER,
PRIVATE ::
alpha=0.6d0
21 interface_h_mu, hn, bn, phin, hn1, bn1, phin1, vel, stab_in, sigma_in, &
22 r_fourier, index_fourier, mu_h_field, mu_phi, time, dt_in, rem, list_mode, &
23 h_phi_per, la_h, la_pmag, la_phi, la_mhd, one_over_sigma_ns_in, jj_v_to_h)
39 #include "petsc/finclude/petscsys.h" 40 #include "petsc/finclude/petscmat.h" 41 #include "petsc/finclude/petscksp.h" 42 #include "petsc/finclude/petscvec.h" 45 TYPE(
mesh_type),
INTENT(IN) :: H_mesh, phi_mesh, pmag_mesh
46 TYPE(
interface_type),
INTENT(IN) :: interface_H_phi, interface_H_mu
47 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
48 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: vel
49 REAL(KIND=8),
DIMENSION(H_mesh%np,6,SIZE(list_mode)),
INTENT(INOUT) :: Hn, Hn1
50 REAL(KIND=8),
DIMENSION(H_mesh%np,6,SIZE(list_mode)),
INTENT(INOUT) :: Bn, Bn1
51 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: phin, phin1
52 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab_in
53 REAL(KIND=8),
INTENT(IN) :: R_fourier
54 INTEGER,
INTENT(IN) :: index_fourier
55 REAL(KIND=8),
INTENT(IN) :: mu_phi, time, dt_in, Rem
56 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_in, mu_H_field
60 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_ns_in
61 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
63 REAL(KIND=8),
SAVE :: dt
64 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE,
SAVE :: sigma_ns_bar
65 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE,
SAVE :: sigma_np
66 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE,
SAVE :: sigma
67 REAL(KIND=8),
SAVE :: sigma_min
69 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: H_mode_global_js_D
70 TYPE(
dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: H_global_D
71 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: pmag_mode_global_js_D
72 TYPE(
dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: pmag_global_D
73 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: phi_mode_global_js_D
74 TYPE(
dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: phi_global_D
75 INTEGER,
DIMENSION(:),
POINTER,
SAVE :: pmag_js_D, phi_js_D
76 INTEGER,
DIMENSION(:),
ALLOCATABLE,
SAVE :: Dirichlet_bdy_H_sides
77 LOGICAL,
SAVE :: once=.true.
78 INTEGER,
SAVE :: m_max_c
80 REAL(KIND=8),
DIMENSION(3),
SAVE :: stab
81 INTEGER,
SAVE :: my_petscworld_rank
82 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: sigma_curl_gauss_bdy
83 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: J_over_sigma_gauss_bdy
84 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: sigma_curl_gauss_inter_mu
85 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: J_over_sigma_gauss_inter_mu
86 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: sigma_tot_gauss_Neumann
87 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:),
SAVE :: sigma_nj_m
88 REAL(KIND=8),
DIMENSION(H_mesh%gauss%l_G*H_mesh%me,6,SIZE(list_mode)) :: sigma_curl_gauss
89 REAL(KIND=8),
DIMENSION(H_mesh%gauss%l_G*H_mesh%me,6,SIZE(list_mode)) :: J_over_sigma_gauss
90 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),6,SIZE(Hn,3)) :: H_ns
91 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),2,SIZE(Hn,3)) :: one_over_sigma_tot
92 LOGICAL,
ALLOCATABLE,
DIMENSION(:) :: Dir_pmag
93 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: rhs_H
94 REAL(KIND=8),
DIMENSION(phi_mesh%np,2) :: rhs_phi
95 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),6,SIZE(Hn,3)) :: NL, H_ext, B_ext
96 REAL(KIND=8),
DIMENSION(3) :: temps_par
97 INTEGER,
POINTER,
DIMENSION(:) :: H_ifrom, pmag_ifrom, phi_ifrom, H_p_phi_ifrom
98 REAL(KIND=8),
DIMENSION(phi_mesh%np, 2) :: phin_p1
99 REAL(KIND=8),
DIMENSION(H_mesh%np, 6) :: Hn_p1
100 LOGICAL,
DIMENSION(H_mesh%mes) :: virgin1
101 LOGICAL,
DIMENSION(phi_mesh%mes) :: virgin2
102 INTEGER :: mode, k, i, n, m, ms, code, nj, j, count
103 INTEGER :: nb_procs, bloc_size, m_max_pad
104 REAL(KIND=8) :: tps, nr_vel, tps_tot, tps_cumul, norm
106 REAL(KIND=8) :: one_and_half
107 DATA one_and_half/1.5d0/
109 petscerrorcode :: ierr
110 mpi_comm,
DIMENSION(:),
POINTER :: comm_one_d
111 mat,
DIMENSION(:),
POINTER,
SAVE :: h_p_phi_mat1, h_p_phi_mat2
112 mat :: tampon1, tampon2, precond1, precond2
113 ksp,
DIMENSION(:),
POINTER,
SAVE :: h_p_phi_ksp1, h_p_phi_ksp2
114 vec,
SAVE :: vx_1, vb_1, vx_1_ghost, vx_2, vb_2, vx_2_ghost
127 CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
134 IF (
inputs%type_pb==
'fhd')
THEN 140 n =
SIZE(h_ifrom)+
SIZE(pmag_ifrom)+
SIZE(phi_ifrom)
141 ALLOCATE(h_p_phi_ifrom(n))
142 IF (
SIZE(h_ifrom)/=0)
THEN 143 h_p_phi_ifrom(1:
SIZE(h_ifrom)) = h_ifrom
145 IF (
SIZE(pmag_ifrom)/=0)
THEN 146 h_p_phi_ifrom(
SIZE(h_ifrom)+1:
SIZE(h_ifrom)+
SIZE(pmag_ifrom)) = pmag_ifrom
148 IF (
SIZE(phi_ifrom)/=0)
THEN 149 h_p_phi_ifrom(
SIZE(h_ifrom)+
SIZE(pmag_ifrom)+1:)=phi_ifrom
152 n = 3*h_mesh%dom_np + pmag_mesh%dom_np + phi_mesh%dom_np
153 CALL veccreateghost(comm_one_d(1), n, &
154 petsc_determine,
SIZE(h_p_phi_ifrom), h_p_phi_ifrom, vx_1, ierr)
155 CALL vecghostgetlocalform(vx_1, vx_1_ghost, ierr)
156 CALL vecduplicate(vx_1, vb_1, ierr)
157 CALL veccreateghost(comm_one_d(1), n, &
158 petsc_determine,
SIZE(h_p_phi_ifrom), h_p_phi_ifrom, vx_2, ierr)
159 CALL vecghostgetlocalform(vx_2, vx_2_ghost, ierr)
160 CALL vecduplicate(vx_2, vb_2, ierr)
164 ALLOCATE(sigma(
SIZE(sigma_in)))
165 sigma = sigma_in * rem
168 CALL mpi_allreduce(minval(sigma),sigma_min,1,mpi_double_precision, mpi_min,comm_one_d(1), ierr)
213 m_max_c =
SIZE(list_mode)
217 ALLOCATE(sigma_nj_m(h_mesh%gauss%n_w,h_mesh%me))
218 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 219 ALLOCATE(sigma_ns_bar(
SIZE(hn,1)))
220 sigma_ns_bar = sigma_bar_in_fourier_space(h_mesh)*rem
224 DO nj = 1, h_mesh%gauss%n_w
226 IF (jj_v_to_h(j) == -1)
THEN 227 sigma_nj_m(nj,m) = sigma(m)
229 sigma_nj_m(nj,m) = sigma_ns_bar(j)
235 sigma_nj_m(:,m) = sigma(m)
239 ALLOCATE(sigma_np(
SIZE(hn,1)))
242 DO nj = 1, h_mesh%gauss%n_w
243 sigma_np(h_mesh%jj(nj,m)) = sigma_nj_m(nj,m)
251 ALLOCATE (dir_pmag(maxval(pmag_mesh%sides)))
253 DO ms = 1,
SIZE(dir_pmag)
254 IF (minval(abs(
inputs%list_dirichlet_sides_H-ms)) == 0)
THEN 255 dir_pmag(ms) = .true.
257 IF (minval(abs(
inputs%list_inter_H_phi-ms)) == 0)
THEN 258 dir_pmag(ms) = .true.
262 CALL dirichlet_nodes(pmag_mesh%jjs, pmag_mesh%sides, dir_pmag, pmag_js_d)
267 ALLOCATE(pmag_global_d(m_max_c))
269 ALLOCATE(pmag_global_d(i)%DRL(
SIZE(pmag_mode_global_js_d(i)%DIL)))
270 pmag_global_d(i)%DRL = 0.d0
278 IF (interface_h_phi%mes/=0)
THEN 279 virgin1(interface_h_phi%mesh1) = .false.
280 virgin2(interface_h_phi%mesh2) = .false.
282 IF (interface_h_mu%mes/=0)
THEN 283 virgin1(interface_h_mu%mesh1) = .false.
284 virgin1(interface_h_mu%mesh2) = .false.
288 DO ms = 1, h_mesh%mes
289 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))).LT.1d-12*h_mesh%global_diameter) cycle
290 IF (.NOT.virgin1(ms)) cycle
291 IF(minval(abs(h_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))==0) cycle
293 IF (
inputs%my_periodic%nb_periodic_pairs /=0)
THEN 294 IF (minval(abs(
inputs%my_periodic%list_periodic-h_mesh%sides(ms))) == 0) cycle
301 DO ms = 1, h_mesh%mes
302 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))).LT.1d-12*h_mesh%global_diameter) cycle
303 IF (.NOT.virgin1(ms)) cycle
304 IF(minval(abs(h_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))==0) cycle
306 IF (
inputs%my_periodic%nb_periodic_pairs /=0)
THEN 307 IF (minval(abs(
inputs%my_periodic%list_periodic-h_mesh%sides(ms))) == 0) cycle
315 DO ms = 1, pmag_mesh%mes
316 IF (maxval(abs(pmag_mesh%rr(1,pmag_mesh%jjs(:,ms)))).LT.1d-12*pmag_mesh%global_diameter) cycle
317 IF(minval(abs(pmag_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))==0) cycle
318 IF(minval(abs(pmag_mesh%sides(ms)-
inputs%list_inter_H_phi))==0) cycle
320 IF (
inputs%my_periodic%nb_periodic_pairs /=0)
THEN 321 IF (minval(abs(
inputs%my_periodic%list_periodic-pmag_mesh%sides(ms))) == 0) cycle
328 DO ms = 1, pmag_mesh%mes
329 IF (maxval(abs(pmag_mesh%rr(1,pmag_mesh%jjs(:,ms)))).LT.1d-12*pmag_mesh%global_diameter) cycle
330 IF(minval(abs(pmag_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))==0) cycle
331 IF(minval(abs(pmag_mesh%sides(ms)-
inputs%list_inter_H_phi))==0) cycle
333 IF (
inputs%my_periodic%nb_periodic_pairs /=0)
THEN 334 IF (minval(abs(
inputs%my_periodic%list_periodic-pmag_mesh%sides(ms))) == 0) cycle
342 DO ms = 1, phi_mesh%mes
344 IF (phi_mesh%sides(ms)==index_fourier) cycle
346 IF (.NOT.virgin2(ms)) cycle
347 IF (maxval(abs(phi_mesh%rr(1,phi_mesh%jjs(:,ms)))).LT.1d-12*phi_mesh%global_diameter) cycle
348 IF (minval(abs(phi_mesh%sides(ms)-
inputs%phi_list_dirichlet_sides))==0) cycle
353 DO ms = 1, phi_mesh%mes
355 IF (phi_mesh%sides(ms)==index_fourier) cycle
357 IF (.NOT.virgin2(ms)) cycle
358 IF (maxval(abs(phi_mesh%rr(1,phi_mesh%jjs(:,ms)))).LT.1d-12*phi_mesh%global_diameter) cycle
359 IF (minval(abs(phi_mesh%sides(ms)-
inputs%phi_list_dirichlet_sides))==0) cycle
368 DO ms = 1, h_mesh%mes
369 IF (minval(abs(h_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))/=0) cycle
370 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))) .LT.1d-12*h_mesh%global_diameter) cycle
373 ALLOCATE(dirichlet_bdy_h_sides(n))
375 DO ms = 1, h_mesh%mes
376 IF (minval(abs(h_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))/=0) cycle
377 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))) .LT.1d-12*h_mesh%global_diameter) cycle
379 dirichlet_bdy_h_sides(n) = ms
383 ALLOCATE(h_global_d(m_max_c))
385 ALLOCATE(h_global_d(i)%DRL(
SIZE(h_mode_global_js_d(i)%DIL)))
394 ALLOCATE(phi_global_d(m_max_c))
396 ALLOCATE(phi_global_d(i)%DRL(
SIZE(phi_mode_global_js_d(i)%DIL)))
397 phi_global_d(i)%DRL = 0.d0
402 ALLOCATE(h_p_phi_mat1(m_max_c),h_p_phi_ksp1(m_max_c))
403 ALLOCATE(h_p_phi_mat2(m_max_c),h_p_phi_ksp2(m_max_c))
405 IF (
SIZE(dirichlet_bdy_h_sides).GE.1)
THEN 406 ALLOCATE(sigma_curl_gauss_bdy(h_mesh%gauss%l_Gs*
SIZE(dirichlet_bdy_h_sides),6,
SIZE(list_mode)))
407 ALLOCATE(j_over_sigma_gauss_bdy(h_mesh%gauss%l_Gs*
SIZE(dirichlet_bdy_h_sides),6,
SIZE(list_mode)))
409 ALLOCATE(sigma_curl_gauss_bdy(0,0,0))
410 ALLOCATE(j_over_sigma_gauss_bdy(0,0,0))
411 sigma_curl_gauss_bdy = 0.d0
412 j_over_sigma_gauss_bdy = 0.d0
415 IF (interface_h_mu%mes.GE.1)
THEN 416 ALLOCATE(sigma_curl_gauss_inter_mu(2*h_mesh%gauss%l_Gs*interface_h_mu%mes,6,
SIZE(list_mode)))
417 ALLOCATE(j_over_sigma_gauss_inter_mu(2*h_mesh%gauss%l_Gs*interface_h_mu%mes,6,
SIZE(list_mode)))
419 ALLOCATE(sigma_curl_gauss_inter_mu(0,0,0))
420 ALLOCATE(j_over_sigma_gauss_inter_mu(0,0,0))
421 sigma_curl_gauss_inter_mu = 0.d0
422 j_over_sigma_gauss_inter_mu = 0.d0
426 IF(my_petscworld_rank==0)
THEN 427 WRITE(*,*)
"WARNING, Neumann BC: either sigma or CURL(H)_Neumann is axisymmetric." 429 ALLOCATE(sigma_tot_gauss_neumann(h_mesh%gauss%l_Gs*
SIZE(
neumann_bdy_h_sides),2,
SIZE(list_mode)))
431 ALLOCATE(sigma_tot_gauss_neumann(0,0,0))
432 sigma_tot_gauss_neumann = 0.d0
454 mode,mu_h_field, mu_phi, one_and_half/dt, stab, r_fourier, index_fourier, &
455 la_h, la_pmag, la_phi, h_p_phi_mat1(i), h_p_phi_mat2(i), sigma_nj_m, sigma)
462 CALL mat_maxwell_mu(h_mesh, jj_v_to_h, interface_h_mu, mode, stab, &
463 mu_h_field, sigma, la_h, h_p_phi_mat1(i), h_p_phi_mat2(i), sigma_np)
470 mode, stab, la_h, h_p_phi_mat1(i), h_p_phi_mat2(i), sigma_np, sigma)
476 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 478 h_phi_per%perlist, h_p_phi_mat1(i), la_mhd)
480 h_phi_per%perlist, h_p_phi_mat2(i), la_mhd)
498 CALL init_solver(
inputs%my_par_H_p_phi,h_p_phi_ksp1(i),h_p_phi_mat1(i),comm_one_d(1),&
500 CALL init_solver(
inputs%my_par_H_p_phi,h_p_phi_ksp2(i),h_p_phi_mat2(i),comm_one_d(1),&
506 CALL matdestroy(h_p_phi_mat1(i),ierr)
507 CALL matdestroy(h_p_phi_mat2(i),ierr)
515 CALL mpi_comm_rank(petsc_comm_world, my_petscworld_rank, code)
519 nr_vel =
norm_sf(comm_one_d,
'L2', h_mesh, list_mode, vel)
522 IF (nr_vel .LE. 1.d-10)
THEN 524 ELSE IF (
inputs%type_pb==
"fhd")
THEN 527 CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
528 bloc_size =
SIZE(vel,1)/nb_procs+1
529 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
533 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 535 one_over_sigma_tot = 0.d0
537 DO nj = 1, h_mesh%gauss%n_w
540 IF (jj_v_to_h(j) /= -1)
THEN 541 h_ns(j,:,:) = 2*hn(j,:,:)- hn1(j,:,:)
542 one_over_sigma_tot(j,:,:) = one_over_sigma_ns_in(jj_v_to_h(j),:,:)/rem
544 DO i = 1,
SIZE(list_mode)
547 one_over_sigma_tot(j,1,i) = 1.d0/sigma(m)
556 one_over_sigma_tot, sigma_nj_m, sigma, sigma_curl_gauss)
557 IF (
SIZE(dirichlet_bdy_h_sides).GE.1)
THEN 559 one_over_sigma_tot, sigma_np, sigma, sigma_curl_gauss_bdy)
561 sigma_curl_gauss_bdy = 0.d0
563 IF (interface_h_mu%mes.GE.1)
THEN 565 one_over_sigma_tot, sigma_np, sigma, sigma_curl_gauss_inter_mu)
567 sigma_curl_gauss_inter_mu = 0.d0
572 mu_h_field, mu_phi, one_over_sigma_tot, time, sigma, j_over_sigma_gauss)
573 IF (
SIZE(dirichlet_bdy_h_sides).GE.1)
THEN 575 list_mode, b_ext, mu_h_field, mu_phi, one_over_sigma_tot, time, sigma,&
576 j_over_sigma_gauss_bdy)
578 j_over_sigma_gauss_bdy = 0.d0
580 IF (interface_h_mu%mes.GE.1)
THEN 582 list_mode, b_ext, mu_h_field, mu_phi, one_over_sigma_tot, time, sigma,&
583 j_over_sigma_gauss_inter_mu)
585 j_over_sigma_gauss_inter_mu=0.d0
591 list_mode, one_over_sigma_tot, sigma_tot_gauss_neumann)
593 sigma_tot_gauss_neumann = 0.d0
596 sigma_curl_gauss = 0.d0
597 sigma_curl_gauss_bdy = 0.d0
598 sigma_curl_gauss_inter_mu = 0.d0
599 j_over_sigma_gauss = 0.d0
600 j_over_sigma_gauss_bdy = 0.d0
601 j_over_sigma_gauss_inter_mu= 0.d0
602 sigma_tot_gauss_neumann = 0.d0
605 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
618 rhs_h(:,k) = (4*bn(:,k,i)-bn1(:,k,i))/(2*dt)
621 rhs_phi(:,k) = mu_phi*(4*phin(:,k,i)-phin1(:,k,i))/(2*dt)
625 rhs_h, nl(:,:,i), la_h, la_phi, vb_1, vb_2, b_ext(:,:,i),&
626 sigma_curl_gauss(:,:,i), j_over_sigma_gauss(:,:,i))
633 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
637 CALL courant_mu(h_mesh, interface_h_mu, sigma, mu_h_field, time, mode, nl(:,:,i), &
638 la_h, vb_1, vb_2, b_ext(:,:,i), j_over_sigma_gauss_inter_mu(:,:,i), &
639 sigma_curl_gauss_inter_mu(:,:,i))
640 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
646 sigma, mu_h_field, time, mode, nl(:,:,i), stab, la_h, vb_1,vb_2, b_ext(:,:,i), &
647 j_over_sigma_gauss_bdy(:,:,i), sigma_curl_gauss_bdy(:,:,i))
655 CALL surf_int(h_mesh, phi_mesh, pmag_mesh, interface_h_phi, interface_h_mu,
inputs%list_dirichlet_sides_H, &
656 sigma, mu_phi, mu_h_field, time, mode, la_h, la_phi, la_pmag, vb_1, vb_2, &
657 sigma_tot_gauss_neumann(:,:,i), r_fourier, index_fourier)
659 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
665 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 666 CALL periodic_rhs_petsc(h_phi_per%n_bord, h_phi_per%list, h_phi_per%perlist, vb_1, la_mhd)
667 CALL periodic_rhs_petsc(h_phi_per%n_bord, h_phi_per%list, h_phi_per%perlist, vb_2, la_mhd)
674 pmag_global_d(i)%DRL = 0.d0
675 CALL dirichlet_rhs(pmag_mode_global_js_d(i)%DIL-1,pmag_global_d(i)%DRL,vb_1)
676 CALL dirichlet_rhs(pmag_mode_global_js_d(i)%DIL-1,pmag_global_d(i)%DRL,vb_2)
678 IF (
SIZE(phi_js_d)>0)
THEN 685 phi_global_d(i)%DRL(1:n) = phiexact(1,phi_mesh%rr(1:2,phi_js_d), mode, mu_phi, time)
686 IF (
SIZE(phi_global_d(i)%DRL)>n) phi_global_d(i)%DRL(n+1:)=0.d0
687 CALL dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_1)
688 phi_global_d(i)%DRL(1:n) = phiexact(2,phi_mesh%rr(1:2,phi_js_d), mode, mu_phi, time)
689 IF (
SIZE(phi_global_d(i)%DRL)>n) phi_global_d(i)%DRL(n+1:)=0.d0
690 CALL dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_2)
696 phi_global_d(i)%DRL=0.d0
697 CALL dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_1)
698 CALL dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_2)
702 h_global_d(i)%DRL = 0.d0
703 CALL dirichlet_rhs(h_mode_global_js_d(i)%DIL-1,h_global_d(i)%DRL,vb_1)
704 CALL dirichlet_rhs(h_mode_global_js_d(i)%DIL-1,h_global_d(i)%DRL,vb_2)
706 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
711 IF (
inputs%my_par_H_p_phi%verbose .AND. (i==1))
WRITE(*,*)
'start solving' 714 CALL solver(h_p_phi_ksp1(i),vb_1,vx_1,reinit=.false.,
verbose=
inputs%my_par_H_p_phi%verbose)
716 CALL vecghostupdatebegin(vx_1,insert_values,scatter_forward,ierr)
717 CALL vecghostupdateend(vx_1,insert_values,scatter_forward,ierr)
718 IF (h_mesh%me/=0)
THEN 719 CALL extract(vx_1_ghost,1,1,la_mhd,hn_p1(:,1))
720 CALL extract(vx_1_ghost,2,2,la_mhd,hn_p1(:,4))
721 CALL extract(vx_1_ghost,3,3,la_mhd,hn_p1(:,5))
723 IF (phi_mesh%me/=0)
THEN 724 CALL extract(vx_1_ghost,5,5,la_mhd,phin_p1(:,1))
727 CALL solver(h_p_phi_ksp2(i),vb_2,vx_2,reinit=.false.,
verbose=
inputs%my_par_H_p_phi%verbose)
728 CALL vecghostupdatebegin(vx_2,insert_values,scatter_forward,ierr)
729 CALL vecghostupdateend(vx_2,insert_values,scatter_forward,ierr)
730 IF (h_mesh%me/=0)
THEN 731 CALL extract(vx_2_ghost,1,1,la_mhd,hn_p1(:,2))
732 CALL extract(vx_2_ghost,2,2,la_mhd,hn_p1(:,3))
733 CALL extract(vx_2_ghost,3,3,la_mhd,hn_p1(:,6))
735 IF (phi_mesh%me/=0)
THEN 736 CALL extract(vx_2_ghost,5,5,la_mhd,phin_p1(:,2))
739 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
747 IF (h_mesh%me /=0)
THEN 752 IF (phi_mesh%me /=0 )
THEN 759 IF (h_mesh%me /=0)
THEN 760 hn1(:,:,i) = hn(:,:,i)
762 hn(:,1,i) = hn_p1(:,1)
763 hn(:,4,i) = hn_p1(:,4)
764 hn(:,5,i) = hn_p1(:,5)
766 hn(:,2,i) = hn_p1(:,2)
767 hn(:,3,i) = hn_p1(:,3)
768 hn(:,6,i) = hn_p1(:,6)
771 bn1(:,k,i) = bn(:,k,i)
772 bn(:,k,i) = mu_h_field*hn(:,k,i)
777 IF (phi_mesh%me /= 0)
THEN 778 phin1(:,:,i) = phin(:,:,i)
780 phin(:,1,i) = phin_p1(:,1)
782 phin(:,2,i) = phin_p1(:,2)
784 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
791 IF (
inputs%verbose_divergence)
THEN 792 norm =
norm_sf(comm_one_d,
'L2', h_mesh, list_mode, bn)
804 mode, mu_h_field, mu_phi, c_mass, stab, r_fourier, index_fourier, &
805 la_h, la_pmag, la_phi, h_p_phi_mat1, h_p_phi_mat2, sigma_nj_m, sigma)
811 #include "petsc/finclude/petsc.h" 818 INTEGER,
INTENT(IN) :: mode
819 REAL(KIND=8),
INTENT(IN) :: mu_phi, c_mass
820 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
821 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
822 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_w,H_mesh%me),
INTENT(IN) :: sigma_nj_m
823 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN):: sigma
824 REAL(KIND=8),
OPTIONAL :: R_fourier
825 INTEGER,
OPTIONAL :: index_fourier
826 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_ws,phi_mesh%gauss%l_Gs) :: w_cs
827 REAL(KIND=8),
DIMENSION(2, H_mesh%gauss%n_w, phi_mesh%gauss%l_Gs, H_mesh%mes) :: dw_cs
828 INTEGER :: m, l, ms, ls, ni, nj, k, i, j, &
829 n_ws1, n_ws2, n_w2, n_w1, m1, m2, ki, kj,ib,jb, ms1, ms2
830 REAL(KIND=8) :: x, y, hm1, stab_div, stab_colle_H_phi
831 REAL(KIND=8) :: ray, error
832 LOGICAL :: mark=.false.
833 REAL(KIND=8),
DIMENSION(3,H_mesh%gauss%n_w,pmag_mesh%gauss%n_w) :: THpmag
834 REAL(KIND=8),
DIMENSION(pmag_mesh%gauss%n_w,pmag_mesh%gauss%n_w) :: Tpmag
835 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: TH
836 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_w,phi_mesh%gauss%n_w):: TPhi
864 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: Hsij
865 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_w,phi_mesh%gauss%n_w) :: Phisij
866 REAL(KIND=8),
DIMENSION(6,phi_mesh%gauss%n_w,phi_mesh%gauss%n_w) :: Sij
888 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: ww_H
890 REAL(KIND=8),
DIMENSION(:,:,:,:),
POINTER :: dw_H
892 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: rj_H
894 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: ww_phi
896 REAL(KIND=8),
DIMENSION(:,:,:,:),
POINTER :: dw_phi
899 REAL(KIND=8),
DIMENSION(2,pmag_mesh%gauss%n_w,H_mesh%gauss%l_G) :: dwp
900 REAL(KIND=8),
DIMENSION(pmag_mesh%gauss%n_w,H_mesh%gauss%l_G) :: wwp
902 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: rj_phi
904 REAL(KIND=8),
DIMENSION(2,phi_mesh%gauss%l_Gs) :: gauss1, gauss2
906 REAL(KIND=8) :: ref, diff, mu_H, c_mu_H, c_mu_phi, muhl, &
907 dzmuhl, drmuhl, c_div, hloc, viscolm, xij, eps
909 REAL(KIND=8) :: c_sym=.0d0
912 REAL(KIND=8) :: c_lap
917 REAL(KIND=8),
DIMENSION(3*H_mesh%gauss%n_w+pmag_mesh%gauss%n_w+ &
phi_mesh%gauss%n_w , 3*H_mesh%gauss%n_w+pmag_mesh%gauss%n_w+ &
phi_mesh%gauss%n_w) :: mat_loc1, mat_loc2
918 INTEGER ,
DIMENSION(3*H_mesh%gauss%n_w+pmag_mesh%gauss%n_w+ &
phi_mesh%gauss%n_w) :: idxn, jdxn
921 INTEGER :: n_wpmag, n_wH, n_wphi, ix, jx
923 REAL(KIND=8) :: sigma_np_gauss
925 mat :: h_p_phi_mat1, h_p_phi_mat2
926 petscerrorcode :: ierr
928 CALL matzeroentries (h_p_phi_mat1, ierr)
929 CALL matzeroentries (h_p_phi_mat2, ierr)
930 CALL matsetoption (h_p_phi_mat1, mat_row_oriented, petsc_false, ierr)
931 CALL matsetoption (h_p_phi_mat2, mat_row_oriented, petsc_false, ierr)
935 stab_colle_h_phi = stab(2)
939 c_mu_phi = c_mass*mu_phi
941 ww_h => h_mesh%gauss%ww
942 dw_h => h_mesh%gauss%dw
943 rj_h => h_mesh%gauss%rj
944 ww_phi => phi_mesh%gauss%ww
945 dw_phi => phi_mesh%gauss%dw
946 rj_phi => phi_mesh%gauss%rj
948 n_wh = h_mesh%gauss%n_w
949 n_wpmag = pmag_mesh%gauss%n_w
950 n_wphi = phi_mesh%gauss%n_w
960 DO l = 1, h_mesh%gauss%l_G
962 hloc = (sqrt(sum(h_mesh%gauss%rj(:,m)))/h_mesh%global_diameter)**(2*
alpha)
965 muhl = sum(mu_h_field(h_mesh%jj(:,m))*ww_h(:,l))
966 drmuhl = sum(mu_h_field(h_mesh%jj(:,m))*dw_h(1,:,l,m))
967 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m))*dw_h(2,:,l,m))
970 sigma_np_gauss = sum(sigma_nj_m(:,m)*ww_h(:,l))
980 c_div = stab_div*hloc/(
inputs%mu_min**2*
inputs%sigma_min)
983 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
984 ray = ray + h_mesh%rr(1,i)*ww_h(ni,l)
987 DO ni = 1, h_mesh%gauss%n_w
988 DO nj = 1, h_mesh%gauss%n_w
993 th(1,ni,nj) = th(1,ni,nj) + rj_h(l,m) * ray* ( &
996 c_mass*mu_h_field(j)*ww_h(nj,l)*ww_h(ni,l) &
997 + (dw_h(2,ni,l,m)*dw_h(2,nj,l,m) + mode**2/ray**2*ww_h(ni,l)*ww_h(nj,l))/sigma_np_gauss &
999 + c_div*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) + ww_h(ni,l)*drmuhl) &
1000 *(muhl*(ww_h(nj,l)/ray+dw_h(1,nj,l,m)) + ww_h(nj,l)*drmuhl))
1005 th(2,ni,nj) = th(2,ni,nj)+ rj_h(l,m) * ray* ( &
1006 mode/ray**2 * ww_h(ni,l)*(ww_h(nj,l)+ray*dw_h(1,nj,l,m))/sigma_np_gauss &
1008 + c_div*mode/ray*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) &
1009 + ww_h(ni,l)*drmuhl)*muhl*ww_h(nj,l))
1013 th(8,ni,nj) = th(8,ni,nj)+ rj_h(l,m) * ray* ( &
1014 - mode/ray**2 * ww_h(ni,l)*(ww_h(nj,l)+ray*dw_h(1,nj,l,m))/sigma_np_gauss &
1016 - c_div*mode/ray*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) &
1017 + ww_h(ni,l)*drmuhl)*muhl*ww_h(nj,l))
1021 th(3,ni,nj) = th(3,ni,nj)+ rj_h(l,m) * ray* ( &
1022 - dw_h(2,ni,l,m)*dw_h(1,nj,l,m)/sigma_np_gauss &
1024 + c_div*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) + ww_h(ni,l)*drmuhl)*&
1025 (muhl*dw_h(2,nj,l,m) + ww_h(nj,l)*dzmuhl))
1029 th(4,ni,nj) = th(4,ni,nj) + rj_h(l,m) * ray* ( &
1032 c_mass*mu_h_field(j)*ww_h(nj,l)*ww_h(ni,l) &
1033 + (dw_h(2,ni,l,m)*dw_h(2,nj,l,m) &
1034 + 1/ray**2 *(ww_h(ni,l)+ray*dw_h(1,ni,l,m))*(ww_h(nj,l)&
1035 +ray*dw_h(1,nj,l,m)))/sigma_np_gauss &
1037 +c_div*muhl**2*mode**2/ray**2*ww_h(ni,l)*ww_h(nj,l))
1041 th(5,ni,nj) = th(5,ni,nj) + rj_h(l,m) * ray* (&
1042 + mode/ray*dw_h(2,ni,l,m)*ww_h(nj,l)/sigma_np_gauss &
1044 +c_div*mode/ray*muhl*ww_h(ni,l)*(muhl*dw_h(2,nj,l,m) + ww_h(nj,l)*dzmuhl))
1048 th(9,ni,nj) = th(9,ni,nj) + rj_h(l,m) * ray* (&
1049 - mode/ray*dw_h(2,ni,l,m)*ww_h(nj,l)/sigma_np_gauss &
1051 - c_div*mode/ray*muhl*ww_h(ni,l)*(muhl*dw_h(2,nj,l,m) + ww_h(nj,l)*dzmuhl))
1055 th(6,ni,nj) = th(6,ni,nj) + rj_h(l,m) * ray* ( &
1058 c_mass*mu_h_field(j)*ww_h(nj,l)*ww_h(ni,l) &
1059 + (mode**2/ray**2*ww_h(ni,l)*ww_h(nj,l) + dw_h(1,ni,l,m)*dw_h(1,nj,l,m))/sigma_np_gauss &
1061 + c_div*(muhl*dw_h(2,ni,l,m) + ww_h(ni,l)*dzmuhl) &
1062 *(muhl*dw_h(2,nj,l,m) + ww_h(nj,l)*dzmuhl))
1143 i = h_mesh%jj(ni, m)
1144 ib = la_h%loc_to_glob(ki,i)
1149 j = h_mesh%jj(nj, m)
1150 jb = la_h%loc_to_glob(kj,j)
1154 IF ((ki == 1) .AND. (kj == 1))
THEN 1155 mat_loc1(ix,jx) = th(1,ni,nj)
1156 mat_loc2(ix,jx) = th(1,ni,nj)
1157 ELSEIF ((ki == 1) .AND. (kj == 2))
THEN 1158 mat_loc1(ix,jx) = th(2,ni,nj)
1159 mat_loc2(ix,jx) = th(8,ni,nj)
1160 ELSEIF ((ki == 2) .AND. (kj == 1))
THEN 1161 mat_loc1(ix,jx) = th(2,nj,ni)
1162 mat_loc2(ix,jx) = th(8,nj,ni)
1163 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 1164 mat_loc1(ix,jx) = th(3,ni,nj)
1165 mat_loc2(ix,jx) = th(3,ni,nj)
1166 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 1167 mat_loc1(ix,jx) = th(3,nj,ni)
1168 mat_loc2(ix,jx) = th(3,nj,ni)
1169 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 1170 mat_loc1(ix,jx) = th(4,ni,nj)
1171 mat_loc2(ix,jx) = th(4,ni,nj)
1172 ELSEIF ((ki == 2) .AND. (kj == 3))
THEN 1173 mat_loc1(ix,jx) = th(5,ni,nj)
1174 mat_loc2(ix,jx) = th(9,ni,nj)
1175 ELSEIF ((ki == 3) .AND. (kj == 2))
THEN 1176 mat_loc1(ix,jx) = th(5,nj,ni)
1177 mat_loc2(ix,jx) = th(9,nj,ni)
1178 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 1179 mat_loc1(ix,jx) = th(6,ni,nj)
1180 mat_loc2(ix,jx) = th(6,ni,nj)
1187 CALL matsetvalues(h_p_phi_mat1, 3*n_wh, idxn(1:3*n_wh), 3*n_wh, jdxn(1:3*n_wh), &
1188 mat_loc1(1:3*n_wh,1:3*n_wh), add_values, ierr)
1189 CALL matsetvalues(h_p_phi_mat2, 3*n_wh, idxn(1:3*n_wh), 3*n_wh, jdxn(1:3*n_wh), &
1190 mat_loc2(1:3*n_wh,1:3*n_wh), add_values, ierr)
1194 DO m = 1, pmag_mesh%me
1196 hloc = stab_div*(sqrt(sum(pmag_mesh%gauss%rj(:,m)))/pmag_mesh%global_diameter)**(2*(1-
alpha))
1198 DO l = 1, pmag_mesh%gauss%l_G
1207 DO ni = 1, pmag_mesh%gauss%n_w
1208 i = pmag_mesh%jj(ni,m)
1209 ray = ray + pmag_mesh%rr(1,i)*pmag_mesh%gauss%ww(ni,l)
1212 viscolm = (pmag_mesh%global_diameter)**2*
inputs%mu_min**2*
inputs%sigma_min*hloc*pmag_mesh%gauss%rj(l,m)
1213 DO nj = 1, pmag_mesh%gauss%n_w
1214 j = pmag_mesh%jj(nj, m)
1215 DO ni = 1, pmag_mesh%gauss%n_w
1216 i = pmag_mesh%jj(ni, m)
1220 xij = xij + pmag_mesh%gauss%dw(k,nj,l,m) * pmag_mesh%gauss%dw(k,ni,l,m)
1223 tpmag(ni,nj) = tpmag(ni,nj) + ray * viscolm* xij &
1224 + viscolm*mode**2*pmag_mesh%gauss%ww(ni,l)*pmag_mesh%gauss%ww(nj,l)/ray
1229 DO ni = 1, pmag_mesh%gauss%n_w
1230 i = pmag_mesh%jj(ni, m)
1231 ib = la_pmag%loc_to_glob(1,i)
1233 DO nj = 1, pmag_mesh%gauss%n_w
1234 j = pmag_mesh%jj(nj, m)
1235 jb = la_pmag%loc_to_glob(1,j)
1239 CALL matsetvalues(h_p_phi_mat1, n_wpmag, idxn(1:n_wpmag), n_wpmag, jdxn(1:n_wpmag), &
1240 tpmag(1:n_wpmag,1:n_wpmag), add_values, ierr)
1241 CALL matsetvalues(h_p_phi_mat2, n_wpmag, idxn(1:n_wpmag), n_wpmag, jdxn(1:n_wpmag), &
1242 tpmag(1:n_wpmag,1:n_wpmag), add_values, ierr)
1247 DO m = 1, pmag_mesh%me
1248 IF (h_mesh%gauss%n_w==3)
THEN 1249 dwp=h_mesh%gauss%dw(:,:,:,m)
1252 dwp(:,1,:) = h_mesh%gauss%dw(:,1,:,m) + 0.5d0*(h_mesh%gauss%dw(:,5,:,m)+h_mesh%gauss%dw(:,6,:,m))
1253 dwp(:,2,:) = h_mesh%gauss%dw(:,2,:,m) + 0.5d0*(h_mesh%gauss%dw(:,6,:,m)+h_mesh%gauss%dw(:,4,:,m))
1254 dwp(:,3,:) = h_mesh%gauss%dw(:,3,:,m) + 0.5d0*(h_mesh%gauss%dw(:,4,:,m)+h_mesh%gauss%dw(:,5,:,m))
1255 wwp(1,:) = h_mesh%gauss%ww(1,:) + 0.5d0*(h_mesh%gauss%ww(5,:)+h_mesh%gauss%ww(6,:))
1256 wwp(2,:) = h_mesh%gauss%ww(2,:) + 0.5d0*(h_mesh%gauss%ww(6,:)+h_mesh%gauss%ww(4,:))
1257 wwp(3,:) = h_mesh%gauss%ww(3,:) + 0.5d0*(h_mesh%gauss%ww(4,:)+h_mesh%gauss%ww(5,:))
1261 DO l = 1, h_mesh%gauss%l_G
1263 DO ni = 1, h_mesh%gauss%n_w
1265 ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
1267 muhl = stab_div*ray*h_mesh%gauss%rj(l,m)*sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
1268 DO nj = 1, pmag_mesh%gauss%n_w
1269 j = pmag_mesh%jj(nj, m)
1270 DO ni = 1, h_mesh%gauss%n_w
1271 i = h_mesh%jj(ni, m)
1272 thpmag(1,ni,nj) = thpmag(1,ni,nj) + muhl*dwp(1,nj,l)*h_mesh%gauss%ww(ni,l)
1273 thpmag(2,ni,nj) = thpmag(2,ni,nj) - muhl*mode*wwp(nj,l)*h_mesh%gauss%ww(ni,l)/ray
1274 thpmag(3,ni,nj) = thpmag(3,ni,nj) + muhl*dwp(2,nj,l)*h_mesh%gauss%ww(ni,l)
1284 i = h_mesh%jj(ni, m)
1291 ib = la_h%loc_to_glob(k,i)
1292 ix = (k-1)*n_wh + ni
1295 j = pmag_mesh%jj(nj, m)
1296 jb = la_pmag%loc_to_glob(1,j)
1299 mat_loc1(ix,jx) = thpmag(k,ni,nj)
1300 mat_loc2(ix,jx) = eps*thpmag(k,ni,nj)
1305 CALL matsetvalues(h_p_phi_mat1, 3*n_wh, idxn(1:3*n_wh), n_wpmag, jdxn(1:n_wpmag), &
1306 mat_loc1(1:3*n_wh,1:n_wpmag), add_values, ierr)
1307 CALL matsetvalues(h_p_phi_mat2, 3*n_wh, idxn(1:3*n_wh), n_wpmag, jdxn(1:n_wpmag), &
1308 mat_loc2(1:3*n_wh,1:n_wpmag), add_values, ierr)
1313 i = pmag_mesh%jj(ni, m)
1314 ib = la_pmag%loc_to_glob(1,i)
1324 j = h_mesh%jj(nj, m)
1325 jb = la_h%loc_to_glob(k,j)
1326 jx = (k-1)*n_wh + nj
1328 mat_loc1(ix,jx) = - thpmag(k,nj,ni)
1329 mat_loc2(ix,jx) = - eps*thpmag(k,nj,ni)
1333 CALL matsetvalues(h_p_phi_mat1, n_wpmag, idxn(1:n_wpmag), 3*n_wh, jdxn(1:3*n_wh), &
1334 mat_loc1(1:n_wpmag,1:3*n_wh), add_values, ierr)
1335 CALL matsetvalues(h_p_phi_mat2, n_wpmag, idxn(1:n_wpmag), 3*n_wh, jdxn(1:3*n_wh), &
1336 mat_loc2(1:n_wpmag,1:3*n_wh), add_values, ierr)
1341 DO m = 1,phi_mesh%me
1345 DO l = 1, phi_mesh%gauss%l_G
1349 DO ni = 1, phi_mesh%gauss%n_w; i = phi_mesh%jj(ni,m)
1350 ray = ray + phi_mesh%rr(1,i)*ww_phi(ni,l)
1353 DO ni = 1, phi_mesh%gauss%n_w
1354 DO nj = 1, phi_mesh%gauss%n_w
1362 tphi(ni,nj) = tphi(ni,nj) + rj_phi(l,m) * ray* (c_mass+c_lap)*mu_phi &
1363 *(dw_phi(1,ni,l,m)*dw_phi(1,nj,l,m)+dw_phi(2,ni,l,m)*dw_phi(2,nj,l,m) &
1364 +mode**2/ray**2*ww_phi(ni,l)*ww_phi(nj,l))
1374 DO ni = 1, phi_mesh%gauss%n_w
1375 i = phi_mesh%jj(ni, m)
1376 ib = la_phi%loc_to_glob(1,i)
1378 DO nj = 1, phi_mesh%gauss%n_w
1379 j = phi_mesh%jj(nj, m)
1380 jb = la_phi%loc_to_glob(1,j)
1384 CALL matsetvalues(h_p_phi_mat1, n_wphi, idxn(1:n_wphi), n_wphi, jdxn(1:n_wphi), &
1385 tphi(1:n_wphi,1:n_wphi), add_values, ierr)
1386 CALL matsetvalues(h_p_phi_mat2, n_wphi, idxn(1:n_wphi), n_wphi, jdxn(1:n_wphi), &
1387 tphi(1:n_wphi,1:n_wphi), add_values, ierr)
1395 CALL gauss(phi_mesh)
1396 n_ws1 = h_mesh%gauss%n_ws
1397 n_ws2 = phi_mesh%gauss%n_ws
1398 n_w1 = h_mesh%gauss%n_w
1399 n_w2 = phi_mesh%gauss%n_w
1401 IF (h_mesh%gauss%n_ws ==
n_ws)
THEN 1403 DO ms = 1, interface_h_phi%mes
1405 ms2 = interface_h_phi%mesh2(ms)
1406 m2 = phi_mesh%neighs(ms2)
1407 ms1 = interface_h_phi%mesh1(ms)
1408 m1 = h_mesh%neighs(ms1)
1410 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
1411 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - phi_mesh%rr(:,phi_mesh%jjs(1,ms2)))**2)
1412 IF (diff/ref .LT. 1.d-10)
THEN 1416 w_cs(1,ls)=
wws(2,ls)
1417 w_cs(2,ls)=
wws(1,ls)
1418 w_cs(3,ls)=
wws(3,ls)
1419 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 1424 gauss2(1,ls) = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))*phi_mesh%gauss%wws(:,ls))
1425 gauss2(2,ls) = sum(phi_mesh%rr(2,phi_mesh%jjs(:,ms2))*phi_mesh%gauss%wws(:,ls))
1426 gauss1(1,ls) = sum( h_mesh%rr(1, h_mesh%jjs(:,ms1))* h_mesh%gauss%wws(:,ls))
1427 gauss1(2,ls) = sum( h_mesh%rr(2, h_mesh%jjs(:,ms1))* h_mesh%gauss%wws(:,ls))
1431 ref = sqrt(1.d-8+sum(gauss2(:,ls2)**2))
1434 diff = sqrt(sum((gauss2(:,ls2)-gauss1(:,ls1))**2))
1435 IF (diff .LT. 1.d-10)
THEN 1436 dw_cs(:,:,ls2,ms1) = h_mesh%gauss%dw_s(:,:,ls1,ms1)
1441 IF (.NOT.mark)
WRITE(*,*)
' BUG ' 1447 DO ms = 1, interface_h_phi%mes
1449 ms2 = interface_h_phi%mesh2(ms)
1450 m2 = phi_mesh%neighs(ms2)
1451 ms1 = interface_h_phi%mesh1(ms)
1452 m1 = h_mesh%neighs(ms1)
1454 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
1455 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - phi_mesh%rr(:,phi_mesh%jjs(1,ms2)))**2)
1456 IF (diff/ref .LT. 1.d-10)
THEN 1458 w_cs(1,ls)=
wws(1,ls)+0.5*
wws(3,ls)
1459 w_cs(2,ls)=
wws(2,ls)+0.5*
wws(3,ls)
1464 w_cs(1,ls)=
wws(2,ls)+0.5*
wws(3,ls)
1465 w_cs(2,ls)=
wws(1,ls)+0.5*
wws(3,ls)
1467 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 1472 dw_cs(1,:,ls,ms1) = h_mesh%gauss%dw(1,:,1,m1)
1473 dw_cs(2,:,ls,ms1) = h_mesh%gauss%dw(2,:,1,m1)
1480 DO ms = 1, interface_h_phi%mes
1482 ms2 = interface_h_phi%mesh2(ms)
1483 ms1 = interface_h_phi%mesh1(ms)
1484 m2 = phi_mesh%neighs(ms2)
1485 m1 = h_mesh%neighs(ms1)
1486 mu_h = sum(mu_h_field(h_mesh%jj(:,m1)))/h_mesh%gauss%n_w
1491 hm1 = stab_colle_h_phi/(sum(
rjs(:,ms2))*
inputs%sigma_min)
1503 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1504 x = hm1*
rjs(ls,ms2)*ray
1508 y = x * w_cs(ni,ls)*w_cs(nj,ls)
1509 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(
rnorms(2,ls,ms2)**2)
1510 hsij(4,ni,nj) = hsij(4,ni,nj) - y*
rnorms(1,ls,ms2)*
rnorms(2,ls,ms2)
1511 hsij(5,ni,nj) = hsij(5,ni,nj) + y
1512 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(
rnorms(1,ls,ms2)**2)
1527 i = interface_h_phi%jjs1(ni,ms)
1528 ib = la_h%loc_to_glob(ki,i)
1529 ix = (ki-1)*n_ws1+ni
1533 j = interface_h_phi%jjs1(nj,ms)
1534 jb = la_h%loc_to_glob(kj,j)
1535 jx = (kj-1)*n_ws1+nj
1537 IF ((ki == 1) .AND. (kj == 1))
THEN 1538 mat_loc1(ix,jx) = hsij(1,ni,nj)
1539 mat_loc2(ix,jx) = hsij(1,ni,nj)
1540 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 1541 mat_loc1(ix,jx) = hsij(4,ni,nj)
1542 mat_loc2(ix,jx) = hsij(4,ni,nj)
1543 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 1544 mat_loc1(ix,jx) = hsij(4,nj,ni)
1545 mat_loc2(ix,jx) = hsij(4,nj,ni)
1546 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 1547 mat_loc1(ix,jx) = hsij(5,ni,nj)
1548 mat_loc2(ix,jx) = hsij(5,ni,nj)
1549 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 1550 mat_loc1(ix,jx) = hsij(6,ni,nj)
1551 mat_loc2(ix,jx) = hsij(6,ni,nj)
1557 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1558 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1559 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1560 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1567 DO ls = 1, phi_mesh%gauss%l_Gs
1570 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1571 x =
rjs(ls,ms2)*ray/sigma(m1)
1578 y = x*w_cs(ni,ls)*w_cs(nj,ls)
1579 hsij(2,ni,nj) = hsij(2,ni,nj) + y * (-mode/ray)*(-
rnorms(1,ls,ms2))
1580 hsij(3,ni,nj) = hsij(3,ni,nj) + y * mode/ray *(-
rnorms(1,ls,ms2))
1581 hsij(5,ni,nj) = hsij(5,ni,nj) + y * (-1/ray) *(-
rnorms(1,ls,ms2))
1582 hsij(8,ni,nj) = hsij(8,ni,nj) + y * (-mode/ray)*(-
rnorms(2,ls,ms2))
1583 hsij(9,ni,nj) = hsij(9,ni,nj) + y * mode/ray *(-
rnorms(2,ls,ms2))
1597 i = interface_h_phi%jjs1(ni,ms)
1598 ib = la_h%loc_to_glob(ki,i)
1599 ix = (ki-1)*n_ws1 + ni
1603 j = interface_h_phi%jjs1(nj,ms)
1604 jb = la_h%loc_to_glob(kj,j)
1605 jx = (kj-1)*n_ws1 + nj
1607 IF ( (ki == 2) .AND. (kj == 1))
THEN 1608 mat_loc1(ix,jx) = hsij(2,ni,nj)
1609 mat_loc2(ix,jx) = hsij(3,ni,nj)
1610 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 1611 mat_loc1(ix,jx) = hsij(5,ni,nj)
1612 mat_loc2(ix,jx) = hsij(5,ni,nj)
1613 ELSEIF ( (ki == 2) .AND. (kj == 3))
THEN 1614 mat_loc1(ix,jx) = hsij(8,ni,nj)
1615 mat_loc2(ix,jx) = hsij(9,ni,nj)
1621 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1622 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1623 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1624 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1632 i = interface_h_phi%jjs1(ni,ms)
1633 ib = la_h%loc_to_glob(ki,i)
1634 ix = (ki-1)*n_ws1 + ni
1638 j = interface_h_phi%jjs1(nj,ms)
1639 jb = la_h%loc_to_glob(kj,j)
1640 jx = (kj-1)*n_ws1 + nj
1642 IF ( (kj == 2) .AND. (ki == 1))
THEN 1643 mat_loc1(ix,jx) = hsij(2,nj,ni)
1644 mat_loc2(ix,jx) = hsij(3,nj,ni)
1645 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN 1646 mat_loc1(ix,jx) = hsij(5,nj,ni)
1647 mat_loc2(ix,jx) = hsij(5,nj,ni)
1648 ELSEIF ( (kj == 2) .AND. (ki == 3))
THEN 1649 mat_loc1(ix,jx) = hsij(8,nj,ni)
1650 mat_loc2(ix,jx) = hsij(9,nj,ni)
1657 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1658 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1659 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1660 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1663 DO ls = 1, phi_mesh%gauss%l_Gs
1666 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1667 x =
rjs(ls,ms2)*ray /sigma(m1)
1675 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(-dw_cs(2,nj,ls,ms1))*(-
rnorms(2,ls,ms2))
1676 hsij(4,ni,nj) = hsij(4,ni,nj) + y* dw_cs(1,nj,ls,ms1) *(-
rnorms(2,ls,ms2))
1677 hsij(5,ni,nj) = hsij(5,ni,nj) + &
1678 y*(-dw_cs(2,nj,ls,ms1)*(-
rnorms(2,ls,ms2))-dw_cs(1,nj,ls,ms1)*(-
rnorms(1,ls,ms2)))
1679 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(-dw_cs(1,nj,ls,ms1))*(-
rnorms(1,ls,ms2))
1680 hsij(7,ni,nj) = hsij(7,ni,nj) + y* dw_cs(2,nj,ls,ms1) *(-
rnorms(1,ls,ms2))
1693 i = interface_h_phi%jjs1(ni,ms)
1694 ib = la_h%loc_to_glob(ki,i)
1695 ix = (ki-1)*n_ws1 + ni
1699 j = h_mesh%jj(nj,m1)
1700 jb = la_h%loc_to_glob(kj,j)
1701 jx = (kj-1)*n_w1 + nj
1703 IF ((ki == 1) .AND. (kj == 1))
THEN 1704 mat_loc1(ix,jx) = hsij(1,ni,nj)
1705 mat_loc2(ix,jx) = hsij(1,ni,nj)
1706 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 1707 mat_loc1(ix,jx) = hsij(4,ni,nj)
1708 mat_loc2(ix,jx) = hsij(4,ni,nj)
1709 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 1710 mat_loc1(ix,jx) = hsij(5,ni,nj)
1711 mat_loc2(ix,jx) = hsij(5,ni,nj)
1712 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 1713 mat_loc1(ix,jx) = hsij(6,ni,nj)
1714 mat_loc2(ix,jx) = hsij(6,ni,nj)
1715 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 1716 mat_loc1(ix,jx) = hsij(7,ni,nj)
1717 mat_loc2(ix,jx) = hsij(7,ni,nj)
1724 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
1725 mat_loc1(1:3*n_ws1,1:3*n_w1), add_values, ierr)
1726 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
1727 mat_loc2(1:3*n_ws1,1:3*n_w1), add_values, ierr)
1735 i = h_mesh%jj(ni,m1)
1736 ib = la_h%loc_to_glob(ki,i)
1737 ix = (ki-1)*n_w1 + ni
1741 j = interface_h_phi%jjs1(nj,ms)
1742 jb = la_h%loc_to_glob(kj,j)
1743 jx = (kj-1)*n_ws1 + nj
1745 IF ((kj == 1) .AND. (ki == 1))
THEN 1746 mat_loc1(ix,jx) = hsij(1,nj,ni)
1747 mat_loc2(ix,jx) = hsij(1,nj,ni)
1748 ELSEIF ((kj == 1) .AND. (ki == 3))
THEN 1749 mat_loc1(ix,jx) = hsij(4,nj,ni)
1750 mat_loc2(ix,jx) = hsij(4,nj,ni)
1751 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN 1752 mat_loc1(ix,jx) = hsij(5,nj,ni)
1753 mat_loc2(ix,jx) = hsij(5,nj,ni)
1754 ELSEIF ((kj == 3) .AND. (ki == 3))
THEN 1755 mat_loc1(ix,jx) = hsij(6,nj,ni)
1756 mat_loc2(ix,jx) = hsij(6,nj,ni)
1757 ELSEIF ((kj == 3) .AND. (ki == 1))
THEN 1758 mat_loc1(ix,jx) = hsij(7,nj,ni)
1759 mat_loc2(ix,jx) = hsij(7,nj,ni)
1765 CALL matsetvalues(h_p_phi_mat1, 3*n_w1, idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
1766 mat_loc1(1:3*n_w1,1:3*n_ws1), add_values, ierr)
1767 CALL matsetvalues(h_p_phi_mat2, 3*n_w1, idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
1768 mat_loc2(1:3*n_w1,1:3*n_ws1), add_values, ierr)
1781 DO ls = 1, phi_mesh%gauss%l_Gs
1784 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1785 x = hm1*
rjs(ls,ms2)*ray
1790 phisij(ni,nj) = phisij(ni,nj) + x*mode**2/ray**2*
wws(ni,ls)*
wws(nj,ls)
1801 i = interface_h_phi%jjs2(ni,ms)
1802 ib = la_phi%loc_to_glob(1,i)
1805 j = interface_h_phi%jjs2(nj,ms)
1806 jb = la_phi%loc_to_glob(1,j)
1810 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
1811 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
1812 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
1813 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
1819 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1820 x = hm1*
rjs(ls,ms2)*ray
1825 phisij(ni,nj) = phisij(ni,nj) + x*( &
1826 (
dw_s(2,ni,ls,ms2)*
rnorms(1,ls,ms2) -
dw_s(1,ni,ls,ms2)*
rnorms(2,ls,ms2))* &
1827 (
dw_s(2,nj,ls,ms2)*
rnorms(1,ls,ms2) -
dw_s(1,nj,ls,ms2)*
rnorms(2,ls,ms2)))
1837 i = phi_mesh%jj(ni, m2)
1838 ib = la_phi%loc_to_glob(1,i)
1841 j = phi_mesh%jj(nj, m2)
1842 jb = la_phi%loc_to_glob(1,j)
1846 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), n_w2, jdxn(1:n_w2), &
1847 phisij(1:n_w2,1:n_w2), add_values, ierr)
1848 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), n_w2, jdxn(1:n_w2), &
1849 phisij(1:n_w2,1:n_w2), add_values, ierr)
1863 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1864 x = hm1*
rjs(ls,ms2)*ray
1869 sij(3,ni,nj) = sij(3,ni,nj) + x*(mode/ray)*w_cs(ni,ls)*
wws(nj,ls)
1873 sij(4,:,:) = -sij(3,:,:)
1882 i = interface_h_phi%jjs1(ni,ms)
1883 ib = la_h%loc_to_glob(ki,i)
1886 j = interface_h_phi%jjs2(nj,ms)
1887 jb = la_phi%loc_to_glob(1,j)
1891 CALL matsetvalues(h_p_phi_mat1, n_ws1, idxn(1:n_ws1), n_ws2, jdxn(1:n_ws2), &
1892 sij(3,1:n_ws1,1:n_ws2), add_values, ierr)
1893 CALL matsetvalues(h_p_phi_mat2, n_ws1, idxn(1:n_ws1), n_ws2, jdxn(1:n_ws2), &
1894 sij(4,1:n_ws1,1:n_ws2), add_values, ierr)
1903 i = interface_h_phi%jjs2(ni,ms)
1904 ib = la_phi%loc_to_glob(1,i)
1907 j = interface_h_phi%jjs1(nj,ms)
1908 jb = la_h%loc_to_glob(kj,j)
1910 mat_loc1(ni,nj) = sij(3,nj,ni)
1911 mat_loc2(ni,nj) = sij(4,nj,ni)
1914 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws1, jdxn(1:n_ws1), &
1915 mat_loc1(1:n_ws2,1:n_ws1), add_values, ierr)
1916 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws1, jdxn(1:n_ws1), &
1917 mat_loc2(1:n_ws2,1:n_ws1), add_values, ierr)
1925 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1926 x = hm1*
rjs(ls,ms2)*ray
1932 sij(1,ni,nj) = sij(1,ni,nj) + &
1933 y*(-
dw_s(1,nj,ls,ms2)*
rnorms(2,ls,ms2)**2 +
dw_s(2,nj,ls,ms2)*
rnorms(1,ls,ms2)*
rnorms(2,ls,ms2))
1934 sij(5,ni,nj) = sij(5,ni,nj) + &
1935 y*(-
dw_s(2,nj,ls,ms2)*
rnorms(1,ls,ms2)**2 +
dw_s(1,nj,ls,ms2)*
rnorms(1,ls,ms2)*
rnorms(2,ls,ms2))
1948 i = interface_h_phi%jjs1(ni,ms)
1949 ib = la_h%loc_to_glob(ki,i)
1950 ix = (ki-1)*n_ws1 + ni
1953 j = phi_mesh%jj(nj,m2)
1954 jb = la_phi%loc_to_glob(1,j)
1958 mat_loc1(ix,jx) = sij(1,ni,nj)
1959 mat_loc2(ix,jx) = sij(1,ni,nj)
1960 ELSEIF (ki == 3)
THEN 1961 mat_loc1(ix,jx) = sij(5,ni,nj)
1962 mat_loc2(ix,jx) = sij(5,ni,nj)
1967 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), n_w2, jdxn(1:n_w2), &
1968 mat_loc1(1:3*n_ws1,1:n_w2), add_values, ierr)
1969 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), n_w2, jdxn(1:n_w2), &
1970 mat_loc2(1:3*n_ws1,1:n_w2), add_values, ierr)
1978 i = phi_mesh%jj(ni,m2)
1979 ib = la_phi%loc_to_glob(1,i)
1984 j = interface_h_phi%jjs1(nj,ms)
1985 jb = la_h%loc_to_glob(kj,j)
1986 jx = (kj-1)*n_ws1 + nj
1989 mat_loc1(ix,jx) = sij(1,nj,ni)
1990 mat_loc2(ix,jx) = sij(1,nj,ni)
1991 ELSEIF (kj == 3)
THEN 1992 mat_loc1(ix,jx) = sij(5,nj,ni)
1993 mat_loc2(ix,jx) = sij(5,nj,ni)
1998 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
1999 mat_loc1(1:n_w2,1:3*n_ws1), add_values, ierr)
2000 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
2001 mat_loc2(1:n_w2,1:3*n_ws1), add_values, ierr)
2016 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2017 x =
rjs(ls,ms2)*ray/sigma(m1)
2024 y = x *
wws(ni,ls)*w_cs(nj,ls)
2025 sij(1,ni,nj) = sij(1,ni,nj) + y*( mode/ray)**2*
rnorms(1,ls,ms2)
2026 sij(3,ni,nj) = sij(3,ni,nj) + y*( mode/ray**2)*
rnorms(1,ls,ms2)
2027 sij(4,ni,nj) = sij(4,ni,nj) + y*(-mode/ray**2)*
rnorms(1,ls,ms2)
2028 sij(5,ni,nj) = sij(5,ni,nj) + y*( mode/ray)**2*
rnorms(2,ls,ms2)
2040 i = interface_h_phi%jjs2(ni,ms)
2041 ib = la_phi%loc_to_glob(1,i)
2046 j = interface_h_phi%jjs1(nj,ms)
2047 jb = la_h%loc_to_glob(kj,j)
2048 jx = (kj-1)*n_ws1 + nj
2051 mat_loc1(ix,jx) = sij(1,ni,nj)
2052 mat_loc2(ix,jx) = sij(1,ni,nj)
2053 ELSEIF (kj == 2)
THEN 2054 mat_loc1(ix,jx) = sij(3,ni,nj)
2055 mat_loc2(ix,jx) = sij(4,ni,nj)
2056 ELSEIF (kj == 3)
THEN 2057 mat_loc1(ix,jx) = sij(5,ni,nj)
2058 mat_loc2(ix,jx) = sij(5,ni,nj)
2063 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), 3*n_ws1, jdxn(1:3*n_ws1), &
2064 mat_loc1(1:n_ws2,1:3*n_ws1), add_values, ierr)
2065 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), 3*n_ws1, jdxn(1:3*n_ws1), &
2066 mat_loc2(1:n_ws2,1:3*n_ws1), add_values, ierr)
2074 i = interface_h_phi%jjs1(ni,ms)
2075 ib = la_h%loc_to_glob(ki,i)
2076 ix = (ki-1)*n_ws1 + ni
2079 j = interface_h_phi%jjs2(nj,ms)
2080 jb = la_phi%loc_to_glob(1,j)
2084 mat_loc1(ix,jx) = sij(1,nj,ni)
2085 mat_loc2(ix,jx) = sij(1,nj,ni)
2086 ELSEIF (ki == 2)
THEN 2087 mat_loc1(ix,jx) = sij(3,nj,ni)
2088 mat_loc2(ix,jx) = sij(4,nj,ni)
2089 ELSEIF (ki == 3)
THEN 2090 mat_loc1(ix,jx) = sij(5,nj,ni)
2091 mat_loc2(ix,jx) = sij(5,nj,ni)
2096 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), n_ws2, jdxn(1:n_ws2), &
2097 mat_loc1(1:3*n_ws1,1:n_ws2), add_values, ierr)
2098 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), n_ws2, jdxn(1:n_ws2), &
2099 mat_loc2(1:3*n_ws1,1:n_ws2), add_values, ierr)
2107 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2108 x =
rjs(ls,ms2)*ray/sigma(m1)
2114 y = x*
wws(ni,ls)*mode/ray
2116 sij(3,ni,nj) = sij(3,ni,nj) + &
2117 y*(dw_cs(2,nj,ls,ms1)*
rnorms(2,ls,ms2) + dw_cs(1,nj,ls,ms1)*
rnorms(1,ls,ms2))
2121 sij(4,:,:) = -sij(3,:,:)
2127 i = interface_h_phi%jjs2(ni,ms)
2128 ib = la_phi%loc_to_glob(1,i)
2131 j = h_mesh%jj(nj,m1)
2132 jb = la_h%loc_to_glob(kj,j)
2136 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_w1, jdxn(1:n_w1), &
2137 sij(3,1:n_ws2,1:n_w1), add_values, ierr)
2138 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_w1, jdxn(1:n_w1), &
2139 sij(4,1:n_ws2,1:n_w1), add_values, ierr)
2148 i = h_mesh%jj(ni,m1)
2149 ib = la_h%loc_to_glob(ki,i)
2152 j = interface_h_phi%jjs2(nj,ms)
2153 jb = la_phi%loc_to_glob(1,j)
2155 mat_loc1(ix,jx) = sij(3,nj,ni)
2156 mat_loc2(ix,jx) = sij(4,nj,ni)
2159 CALL matsetvalues(h_p_phi_mat1, n_w1, idxn(1:n_w1), n_ws2, jdxn(1:n_ws2), &
2160 mat_loc1(1:n_w1,1:n_ws2), add_values, ierr)
2161 CALL matsetvalues(h_p_phi_mat2, n_w1, idxn(1:n_w1), n_ws2, jdxn(1:n_ws2), &
2162 mat_loc2(1:n_w1,1:n_ws2), add_values, ierr)
2169 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2170 x =
rjs(ls,ms2)*ray/sigma(m1)
2176 y = x*(
dw_s(2,ni,ls,ms2)*
rnorms(1,ls,ms2) -
dw_s(1,ni,ls,ms2)*
rnorms(2,ls,ms2))
2178 sij(1,ni,nj) = sij(1,ni,nj) + y *dw_cs(2,nj,ls,ms1)
2179 sij(5,ni,nj) = sij(5,ni,nj) + (-y)*dw_cs(1,nj,ls,ms1)
2191 i = phi_mesh%jj(ni,m2)
2192 ib = la_phi%loc_to_glob(1,i)
2196 j = h_mesh%jj(nj,m1)
2198 jb = la_h%loc_to_glob(kj,j)
2199 jx = (kj-1)*n_w1 + nj
2202 mat_loc1(ix,jx) = sij(1,ni,nj)
2203 mat_loc2(ix,jx) = sij(1,ni,nj)
2204 ELSEIF (kj == 3)
THEN 2205 mat_loc1(ix,jx) = sij(5,ni,nj)
2206 mat_loc2(ix,jx) = sij(5,ni,nj)
2211 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 3*n_w1, jdxn(1:3*n_w1), &
2212 mat_loc1(1:n_w2,1:3*n_w1), add_values, ierr)
2213 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 3*n_w1, jdxn(1:3*n_w1), &
2214 mat_loc2(1:n_w2,1:3*n_w1), add_values, ierr)
2222 i = h_mesh%jj(ni,m1)
2223 ib = la_h%loc_to_glob(ki,i)
2224 ix = (ki-1)*n_w1 + ni
2227 j = phi_mesh%jj(nj,m2)
2228 jb = la_phi%loc_to_glob(1,j)
2232 mat_loc1(ix,jx) = sij(1,nj,ni)
2233 mat_loc2(ix,jx) = sij(1,nj,ni)
2234 ELSEIF (ki == 3)
THEN 2235 mat_loc1(ix,jx) = sij(5,nj,ni)
2236 mat_loc2(ix,jx) = sij(5,nj,ni)
2241 CALL matsetvalues(h_p_phi_mat1, 3*n_w1, idxn(1:3*n_w1), n_w2, jdxn(1:n_w2), &
2242 mat_loc1(1:3*n_w1,1:n_w2), add_values, ierr)
2243 CALL matsetvalues(h_p_phi_mat2, 3*n_w1, idxn(1:3*n_w1), n_w2, jdxn(1:n_w2), &
2244 mat_loc2(1:3*n_w1,1:n_w2), add_values, ierr)
2251 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2252 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2253 x = c_lap*muhl*
rjs(ls,ms2)*ray
2256 sij(1,ni,nj) = sij(1,ni,nj) - x*w_cs(nj,ls)*
wws(ni,ls)*
rnorms(1,ls,ms2)
2257 sij(5,ni,nj) = sij(5,ni,nj) - x*w_cs(nj,ls)*
wws(ni,ls)*
rnorms(2,ls,ms2)
2265 i = interface_h_phi%jjs2(ni,ms)
2266 ib = la_phi%loc_to_glob(1,i)
2270 j = interface_h_phi%jjs1(nj,ms)
2271 jb = la_h%loc_to_glob(1,j)
2274 mat_loc1(ix,jx) = sij(1,ni,nj)
2275 mat_loc2(ix,jx) = sij(1,ni,nj)
2277 jb = la_h%loc_to_glob(3,j)
2280 mat_loc1(ix,jx) = sij(5,ni,nj)
2281 mat_loc2(ix,jx) = sij(5,ni,nj)
2284 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), 2*n_ws1, jdxn(1:2*n_ws1), &
2285 mat_loc1(1:n_ws2,1:2*n_ws1), add_values, ierr)
2286 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), 2*n_ws1, jdxn(1:2*n_ws1), &
2287 mat_loc2(1:n_ws2,1:2*n_ws1), add_values, ierr)
2298 IF (stab(2) > 1.d-12)
THEN 2307 ms2 = interface_h_phi%mesh2(ms)
2309 hm1 =(sum(
rjs(:,ms2))/h_mesh%global_diameter)**(2*
alpha-1)/(
inputs%sigma_min*
inputs%mu_min**2*h_mesh%global_diameter)
2314 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2317 DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,ms2)
2318 ray = ray + phi_mesh%rr(1,i)* phi_mesh%gauss%wws(ni,ls)
2324 ray = stab_div*ray*hm1*
rjs(ls,ms2)
2330 x = muhl**2*w_cs(ni,ls)*w_cs(nj,ls)*ray
2331 hsij(1,ni,nj) = hsij(1,ni,nj) + x*
rnorms(1,ls,ms2)**2
2332 hsij(4,ni,nj) = hsij(4,ni,nj) + x*
rnorms(1,ls,ms2)*
rnorms(2,ls,ms2)
2333 hsij(6,ni,nj) = hsij(6,ni,nj) + x*
rnorms(2,ls,ms2)**2
2337 x = muhl*mu_phi*w_cs(ni,ls)*(
dw_s(1,nj,ls,ms2)*
rnorms(1,ls,ms2) +&
2339 sij(1,ni,nj) = sij(1,ni,nj) - x*
rnorms(1,ls,ms2)
2340 sij(5,ni,nj) = sij(5,ni,nj) - x*
rnorms(2,ls,ms2)
2346 x = mu_phi**2*(
dw_s(1,ni,ls,ms2)*
rnorms(1,ls,ms2) +
dw_s(2,ni,ls,ms2)*
rnorms(2,ls,ms2))* &
2347 (
dw_s(1,nj,ls,ms2)*
rnorms(1,ls,ms2) +
dw_s(2,nj,ls,ms2)*
rnorms(2,ls,ms2))*ray
2348 phisij(ni,nj) = phisij(ni,nj) + x
2353 sij(2,:,:) = sij(1,:,:)
2354 sij(6,:,:) = sij(5,:,:)
2360 i = h_mesh%jjs(ni,ms1)
2362 ib = la_h%loc_to_glob(ki,i)
2363 ix = (ki/2)*n_ws1 + ni
2366 j = h_mesh%jjs(nj,ms1)
2368 jb = la_h%loc_to_glob(kj,j)
2369 jx = (kj/2)*n_ws1 + nj
2372 mat_loc1(ix,jx) = hsij(1,ni,nj)
2373 mat_loc2(ix,jx) = hsij(1,ni,nj)
2374 ELSE IF (ki*kj==9)
THEN 2375 mat_loc1(ix,jx) = hsij(6,ni,nj)
2376 mat_loc2(ix,jx) = hsij(6,ni,nj)
2377 ELSE IF (ki*kj==3)
THEN 2378 mat_loc1(ix,jx) = hsij(4,ni,nj)
2379 mat_loc2(ix,jx) = hsij(4,ni,nj)
2385 j = phi_mesh%jj(nj,m2)
2386 jb = la_phi%loc_to_glob(1,j)
2389 mat_loc1(ix,jx) = sij(2*ki-1,ni,nj)
2390 mat_loc2(ix,jx) = sij(2*ki-1,ni,nj)
2394 CALL matsetvalues(h_p_phi_mat1, 2*n_ws1, idxn(1:2*n_ws1), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2395 mat_loc1(1:2*n_ws1,1:2*n_ws1+n_w2), add_values, ierr)
2396 CALL matsetvalues(h_p_phi_mat2, 2*n_ws1, idxn(1:2*n_ws1), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2397 mat_loc2(1:2*n_ws1,1:2*n_ws1+n_w2), add_values, ierr)
2402 i = phi_mesh%jj(ni,m2)
2403 ib = la_phi%loc_to_glob(1,i)
2407 j = h_mesh%jjs(nj,ms1)
2409 jb = la_h%loc_to_glob(kj,j)
2410 jx = (kj/2)*n_ws1 + nj
2412 mat_loc1(ix,jx) = sij(2*kj-1,nj,ni)
2413 mat_loc2(ix,jx) = sij(2*kj-1,nj,ni)
2418 j = phi_mesh%jj(nj,m2)
2419 jb = la_phi%loc_to_glob(1,j)
2422 mat_loc1(ix,jx) = phisij(ni,nj)
2423 mat_loc2(ix,jx) = phisij(ni,nj)
2426 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2427 mat_loc1(1:n_w2,1:2*n_ws1+n_w2), add_values, ierr)
2428 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2429 mat_loc2(1:n_w2,1:2*n_ws1+n_w2), add_values, ierr)
2440 IF (.NOT.
PRESENT(index_fourier) .OR. .NOT.
PRESENT(r_fourier))
RETURN 2441 IF (r_fourier.GT.0.d0)
THEN 2443 DO ms = 1, phi_mesh%mes
2444 IF (phi_mesh%sides(ms) /= index_fourier) cycle
2448 DO ls = 1, phi_mesh%gauss%l_Gs
2451 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms))* phi_mesh%gauss%wws(:,ls))
2453 x = c_mu_phi*
rjs(ls,ms)*ray/r_fourier
2455 DO ni=1, phi_mesh%gauss%n_ws
2456 DO nj=1, phi_mesh%gauss%n_ws
2457 phisij(ni,nj) = phisij(ni,nj) + x*
wws(ni,ls)*
wws(nj,ls)
2464 DO ni = 1, phi_mesh%gauss%n_ws
2465 i = phi_mesh%jjs(ni,ms)
2466 ib = la_phi%loc_to_glob(1,i)
2468 DO nj = 1, phi_mesh%gauss%n_ws
2469 j = phi_mesh%jjs(nj,ms)
2470 jb = la_phi%loc_to_glob(1,j)
2474 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
2475 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
2476 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
2477 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
2481 CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
2482 CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
2483 CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
2484 CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
2491 mode, stab, la_h, h_p_phi_mat1, h_p_phi_mat2, sigma_np, sigma)
2498 #include "petsc/finclude/petsc.h" 2502 INTEGER,
DIMENSION(:),
INTENT(IN) :: jj_v_to_H
2503 INTEGER,
DIMENSION(:),
INTENT(IN) :: Dirichlet_bdy_H_sides
2504 INTEGER,
INTENT(IN) :: mode
2505 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
2506 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_np
2507 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN) :: sigma
2509 INTEGER :: ms, ls, ni, nj, i, j, &
2510 n_ws1, n_w1, m1, ki, kj, ib, jb
2511 REAL(KIND=8) :: x, y, hm1
2512 REAL(KIND=8) :: ray, error, stab_colle_h_mu
2513 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: Hsij
2533 REAL(KIND=8) :: c_sym=.0d0
2538 REAL(KIND=8),
DIMENSION(3*H_mesh%gauss%n_w,3*H_mesh%gauss%n_w) :: mat_loc1, mat_loc2
2539 INTEGER ,
DIMENSION(3*H_mesh%gauss%n_w) :: idxn, jdxn
2544 petscerrorcode :: ierr
2545 mat :: h_p_phi_mat1, h_p_phi_mat2
2548 stab_colle_h_mu = stab(3)
2555 n_ws1 = h_mesh%gauss%n_ws
2556 n_w1 = h_mesh%gauss%n_w
2564 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
2565 ms = dirichlet_bdy_h_sides(count)
2567 hm1 = stab_colle_h_mu/(sum(h_mesh%gauss%rjs(:,ms))*
inputs%sigma_min)
2568 m1 = h_mesh%neighs(ms)
2577 DO ls = 1, h_mesh%gauss%l_Gs
2579 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2580 x = hm1*h_mesh%gauss%rjs(ls,ms)*ray
2582 DO ni = 1, h_mesh%gauss%n_ws
2583 DO nj = 1, h_mesh%gauss%n_ws
2584 y = x * h_mesh%gauss%wws(ni,ls)*h_mesh%gauss%wws(nj,ls)
2585 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(h_mesh%gauss%rnorms(2,ls,ms)**2)
2586 hsij(4,ni,nj) = hsij(4,ni,nj) - y*h_mesh%gauss%rnorms(1,ls,ms)*h_mesh%gauss%rnorms(2,ls,ms)
2587 hsij(5,ni,nj) = hsij(5,ni,nj) + y
2588 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(h_mesh%gauss%rnorms(1,ls,ms)**2)
2603 i = h_mesh%jjs(ni,ms)
2604 ib = la_h%loc_to_glob(ki,i)
2605 ix = ni + (ki-1)*n_ws1
2609 j = h_mesh%jjs(nj,ms)
2610 jb = la_h%loc_to_glob(kj,j)
2611 jx = nj + (kj-1)*n_ws1
2613 IF ((ki == 1) .AND. (kj == 1))
THEN 2614 mat_loc1(ix,jx) = hsij(1,ni,nj)
2615 mat_loc2(ix,jx) = hsij(1,ni,nj)
2616 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 2617 mat_loc1(ix,jx) = hsij(4,ni,nj)
2618 mat_loc2(ix,jx) = hsij(4,ni,nj)
2619 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 2620 mat_loc1(ix,jx) = hsij(4,nj,ni)
2621 mat_loc2(ix,jx) = hsij(4,nj,ni)
2622 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 2623 mat_loc1(ix,jx) = hsij(5,ni,nj)
2624 mat_loc2(ix,jx) = hsij(5,ni,nj)
2625 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 2626 mat_loc1(ix,jx) = hsij(6,ni,nj)
2627 mat_loc2(ix,jx) = hsij(6,ni,nj)
2634 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2635 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2636 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2637 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2648 DO ls = 1, h_mesh%gauss%l_Gs
2650 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
2653 IF (jj_v_to_h(h_mesh%jj(1,m1)) == -1)
THEN 2654 x = h_mesh%gauss%rjs(ls,ms)*ray/sigma(m1)
2656 x = h_mesh%gauss%rjs(ls,ms)*ray/sum(sigma_np(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2663 y = x*h_mesh%gauss%wws(ni,ls)*h_mesh%gauss%wws(nj,ls)
2664 hsij(2,ni,nj) = hsij(2,ni,nj) + y * (-mode/ray)*(
rnorms(1,ls,ms))
2665 hsij(3,ni,nj) = hsij(3,ni,nj) + y * mode/ray *(
rnorms(1,ls,ms))
2666 hsij(5,ni,nj) = hsij(5,ni,nj) + y * (-1/ray) *(
rnorms(1,ls,ms))
2667 hsij(8,ni,nj) = hsij(8,ni,nj) + y * (-mode/ray)*(
rnorms(2,ls,ms))
2668 hsij(9,ni,nj) = hsij(9,ni,nj) + y * mode/ray *(
rnorms(2,ls,ms))
2682 i = h_mesh%jjs(ni,ms)
2683 ib = la_h%loc_to_glob(ki,i)
2684 ix = ni + (ki-1)*n_ws1
2688 j = h_mesh%jjs(nj,ms)
2689 jb = la_h%loc_to_glob(kj,j)
2690 jx = nj + (kj-1)*n_ws1
2692 IF ( (ki == 2) .AND. (kj == 1))
THEN 2693 mat_loc1(ix,jx) = hsij(2,ni,nj)
2694 mat_loc2(ix,jx) = hsij(3,ni,nj)
2695 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 2696 mat_loc1(ix,jx) = hsij(5,ni,nj)
2697 mat_loc2(ix,jx) = hsij(5,ni,nj)
2698 ELSEIF ( (ki == 2) .AND. (kj == 3))
THEN 2699 mat_loc1(ix,jx) = hsij(8,ni,nj)
2700 mat_loc2(ix,jx) = hsij(9,ni,nj)
2707 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2708 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2709 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2710 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2718 i = h_mesh%jjs(ni,ms)
2719 ib = la_h%loc_to_glob(ki,i)
2720 ix = ni + (ki-1)*n_ws1
2724 j = h_mesh%jjs(nj,ms)
2725 jb = la_h%loc_to_glob(kj,j)
2726 jx = nj + (kj-1)*n_ws1
2728 IF ( (kj == 2) .AND. (ki == 1))
THEN 2729 mat_loc1(ix,jx) = hsij(2,nj,ni)
2730 mat_loc2(ix,jx) = hsij(3,nj,ni)
2731 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN 2732 mat_loc1(ix,jx) = hsij(5,nj,ni)
2733 mat_loc2(ix,jx) = hsij(5,nj,ni)
2734 ELSEIF ( (kj == 2) .AND. (ki == 3))
THEN 2735 mat_loc1(ix,jx) = hsij(8,nj,ni)
2736 mat_loc2(ix,jx) = hsij(9,nj,ni)
2742 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2743 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2744 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2745 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2751 DO ls = 1, h_mesh%gauss%l_Gs
2753 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
2756 IF (jj_v_to_h(h_mesh%jj(1,m1)) == -1)
THEN 2757 x = h_mesh%gauss%rjs(ls,ms)*ray/(sigma(m1))
2759 x = h_mesh%gauss%rjs(ls,ms)*ray/(sum(sigma_np(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls)))
2765 y = x*h_mesh%gauss%wws(ni,ls)
2767 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(-h_mesh%gauss%dw_s(2,nj,ls,ms))*(
rnorms(2,ls,ms))
2768 hsij(4,ni,nj) = hsij(4,ni,nj) + y* h_mesh%gauss%dw_s(1,nj,ls,ms) *(
rnorms(2,ls,ms))
2769 hsij(5,ni,nj) = hsij(5,ni,nj) + &
2770 y*(-h_mesh%gauss%dw_s(2,nj,ls,ms)*(
rnorms(2,ls,ms))-h_mesh%gauss%dw_s(1,nj,ls,ms)*(
rnorms(1,ls,ms)))
2771 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(-h_mesh%gauss%dw_s(1,nj,ls,ms))*(
rnorms(1,ls,ms))
2772 hsij(7,ni,nj) = hsij(7,ni,nj) + y* h_mesh%gauss%dw_s(2,nj,ls,ms) *(
rnorms(1,ls,ms))
2787 i = h_mesh%jjs(ni,ms)
2788 ib = la_h%loc_to_glob(ki,i)
2789 ix = ni + (ki-1)*n_ws1
2793 j = h_mesh%jj(nj,m1)
2794 jb = la_h%loc_to_glob(kj,j)
2795 jx = nj + (kj-1)*n_w1
2797 IF ((ki == 1) .AND. (kj == 1))
THEN 2798 mat_loc1(ix,jx) = hsij(1,ni,nj)
2799 mat_loc2(ix,jx) = hsij(1,ni,nj)
2800 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 2801 mat_loc1(ix,jx) = hsij(4,ni,nj)
2802 mat_loc2(ix,jx) = hsij(4,ni,nj)
2803 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 2804 mat_loc1(ix,jx) = hsij(5,ni,nj)
2805 mat_loc2(ix,jx) = hsij(5,ni,nj)
2806 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 2807 mat_loc1(ix,jx) = hsij(6,ni,nj)
2808 mat_loc2(ix,jx) = hsij(6,ni,nj)
2809 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 2810 mat_loc1(ix,jx) = hsij(7,ni,nj)
2811 mat_loc2(ix,jx) = hsij(7,ni,nj)
2818 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
2819 mat_loc1(1:3*n_ws1,1:3*n_w1), add_values, ierr)
2820 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
2821 mat_loc2(1:3*n_ws1,1:3*n_w1), add_values, ierr)
2829 i = h_mesh%jj(ni,m1)
2830 ib = la_h%loc_to_glob(ki,i)
2831 ix = ni + (ki-1)*n_w1
2835 j = h_mesh%jjs(nj,ms)
2836 jb = la_h%loc_to_glob(kj,j)
2837 jx = nj + (kj-1)*n_ws1
2839 IF ((kj == 1) .AND. (ki == 1))
THEN 2840 mat_loc1(ix,jx) = hsij(1,nj,ni)
2841 mat_loc2(ix,jx) = hsij(1,nj,ni)
2842 ELSEIF ((kj == 1) .AND. (ki == 3))
THEN 2843 mat_loc1(ix,jx) = hsij(4,nj,ni)
2844 mat_loc2(ix,jx) = hsij(4,nj,ni)
2845 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN 2846 mat_loc1(ix,jx) = hsij(5,nj,ni)
2847 mat_loc2(ix,jx) = hsij(5,nj,ni)
2848 ELSEIF ((kj == 3) .AND. (ki == 3))
THEN 2849 mat_loc1(ix,jx) = hsij(6,nj,ni)
2850 mat_loc2(ix,jx) = hsij(6,nj,ni)
2851 ELSEIF ((kj == 3) .AND. (ki == 1))
THEN 2852 mat_loc1(ix,jx) = hsij(7,nj,ni)
2853 mat_loc2(ix,jx) = hsij(7,nj,ni)
2860 CALL matsetvalues(h_p_phi_mat1, 3*n_w1 , idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
2861 mat_loc1(1:3*n_w1,1:3*n_ws1), add_values, ierr)
2862 CALL matsetvalues(h_p_phi_mat2, 3*n_w1 , idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
2863 mat_loc2(1:3*n_w1,1:3*n_ws1), add_values, ierr)
2867 CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
2868 CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
2869 CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
2870 CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
2880 SUBROUTINE courant_int_by_parts(H_mesh,phi_mesh,interface_H_phi,sigma,mu_phi,mu_H_field,time,mode,&
2881 rhs_h,nl, la_h, la_phi, vb_1, vb_2, b_ext, sigma_curl_gauss, j_over_sigma_gauss)
2890 #include "petsc/finclude/petsc.h" 2893 TYPE(
mesh_type),
INTENT(IN) :: H_mesh, phi_mesh
2895 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
2896 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma
2897 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
2898 INTEGER,
INTENT(IN) :: mode
2899 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: nl
2900 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: B_ext
2901 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rhs_H
2902 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: J_over_sigma_gauss
2903 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: sigma_curl_gauss
2905 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: src_H
2906 REAL(KIND=8),
DIMENSION(phi_mesh%np,2) :: src_phi
2912 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_ws,phi_mesh%gauss%l_Gs) :: w_cs
2913 REAL(KIND=8),
DIMENSION(2) :: gaussp
2915 INTEGER :: m, l, i, ni, k, ms, ls, n_ws1, n_ws2, ms1, ms2, H_bloc_size, n_w2, m1
2917 REAL(KIND=8),
DIMENSION(6) :: JsolH_anal, rhs_Hl
2918 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%n_w) :: dwH
2921 REAL(KIND=8) :: ray_rjl, muhl
2926 INTEGER,
DIMENSION(H_mesh%np) :: idxn_H
2927 INTEGER,
DIMENSION(phi_mesh%np) :: idxn_phi
2930 REAL(KIND=8),
DIMENSION(6) :: B_ext_l
2931 petscerrorcode :: ierr
2936 CALL veczeroentries(vb_1, ierr)
2937 CALL veczeroentries(vb_2, ierr)
2949 mesh_id1 = h_mesh%i_d(m)
2950 DO l = 1, h_mesh%gauss%l_G
2953 muhl=sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
2955 dwh = h_mesh%gauss%dw(:,:,l,m)
2958 b_ext_l(k) = sum(b_ext(h_mesh%jj(:,m),k)*h_mesh%gauss%ww(:,l))
2964 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
2965 gaussp = gaussp + h_mesh%rr(:,i)*h_mesh%gauss%ww(ni,l)
2966 jsolh_anal(:) = jsolh_anal(:) + muhl*nl(i,:)*h_mesh%gauss%ww(ni,l)
2967 rhs_hl(:) = rhs_hl(:) + rhs_h(i,:)*h_mesh%gauss%ww(ni,l)
2970 ray_rjl = h_mesh%gauss%rj(l,m)*ray
2972 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 2974 jsolh_anal(k) = jsolh_anal(k) + j_over_sigma_gauss(index,k) + sigma_curl_gauss(index,k)
2979 jsolh_anal(k) = jsolh_anal(k) + &
2980 jexact_gauss(k, gaussp, mode, mu_phi, sigma(m), muhl, time, mesh_id1, b_ext_l)/sigma(m)
2984 DO ni = 1,h_mesh%gauss%n_w
2989 src_h(i,1) = src_h(i,1)+ ray_rjl &
2990 *(jsolh_anal(3)*dwh(2,ni) &
2991 + mode/ray*jsolh_anal(6)*h_mesh%gauss%ww(ni,l) &
2992 + rhs_hl(1)*h_mesh%gauss%ww(ni,l))
2994 src_h(i,2) = src_h(i,2)+ ray_rjl &
2995 *(jsolh_anal(4)*dwh(2,ni) &
2996 - mode/ray*jsolh_anal(5)*h_mesh%gauss%ww(ni,l) &
2997 + rhs_hl(2)*h_mesh%gauss%ww(ni,l))
3000 src_h(i,3) = src_h(i,3)+ ray_rjl &
3001 * (-jsolh_anal(1)*dwh(2,ni) &
3002 + 1/ray*jsolh_anal(5)*(ray*dwh(1,ni) + h_mesh%gauss%ww(ni,l)) &
3003 + rhs_hl(3)*h_mesh%gauss%ww(ni,l))
3005 src_h(i,4) = src_h(i,4)+ ray_rjl &
3006 * (-jsolh_anal(2)*dwh(2,ni) &
3007 + 1/ray*jsolh_anal(6)*(ray*dwh(1,ni) + h_mesh%gauss%ww(ni,l)) &
3008 + rhs_hl(4)*h_mesh%gauss%ww(ni,l))
3011 src_h(i,5) = src_h(i,5)+ ray_rjl* &
3012 (-mode/ray*jsolh_anal(2)*h_mesh%gauss%ww(ni,l) &
3013 - jsolh_anal(3)*dwh(1,ni) &
3014 + rhs_hl(5)*h_mesh%gauss%ww(ni,l))
3016 src_h(i,6) = src_h(i,6)+ ray_rjl* &
3017 (mode/ray*jsolh_anal(1)*h_mesh%gauss%ww(ni,l) &
3018 - jsolh_anal(4)*dwh(1,ni) &
3019 + rhs_hl(6)*h_mesh%gauss%ww(ni,l))
3076 CALL gauss(phi_mesh)
3078 n_ws1 = h_mesh%gauss%n_ws
3079 n_ws2 = phi_mesh%gauss%n_ws
3080 n_w2 = phi_mesh%gauss%n_w
3082 h_bloc_size = h_mesh%np
3084 IF (interface_h_phi%mes /=0)
THEN 3085 IF (h_mesh%gauss%n_ws ==
n_ws)
THEN 3089 w_cs(1,ls)=
wws(1,ls)+0.5*
wws(3,ls)
3090 w_cs(2,ls)=
wws(2,ls)+0.5*
wws(3,ls)
3098 DO ms = 1, interface_h_phi%mes
3100 ms2 = interface_h_phi%mesh2(ms)
3101 ms1 = interface_h_phi%mesh1(ms)
3102 m = phi_mesh%neighs(ms2)
3103 m1 = h_mesh%neighs(ms1)
3104 mesh_id1 = h_mesh%i_d(m1)
3107 muhl=sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
3110 b_ext_l(k) = sum(b_ext(interface_h_phi%jjs1(1:n_ws1,ms),k)*w_cs(1:n_ws1,ls))
3115 DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,interface_h_phi%mesh2(ms))
3116 ray = ray + phi_mesh%rr(1,i)*
wws(ni,ls)
3121 i=phi_mesh%jjs(ni,ms2)
3122 gaussp = gaussp + phi_mesh%rr(:,i)*phi_mesh%gauss%wws(ni,ls)
3126 jsolh_anal(k) = jexact_gauss(k, gaussp, mode, mu_phi ,sigma(m1), &
3127 muhl, time, mesh_id1, b_ext_l)/sigma(m1) &
3128 + muhl * sum(nl(h_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
3148 i = interface_h_phi%jjs1(ni,ms)
3149 src_h(i,1) = src_h(i,1)+
rjs(ls,ms2)*ray*( &
3150 -jsolh_anal(3)*w_cs(ni,ls)*(-
rnorms(2,ls,ms2)))
3152 src_h(i,2) = src_h(i,2)+
rjs(ls,ms2)*ray*( &
3153 -jsolh_anal(4)*w_cs(ni,ls)*(-
rnorms(2,ls,ms2)))
3155 src_h(i,3) = src_h(i,3)+
rjs(ls,ms2)*ray*( &
3156 jsolh_anal(1)*w_cs(ni,ls)*(-
rnorms(2,ls,ms2)) &
3157 -jsolh_anal(5)*w_cs(ni,ls)*(-
rnorms(1,ls,ms2)))
3159 src_h(i,4) = src_h(i,4)+
rjs(ls,ms2)*ray*( &
3160 jsolh_anal(2)*w_cs(ni,ls)*(-
rnorms(2,ls,ms2)) &
3161 -jsolh_anal(6)*w_cs(ni,ls)*(-
rnorms(1,ls,ms2)))
3163 src_h(i,5) = src_h(i,5)+
rjs(ls,ms2)*ray*( &
3164 jsolh_anal(3)*w_cs(ni,ls)*(-
rnorms(1,ls,ms2)))
3166 src_h(i,6) = src_h(i,6)+
rjs(ls,ms2)*ray*( &
3167 jsolh_anal(4)*w_cs(ni,ls)*(-
rnorms(1,ls,ms2)))
3173 i = interface_h_phi%jjs2(ni,ms)
3176 src_phi(i,1) = src_phi(i,1)+
rjs(ls,ms2)*( &
3177 - mode*jsolh_anal(2)*
wws(ni,ls) *
rnorms(2,ls,ms2) &
3178 + mode*jsolh_anal(6)*
wws(ni,ls) *
rnorms(1,ls,ms2))
3180 src_phi(i,2) = src_phi(i,2)+
rjs(ls,ms2)*( &
3181 + mode*jsolh_anal(1)*
wws(ni,ls) *
rnorms(2,ls,ms2) &
3182 - mode*jsolh_anal(5)*
wws(ni,ls) *
rnorms(1,ls,ms2))
3188 i = phi_mesh%jj(ni,m)
3189 src_phi(i,1) = src_phi(i,1)+
rjs(ls,ms2)*ray*( &
3190 + jsolh_anal(3) *(
dw_s(2,ni,ls,ms2) *
rnorms(1,ls,ms2)&
3193 src_phi(i,2) = src_phi(i,2)+
rjs(ls,ms2)*ray*( &
3194 + jsolh_anal(4)*(
dw_s(2,ni,ls,ms2) *
rnorms(1,ls,ms2)&
3202 i = interface_h_phi%jjs1(ni,ms)
3203 rhs_hl(:) = rhs_hl(:) + rhs_h(i,:)*w_cs(ni,ls)
3207 i = interface_h_phi%jjs2(ni,ms)
3208 src_phi(i,1) = src_phi(i,1)+
rjs(ls,ms2)*ray*
wws(ni,ls)*( &
3209 rhs_hl(1)*
rnorms(1,ls,ms2) + rhs_hl(5)*
rnorms(2,ls,ms2))
3210 src_phi(i,2) = src_phi(i,2)+
rjs(ls,ms2)*ray*
wws(ni,ls)*( &
3211 rhs_hl(2)*
rnorms(1,ls,ms2) + rhs_hl(6)*
rnorms(2,ls,ms2))
3219 IF (h_mesh%np /= 0)
THEN 3221 idxn_h = la_h%loc_to_glob(1,:)-1
3222 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,1), add_values, ierr)
3223 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,2), add_values, ierr)
3224 idxn_h = la_h%loc_to_glob(2,:)-1
3225 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,4), add_values, ierr)
3226 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,3), add_values, ierr)
3227 idxn_h = la_h%loc_to_glob(3,:)-1
3228 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,5), add_values, ierr)
3229 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,6), add_values, ierr)
3232 IF (phi_mesh%np /=0)
THEN 3234 idxn_phi = la_phi%loc_to_glob(1,:)-1
3235 CALL vecsetvalues(vb_1, phi_mesh%np, idxn_phi, src_phi(:,1), add_values, ierr)
3236 CALL vecsetvalues(vb_2, phi_mesh%np, idxn_phi, src_phi(:,2), add_values, ierr)
3240 CALL vecassemblybegin(vb_1,ierr)
3241 CALL vecassemblyend(vb_1,ierr)
3242 CALL vecassemblybegin(vb_2,ierr)
3243 CALL vecassemblyend(vb_2,ierr)
3254 SUBROUTINE surf_int(H_mesh, phi_mesh, pmag_mesh, interface_H_phi, interface_H_mu, &
3255 list_dirichlet_sides_h, sigma,mu_phi, mu_h_field, time, mode, &
3256 la_h, la_phi, la_pmag, vb_1, vb_2, sigma_tot_gauss, r_fourier, index_fourier)
3262 #include "petsc/finclude/petsc.h" 3265 TYPE(
mesh_type),
INTENT(IN) :: H_mesh, phi_mesh, pmag_mesh
3266 TYPE(
interface_type),
INTENT(IN) :: interface_H_phi, interface_H_mu
3267 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_dirichlet_sides_H
3268 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
3269 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN):: sigma
3270 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
3271 INTEGER,
INTENT(IN) :: mode
3272 REAL(KIND=8),
DIMENSION(:,:) :: sigma_tot_gauss
3273 REAL(KIND=8),
OPTIONAL :: R_fourier
3274 INTEGER,
OPTIONAL :: index_fourier
3275 REAL(KIND=8),
DIMENSION(H_mesh%np, 6) :: src_H
3276 REAL(KIND=8),
DIMENSION(phi_mesh%np, 2) :: src_phi
3277 REAL(KIND=8),
DIMENSION(pmag_mesh%np, 2) :: src_pmag
3279 REAL(KIND=8),
DIMENSION(4,H_mesh%gauss%l_Gs):: Banal
3280 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%l_Gs):: rloc
3281 REAL(KIND=8),
DIMENSION(H_mesh%gauss%l_Gs) :: B_dot_n_cos, B_dot_n_sin, muloc
3282 REAL(KIND=8),
DIMENSION(4,pmag_mesh%gauss%l_Gs):: Banal_pmesh
3283 REAL(KIND=8),
DIMENSION(2,pmag_mesh%gauss%l_Gs):: rloc_pmesh
3284 REAL(KIND=8),
DIMENSION(pmag_mesh%gauss%l_Gs) :: B_dot_n_cos_pmesh, B_dot_n_sin_pmesh, muloc_pmesh
3285 REAL(KIND=8) :: ray, muhl, norm, y
3286 INTEGER :: ms, ls, ns, i, k, m, n, ni, count, index
3287 REAL(KIND=8),
DIMENSION(2) :: gaussp
3288 REAL(KIND=8),
DIMENSION(6) :: EsolH_anal, Esolphi_anal
3289 INTEGER,
DIMENSION(H_mesh%np) :: idxn_H
3290 INTEGER,
DIMENSION(phi_mesh%np) :: idxn_phi
3291 INTEGER,
DIMENSION(pmag_mesh%np) :: idxn_pmag
3293 petscerrorcode :: ierr
3316 m = h_mesh%neighs(ms)
3318 DO ls = 1, h_mesh%gauss%l_Gs
3319 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3321 rloc(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3322 rloc(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3324 banal(1,:) = muhl*hexact(h_mesh, 1, rloc, mode, muloc, time)
3325 banal(2,:) = muhl*hexact(h_mesh, 2, rloc, mode, muloc, time)
3326 banal(3,:) = muhl*hexact(h_mesh, 5, rloc, mode, muloc, time)
3327 banal(4,:) = muhl*hexact(h_mesh, 6, rloc, mode, muloc, time)
3328 b_dot_n_cos = banal(1,:)*h_mesh%gauss%rnorms(1,:,ms) + banal(3,:)*h_mesh%gauss%rnorms(2,:,ms)
3329 b_dot_n_sin = banal(2,:)*h_mesh%gauss%rnorms(1,:,ms) + banal(4,:)*h_mesh%gauss%rnorms(2,:,ms)
3332 DO ls = 1, h_mesh%gauss%l_Gs
3335 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3340 DO ni = 1, h_mesh%gauss%n_ws; i = h_mesh%jjs(ni,ms)
3341 ray = ray + h_mesh%rr(1,i)* h_mesh%gauss%wws(ni,ls)
3344 IF (ray.LT.1.d-12*h_mesh%global_diameter) cycle
3347 DO ns=1, h_mesh%gauss%n_ws
3349 gaussp = gaussp + h_mesh%rr(:,i)*h_mesh%gauss%wws(ns,ls)
3357 esolh_anal(k) = eexact_gauss(k,gaussp,mode,mu_phi,sigma_tot_gauss(index,mod(k+1,2)+1),muhl, time)
3375 DO ns=1, h_mesh%gauss%n_ws
3376 i = h_mesh%jjs(ns,ms)
3377 src_h(i,1) = src_h(i,1)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3378 -esolh_anal(3)*h_mesh%gauss%wws(ns,ls)* &
3379 (h_mesh%gauss%rnorms(2,ls,ms)))
3381 src_h(i,2) = src_h(i,2)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3382 -esolh_anal(4)*h_mesh%gauss%wws(ns,ls)* &
3383 (h_mesh%gauss%rnorms(2,ls,ms)))
3385 src_h(i,3) = src_h(i,3)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3386 esolh_anal(1)*h_mesh%gauss%wws(ns,ls)* &
3387 (h_mesh%gauss%rnorms(2,ls,ms)) - &
3388 esolh_anal(5)*h_mesh%gauss%wws(ns,ls) * &
3389 (h_mesh%gauss%rnorms(1,ls,ms)))
3391 src_h(i,4) = src_h(i,4)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3392 esolh_anal(2)*h_mesh%gauss%wws(ns,ls) * &
3393 (h_mesh%gauss%rnorms(2,ls,ms)) - &
3394 esolh_anal(6)*h_mesh%gauss%wws(ns,ls) * &
3395 (h_mesh%gauss%rnorms(1,ls,ms)))
3397 src_h(i,5) = src_h(i,5)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3398 esolh_anal(3)*h_mesh%gauss%wws(ns,ls)* &
3399 (h_mesh%gauss%rnorms(1,ls,ms)))
3401 src_h(i,6) = src_h(i,6)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3402 esolh_anal(4)*h_mesh%gauss%wws(ns,ls) * &
3403 (h_mesh%gauss%rnorms(1,ls,ms)))
3408 norm = (
inputs%stab(1)/
inputs%Rem)*(sum(h_mesh%gauss%rjs(:,ms))/h_mesh%global_diameter)**(2*
alpha-1)&
3409 /(h_mesh%global_diameter*
inputs%sigma_min*
inputs%mu_min**2)
3413 DO ns = 1, h_mesh%gauss%n_ws
3414 i = h_mesh%jjs(ns,ms)
3415 y = norm*muhl*h_mesh%gauss%rjs(ls,ms)*ray
3416 src_h(i,1) = src_h(i,1) + y*b_dot_n_cos(ls)*h_mesh%gauss%wws(ns,ls)*h_mesh%gauss%rnorms(1,ls,ms)
3417 src_h(i,2) = src_h(i,2) + y*b_dot_n_sin(ls)*h_mesh%gauss%wws(ns,ls)*h_mesh%gauss%rnorms(1,ls,ms)
3418 src_h(i,5) = src_h(i,5) + y*b_dot_n_cos(ls)*h_mesh%gauss%wws(ns,ls)*h_mesh%gauss%rnorms(2,ls,ms)
3419 src_h(i,6) = src_h(i,6) + y*b_dot_n_sin(ls)*h_mesh%gauss%wws(ns,ls)*h_mesh%gauss%rnorms(2,ls,ms)
3429 m = pmag_mesh%neighs(ms)
3431 DO ls = 1, pmag_mesh%gauss%l_Gs
3432 muhl = sum(mu_h_field(pmag_mesh%jjs(:,ms))*pmag_mesh%gauss%wws(:,ls))
3433 muloc_pmesh(ls) = muhl
3434 rloc_pmesh(1,ls) = sum(pmag_mesh%rr(1,pmag_mesh%jjs(:,ms))*pmag_mesh%gauss%wws(:,ls))
3435 rloc_pmesh(2,ls) = sum(pmag_mesh%rr(2,pmag_mesh%jjs(:,ms))*pmag_mesh%gauss%wws(:,ls))
3437 banal_pmesh(1,:) = muhl*hexact(pmag_mesh, 1, rloc_pmesh, mode, muloc_pmesh, time)
3438 banal_pmesh(2,:) = muhl*hexact(pmag_mesh, 2, rloc_pmesh, mode, muloc_pmesh, time)
3439 banal_pmesh(3,:) = muhl*hexact(pmag_mesh, 5, rloc_pmesh, mode, muloc_pmesh, time)
3440 banal_pmesh(4,:) = muhl*hexact(pmag_mesh, 6, rloc_pmesh, mode, muloc_pmesh, time)
3441 b_dot_n_cos_pmesh = banal_pmesh(1,:)*pmag_mesh%gauss%rnorms(1,:,ms) &
3442 + banal_pmesh(3,:)*pmag_mesh%gauss%rnorms(2,:,ms)
3443 b_dot_n_sin_pmesh = banal_pmesh(2,:)*pmag_mesh%gauss%rnorms(1,:,ms) &
3444 + banal_pmesh(4,:)*pmag_mesh%gauss%rnorms(2,:,ms)
3447 DO ls = 1, pmag_mesh%gauss%l_Gs
3449 muhl = sum(mu_h_field(pmag_mesh%jjs(:,ms))*pmag_mesh%gauss%wws(:,ls))
3454 DO ni = 1, pmag_mesh%gauss%n_ws; i = pmag_mesh%jjs(ni,ms)
3455 ray = ray + pmag_mesh%rr(1,i)* pmag_mesh%gauss%wws(ni,ls)
3458 IF (ray.LT.1.d-12*h_mesh%global_diameter) cycle
3462 DO ns=1, pmag_mesh%gauss%n_ws
3463 i = pmag_mesh%jjs(ns,ms)
3464 src_pmag(i,1) = src_pmag(i,1) -
inputs%stab(1)*b_dot_n_cos_pmesh(ls)* &
3465 pmag_mesh%gauss%wws(ns,ls)*pmag_mesh%gauss%rjs(ls,ms)*ray
3466 src_pmag(i,2) = src_pmag(i,2) -
inputs%stab(1)*b_dot_n_sin_pmesh(ls)* &
3467 pmag_mesh%gauss%wws(ns,ls)*pmag_mesh%gauss%rjs(ls,ms)*ray
3477 m = phi_mesh%neighs(ms)
3478 DO ls = 1, phi_mesh%gauss%l_Gs
3481 DO ni = 1, phi_mesh%gauss%n_ws; i = phi_mesh%jjs(ni,ms)
3482 ray = ray + phi_mesh%rr(1,i)* phi_mesh%gauss%wws(ni,ls)
3486 DO ns=1, phi_mesh%gauss%n_ws
3487 i=phi_mesh%jjs(ns,ms)
3488 gaussp = gaussp + phi_mesh%rr(:,i)*phi_mesh%gauss%wws(ns,ls)
3494 esolphi_anal(k) = eexact_gauss(k,gaussp,mode,mu_phi,sigma(1),mu_h_field(1),time)
3501 DO ns=1, phi_mesh%gauss%n_ws
3502 i = phi_mesh%jjs(ns,ms)
3503 DO n = 1, phi_mesh%gauss%n_w
3504 IF (phi_mesh%jj(n,m) == i)
EXIT 3507 src_phi(i,1) = src_phi(i,1)-phi_mesh%gauss%rjs(ls,ms)*ray*( &
3508 +esolphi_anal(3)*(phi_mesh%gauss%dw_s(2,n,ls,ms)*phi_mesh%gauss%rnorms(1,ls,ms) &
3509 -phi_mesh%gauss%dw_s(1,n,ls,ms)*phi_mesh%gauss%rnorms(2,ls,ms))) &
3510 -phi_mesh%gauss%rjs(ls,ms)*(&
3511 -mode*esolphi_anal(2)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(2,ls,ms) &
3512 +mode*esolphi_anal(6)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(1,ls,ms))
3514 src_phi(i,2) = src_phi(i,2)-phi_mesh%gauss%rjs(ls,ms)*ray*( &
3515 +esolphi_anal(4)*(phi_mesh%gauss%dw_s(2,n,ls,ms)*phi_mesh%gauss%rnorms(1,ls,ms) &
3516 -phi_mesh%gauss%dw_s(1,n,ls,ms)*phi_mesh%gauss%rnorms(2,ls,ms))) &
3517 -phi_mesh%gauss%rjs(ls,ms)*(&
3518 mode*esolphi_anal(1)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(2,ls,ms) &
3519 -mode*esolphi_anal(5)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(1,ls,ms))
3536 IF (
PRESENT(r_fourier))
THEN 3537 IF (r_fourier.GE.0.d0)
CALL error_petsc(
'maxwell_update_time_with_H: R_fourier should be -1')
3560 IF (h_mesh%mes/=0)
THEN 3562 idxn_h = la_h%loc_to_glob(1,:)-1
3563 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,1), add_values, ierr)
3564 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,2), add_values, ierr)
3565 idxn_h = la_h%loc_to_glob(2,:)-1
3566 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,4), add_values, ierr)
3567 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,3), add_values, ierr)
3568 idxn_h = la_h%loc_to_glob(3,:)-1
3569 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,5), add_values, ierr)
3570 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,6), add_values, ierr)
3573 idxn_pmag = la_pmag%loc_to_glob(1,:)-1
3574 CALL vecsetvalues(vb_1, pmag_mesh%np, idxn_pmag, src_pmag(:,1), add_values, ierr)
3575 CALL vecsetvalues(vb_2, pmag_mesh%np, idxn_pmag, src_pmag(:,2), add_values, ierr)
3580 IF (phi_mesh%mes/=0)
THEN 3582 idxn_phi = la_phi%loc_to_glob(1,:)-1
3583 CALL vecsetvalues(vb_1, phi_mesh%np, idxn_phi, src_phi(:,1), add_values, ierr)
3584 CALL vecsetvalues(vb_2, phi_mesh%np, idxn_phi, src_phi(:,2), add_values, ierr)
3588 CALL vecassemblybegin(vb_1,ierr)
3589 CALL vecassemblyend(vb_1,ierr)
3590 CALL vecassemblybegin(vb_2,ierr)
3591 CALL vecassemblyend(vb_2,ierr)
3596 count=index_fourier; count=interface_h_mu%mes; count=interface_h_phi%mes
3597 count=
SIZE(list_dirichlet_sides_h)
3601 SUBROUTINE mat_maxwell_mu(H_mesh, jj_v_to_H, interface_H_mu, mode, stab, &
3602 mu_h_field, sigma, la_h, h_p_phi_mat1, h_p_phi_mat2, sigma_np)
3607 #include "petsc/finclude/petsc.h" 3611 INTEGER,
DIMENSION(:),
INTENT(IN) :: jj_v_to_H
3613 INTEGER,
INTENT(IN) :: mode
3614 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
3615 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma, mu_H_field
3616 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_np
3617 INTEGER :: ms, ls, ni, nj, i, j, &
3618 n_ws1, n_ws2, n_w2, n_w1, m1, m2, ki, kj,ib,jb, ms1, ms2
3619 REAL(KIND=8) :: x, y, z, norm, hm1
3620 REAL(KIND=8) :: ray, stab_colle_H_mu
3621 LOGICAL :: mark=.false.
3622 REAL(KIND=8),
DIMENSION(9,SIZE(H_mesh%jj,1),SIZE(H_mesh%jj,1),2,2) :: Hsij, Gsij
3637 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_ws,H_mesh%gauss%l_Gs) :: w_cs
3638 REAL(KIND=8),
DIMENSION(2, H_mesh%gauss%n_w, H_mesh%gauss%l_Gs, H_mesh%mes) :: dw_cs
3639 REAL(KIND=8),
DIMENSION(2, H_mesh%gauss%n_w) :: dwsi,dwsj
3640 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%l_Gs) :: gauss1, gauss2
3646 REAL(KIND=8),
DIMENSION(2) :: normi, normj
3647 REAL(KIND=8),
DIMENSION(SIZE(H_mesh%jjs,1)) :: wwsi, wwsj
3648 INTEGER :: n_wsi, n_wsj, ci, cj, n_wi, n_wj
3651 REAL(KIND=8) :: ref, diff, mu_H, muhl1, muhl2, muhi, muhj, sigmai, sigmaj
3652 REAL(KIND=8) :: sigmal1, sigmal2
3654 REAL(KIND=8) :: c_sym =.0d0
3655 REAL(KIND=8) :: wwiwwj, normt, stab_div
3660 REAL(KIND=8),
DIMENSION(6*H_mesh%gauss%n_w,6*H_mesh%gauss%n_w) :: mat_loc1, mat_loc2
3661 INTEGER ,
DIMENSION(6*H_mesh%gauss%n_w) :: idxn, jdxn
3665 mat :: h_p_phi_mat1, h_p_phi_mat2
3666 petscerrorcode :: ierr
3669 stab_colle_h_mu = stab(3)
3679 n_ws1 = h_mesh%gauss%n_ws
3680 n_ws2 = h_mesh%gauss%n_ws
3681 n_w1 = h_mesh%gauss%n_w
3682 n_w2 = h_mesh%gauss%n_w
3694 DO ms = 1, interface_h_mu%mes
3696 ms2 = interface_h_mu%mesh2(ms)
3697 m2 = h_mesh%neighs(ms2)
3698 ms1 = interface_h_mu%mesh1(ms)
3699 m1 = h_mesh%neighs(ms1)
3701 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
3702 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(1,ms2)))**2)
3703 IF (diff/ref .LT. 1.d-10)
THEN 3707 w_cs(1,ls)=
wws(2,ls)
3708 w_cs(2,ls)=
wws(1,ls)
3709 IF (n_ws1==3) w_cs(n_ws1,ls) =
wws(n_ws1,ls)
3710 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 3715 gauss2(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))*h_mesh%gauss%wws(:,ls))
3716 gauss2(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms2))*h_mesh%gauss%wws(:,ls))
3717 gauss1(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms1))*h_mesh%gauss%wws(:,ls))
3718 gauss1(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms1))*h_mesh%gauss%wws(:,ls))
3722 ref = sqrt(1.d-8+sum(gauss2(:,ls2)**2))
3725 diff = sqrt(sum((gauss2(:,ls2)-gauss1(:,ls1))**2))
3726 IF (diff .LT. 1.d-10)
THEN 3727 dw_cs(:,:,ls2,ms1) = h_mesh%gauss%dw_s(:,:,ls1,ms1)
3732 IF (.NOT.mark)
WRITE(*,*)
' BUG ' 3737 DO ms = 1, interface_h_mu%mes
3739 ms2 = interface_h_mu%mesh2(ms)
3740 ms1 = interface_h_mu%mesh1(ms)
3741 m2 = h_mesh%neighs(ms2)
3742 m1 = h_mesh%neighs(ms1)
3743 mu_h = sum(mu_h_field(h_mesh%jj(:,m1)))/h_mesh%gauss%n_w
3746 hm1 = 1/sum(
rjs(:,ms2))
3758 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
3759 x = hm1*
rjs(ls,ms2)*ray
3762 muhl1 = sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
3763 muhl2 = sum(mu_h_field(h_mesh%jjs(:,ms2))*
wws(:,ls))
3766 normt =stab_colle_h_mu/
inputs%sigma_min
3774 norm = stab_div*(sum(
rjs(:,ms2))/h_mesh%global_diameter)**(2*
alpha)/(
inputs%sigma_min*
inputs%mu_min**2)
3803 wwiwwj = x * wwsi(ni)*wwsj(nj)
3806 z = norm * muhi * muhj * wwiwwj
3808 hsij(1,ni,nj,ci,cj) = hsij(1,ni,nj,ci,cj) + y*normi(2)*normj(2) &
3809 + z*normi(1)*normj(1)
3810 hsij(4,ni,nj,ci,cj) = hsij(4,ni,nj,ci,cj) - y*normj(1)*normi(2) &
3811 + z*normi(1)*normj(2)
3812 hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) + y*(normi(1)*normj(1) + normi(2)*normj(2))
3813 hsij(6,ni,nj,ci,cj) = hsij(6,ni,nj,ci,cj) + y*normi(1)*normj(1) &
3814 + z*normi(2)*normj(2)
3827 i = interface_h_mu%jjs1(ni,ms)
3829 i = interface_h_mu%jjs2(ni,ms)
3831 ib = la_h%loc_to_glob(ki,i)
3832 ix = ni + n_ws1*((ki-1) + 3*(ci-1))
3839 j = interface_h_mu%jjs1(nj,ms)
3841 j = interface_h_mu%jjs2(nj,ms)
3843 jb = la_h%loc_to_glob(kj,j)
3844 jx = nj + n_ws2*((kj-1) + 3*(cj-1))
3846 IF ((ki == 1) .AND. (kj == 1))
THEN 3847 mat_loc1(ix,jx) = hsij(1,ni,nj,ci,cj)
3848 mat_loc2(ix,jx) = hsij(1,ni,nj,ci,cj)
3849 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 3850 mat_loc1(ix,jx) = hsij(4,ni,nj,ci,cj)
3851 mat_loc2(ix,jx) = hsij(4,ni,nj,ci,cj)
3852 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 3853 mat_loc1(ix,jx) = hsij(4,nj,ni,cj,ci)
3854 mat_loc2(ix,jx) = hsij(4,nj,ni,cj,ci)
3855 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 3856 mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)
3857 mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)
3858 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 3859 mat_loc1(ix,jx) = hsij(6,ni,nj,ci,cj)
3860 mat_loc2(ix,jx) = hsij(6,ni,nj,ci,cj)
3869 CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
3870 mat_loc1(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
3871 CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
3872 mat_loc2(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
3881 DO ls = 1, h_mesh%gauss%l_Gs
3884 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
3887 IF (jj_v_to_h(h_mesh%jj(1,m1)) == -1)
THEN 3890 sigmal1 = sum(sigma_np(h_mesh%jjs(:,ms1))*w_cs(:,ls))
3892 IF (jj_v_to_h(h_mesh%jj(1,m2)) == -1)
THEN 3895 sigmal2 = sum(sigma_np(h_mesh%jjs(:,ms2))*
wws(:,ls))
3937 y = x*wwsi(ni)*wwsj(nj)/(2*sigmaj)
3938 hsij(2,ni,nj,ci,cj) = hsij(2,ni,nj,ci,cj) + y * (-mode/ray)*normi(1)
3939 hsij(3,ni,nj,ci,cj) = hsij(3,ni,nj,ci,cj) + y * mode/ray *normi(1)
3940 hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) + y * (-1/ray) *normi(1)
3941 hsij(8,ni,nj,ci,cj) = hsij(8,ni,nj,ci,cj) + y * (-mode/ray)*normi(2)
3942 hsij(9,ni,nj,ci,cj) = hsij(9,ni,nj,ci,cj) + y * mode/ray *normi(2)
3943 y = x*wwsi(ni)*wwsj(nj)/(2*sigmai)
3944 gsij(2,ni,nj,ci,cj) = gsij(2,ni,nj,ci,cj) + y * (-mode/ray)*normj(1)
3945 gsij(3,ni,nj,ci,cj) = gsij(3,ni,nj,ci,cj) + y * ( mode/ray)*normj(1)
3946 gsij(5,ni,nj,ci,cj) = gsij(5,ni,nj,ci,cj) + y * (-1/ray) *normj(1)
3947 gsij(8,ni,nj,ci,cj) = gsij(8,ni,nj,ci,cj) + y * (-mode/ray)*normj(2)
3948 gsij(9,ni,nj,ci,cj) = gsij(9,ni,nj,ci,cj) + y * mode/ray *normj(2)
3966 i = interface_h_mu%jjs1(ni,ms)
3968 i = interface_h_mu%jjs2(ni,ms)
3970 ib = la_h%loc_to_glob(ki,i)
3971 ix = ni + n_wsi*((ki-1) + 3*(ci-1))
3977 j = interface_h_mu%jjs1(nj,ms)
3979 j = interface_h_mu%jjs2(nj,ms)
3981 jb = la_h%loc_to_glob(kj,j)
3982 jx = nj + n_wsj*((kj-1) + 3*(cj-1))
3984 IF ((ki == 2) .AND. (kj == 1))
THEN 3985 mat_loc1(ix,jx) = hsij(2,ni,nj,ci,cj)
3986 mat_loc2(ix,jx) = hsij(3,ni,nj,ci,cj)
3987 ELSE IF((ki == 1) .AND. (kj == 2))
THEN 3988 mat_loc1(ix,jx) = gsij(2,ni,nj,ci,cj)
3989 mat_loc2(ix,jx) = gsij(3,ni,nj,ci,cj)
3990 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 3991 mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)+gsij(5,ni,nj,ci,cj)
3992 mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)+gsij(5,ni,nj,ci,cj)
3993 ELSEIF ((ki == 2) .AND. (kj == 3))
THEN 3994 mat_loc1(ix,jx) = hsij(8,ni,nj,ci,cj)
3995 mat_loc2(ix,jx) = hsij(9,ni,nj,ci,cj)
3996 ELSEIF ((ki == 3) .AND. (kj == 2))
THEN 3997 mat_loc1(ix,jx) = gsij(8,ni,nj,ci,cj)
3998 mat_loc2(ix,jx) = gsij(9,ni,nj,ci,cj)
4007 CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
4008 mat_loc1(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
4009 CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
4010 mat_loc2(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
4015 DO ls = 1, h_mesh%gauss%l_Gs
4018 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
4021 IF (jj_v_to_h(h_mesh%jj(1,m1)) == -1)
THEN 4024 sigmal1 = sum(sigma_np(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4026 IF (jj_v_to_h(h_mesh%jj(1,m2)) == -1)
THEN 4029 sigmal2 = sum(sigma_np(h_mesh%jjs(:,ms2))*
wws(:,ls))
4037 dwsi = dw_cs(:,:,ls,ms1)
4047 dwsi =
dw_s(:,:,ls,ms2)
4059 dwsj = dw_cs(:,:,ls,ms1)
4069 dwsj =
dw_s(:,:,ls,ms2)
4081 y = x*wwsi(ni)/(2*sigmaj)
4082 hsij(1,ni,nj,ci,cj) = hsij(1,ni,nj,ci,cj) + y*(-dwsj(2,nj))*normi(2)
4083 hsij(4,ni,nj,ci,cj) = hsij(4,ni,nj,ci,cj) + y* dwsj(1,nj) *normi(2)
4084 hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) + y*(-dwsj(2,nj) *normi(2)-dwsj(1,nj)*normi(1))
4085 hsij(6,ni,nj,ci,cj) = hsij(6,ni,nj,ci,cj) + y*(-dwsj(1,nj))*normi(1)
4086 hsij(7,ni,nj,ci,cj) = hsij(7,ni,nj,ci,cj) + y* dwsj(2,nj) *normi(1)
4091 y = x*wwsj(nj)/(2*sigmai)
4092 gsij(1,ni,nj,ci,cj) = gsij(1,ni,nj,ci,cj) + y*(-dwsi(2,ni))*normj(2)
4093 gsij(4,ni,nj,ci,cj) = gsij(4,ni,nj,ci,cj) + y* dwsi(2,ni) *normj(1)
4094 gsij(5,ni,nj,ci,cj) = gsij(5,ni,nj,ci,cj) + y*(-dwsi(2,ni) *normj(2)-dwsi(1,ni)*normj(1))
4095 gsij(6,ni,nj,ci,cj) = gsij(6,ni,nj,ci,cj) + y*(-dwsi(1,ni))*normj(1)
4096 gsij(7,ni,nj,ci,cj) = gsij(7,ni,nj,ci,cj) + y* dwsi(1,ni) *normj(2)
4114 i = interface_h_mu%jjs1(ni,ms)
4116 i = interface_h_mu%jjs2(ni,ms)
4118 ib = la_h%loc_to_glob(ki,i)
4119 ix = ni + n_wsi*((ki-1) + 3*(ci-1))
4125 j = h_mesh%jj(nj,m1)
4127 j = h_mesh%jj(nj,m2)
4129 jb = la_h%loc_to_glob(kj,j)
4130 jx = nj + n_wj*((kj-1) + 3*(cj-1))
4132 IF ((ki == 1) .AND. (kj == 1))
THEN 4133 mat_loc1(ix,jx) = hsij(1,ni,nj,ci,cj)
4134 mat_loc2(ix,jx) = hsij(1,ni,nj,ci,cj)
4135 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 4136 mat_loc1(ix,jx) = hsij(4,ni,nj,ci,cj)
4137 mat_loc2(ix,jx) = hsij(4,ni,nj,ci,cj)
4138 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 4139 mat_loc1(ix,jx) = hsij(7,ni,nj,ci,cj)
4140 mat_loc2(ix,jx) = hsij(7,ni,nj,ci,cj)
4141 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 4142 mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)
4143 mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)
4144 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 4145 mat_loc1(ix,jx) = hsij(6,ni,nj,ci,cj)
4146 mat_loc2(ix,jx) = hsij(6,ni,nj,ci,cj)
4156 CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_w2, jdxn(1:6*n_w2), &
4157 mat_loc1(1:6*n_ws1,1:6*n_w2), add_values, ierr)
4158 CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_w2, jdxn(1:6*n_w2), &
4159 mat_loc2(1:6*n_ws1,1:6*n_w2), add_values, ierr)
4167 i = h_mesh%jj(ni,m1)
4169 i = h_mesh%jj(ni,m2)
4171 ib = la_h%loc_to_glob(ki,i)
4172 ix = ni + n_wi*((ki-1) + 3*(ci-1))
4178 j = interface_h_mu%jjs1(nj,ms)
4180 j = interface_h_mu%jjs2(nj,ms)
4182 jb = la_h%loc_to_glob(kj,j)
4183 jx = nj + n_wsj*((kj-1) + 3*(cj-1))
4185 IF ((ki == 1) .AND. (kj == 1))
THEN 4186 mat_loc1(ix,jx) = gsij(1,ni,nj,ci,cj)
4187 mat_loc2(ix,jx) = gsij(1,ni,nj,ci,cj)
4188 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 4189 mat_loc1(ix,jx) = gsij(4,ni,nj,ci,cj)
4190 mat_loc2(ix,jx) = gsij(4,ni,nj,ci,cj)
4191 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 4192 mat_loc1(ix,jx) = gsij(7,ni,nj,ci,cj)
4193 mat_loc2(ix,jx) = gsij(7,ni,nj,ci,cj)
4194 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 4195 mat_loc1(ix,jx) = gsij(5,ni,nj,ci,cj)
4196 mat_loc2(ix,jx) = gsij(5,ni,nj,ci,cj)
4197 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 4198 mat_loc1(ix,jx) = gsij(6,ni,nj,ci,cj)
4199 mat_loc2(ix,jx) = gsij(6,ni,nj,ci,cj)
4208 CALL matsetvalues(h_p_phi_mat1, 6*n_w1, idxn(1:6*n_w1), 6*n_ws2, jdxn(1:6*n_ws2), &
4209 mat_loc1(1:6*n_w1,1:6*n_ws2), add_values, ierr)
4210 CALL matsetvalues(h_p_phi_mat2, 6*n_w1, idxn(1:6*n_w1), 6*n_ws2, jdxn(1:6*n_ws2), &
4211 mat_loc2(1:6*n_w1,1:6*n_ws2), add_values, ierr)
4215 CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
4216 CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
4217 CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
4218 CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
4225 SUBROUTINE courant_mu(H_mesh,interface_H_mu,sigma,mu_H_field,time,mode,nl, &
4226 la_h, vb_1, vb_2, b_ext, j_over_sigma_gauss, sigma_curl_gauss)
4233 #include "petsc/finclude/petsc.h" 4238 REAL(KIND=8),
INTENT(IN) :: time
4239 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN) :: sigma
4240 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
4241 INTEGER,
INTENT(IN) :: mode
4242 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: nl
4243 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: B_ext
4244 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: J_over_sigma_gauss
4245 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: sigma_curl_gauss
4246 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: src_H
4248 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_ws,H_mesh%gauss%l_Gs) :: w_cs
4249 REAL(KIND=8),
DIMENSION(2) :: normi, gaussp1, gaussp2
4250 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_ws) :: wwsi
4251 REAL(KIND=8) :: x, ray
4252 INTEGER :: i, ni, ms, k, ls, n_ws1, n_ws2, ms1, ms2, n_w1, n_w2, m1, m2, ci, n_wsi
4253 INTEGER :: mesh_id1, mesh_id2, index
4254 REAL(KIND=8),
DIMENSION(6) :: JsolH_anal, test, B_ext_l, J_over_sigma_l
4255 REAL(KIND=8) :: muhl1, muhl2, ref, diff
4262 INTEGER,
DIMENSION(H_mesh%np) :: idxn
4265 petscerrorcode :: ierr
4275 n_ws1 = h_mesh%gauss%n_ws
4276 n_ws2 = h_mesh%gauss%n_ws
4277 n_w1 = h_mesh%gauss%n_w
4278 n_w2 = h_mesh%gauss%n_w
4280 DO ms = 1, interface_h_mu%mes
4281 ms1 = interface_h_mu%mesh1(ms)
4282 ms2 = interface_h_mu%mesh2(ms)
4283 m1 = h_mesh%neighs(ms1)
4284 m2 = h_mesh%neighs(ms2)
4286 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
4287 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(1,ms2)))**2)
4288 IF (diff/ref .LT. 1.d-10)
THEN 4292 w_cs(1,ls)=
wws(2,ls)
4293 w_cs(2,ls)=
wws(1,ls)
4294 IF (n_ws1==3) w_cs(n_ws1,ls) =
wws(n_ws1,ls)
4295 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 4301 DO ms = 1, interface_h_mu%mes
4302 ms2 = interface_h_mu%mesh2(ms)
4303 ms1 = interface_h_mu%mesh1(ms)
4304 m2 = h_mesh%neighs(ms2)
4305 m1 = h_mesh%neighs(ms1)
4306 mesh_id1 = h_mesh%i_d(m1)
4307 mesh_id2 = h_mesh%i_d(m2)
4310 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
4316 b_ext_l(k) = sum(b_ext(h_mesh%jjs(:,ms1),k)*h_mesh%gauss%wws(:,ls))
4317 j_over_sigma_l(k) = j_over_sigma_gauss(index,k) + sigma_curl_gauss(index,k)
4319 muhl1=sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4320 gaussp1(1) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms1))*w_cs(:,ls))
4321 gaussp1(2) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms1))*w_cs(:,ls))
4323 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 4325 jsolh_anal(k) = j_over_sigma_l(k) &
4326 + muhl1 *sum(nl(h_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
4330 jsolh_anal(k) = jexact_gauss(k, gaussp1, mode, one ,sigma(m1), &
4331 muhl1, time, mesh_id1, b_ext_l)/sigma(m1) &
4332 + muhl1 * sum(nl(h_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
4339 b_ext_l(k) = sum(b_ext(h_mesh%jjs(:,ms2),k)*h_mesh%gauss%wws(:,ls))
4340 j_over_sigma_l(k) = j_over_sigma_gauss(index,k) + sigma_curl_gauss(index,k)
4342 muhl2=sum(mu_h_field(h_mesh%jjs(:,ms2))*
wws(:,ls))
4343 gaussp2(1) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))*
wws(:,ls))
4344 gaussp2(2) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms2))*
wws(:,ls))
4345 IF (maxval(abs(gaussp1-gaussp2)) > 1.d-11)
THEN 4346 WRITE(*,*)
' BUG courant_mu ' 4349 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 4351 test(k) = j_over_sigma_l(k) &
4352 + muhl2 * sum(nl(h_mesh%jjs(1:n_ws2,ms2),k)*
wws(1:n_ws2,ls))
4353 jsolh_anal(k) = jsolh_anal(k) + test(k)
4357 test(k) = jexact_gauss(k, gaussp2, mode, one ,sigma(m2), &
4358 muhl2, time, mesh_id2, b_ext_l)/sigma(m2) &
4359 + muhl2 * sum(nl(h_mesh%jjs(1:n_ws2,ms2),k)*
wws(1:n_ws2,ls))
4360 jsolh_anal(k) = jsolh_anal(k) + test(k)
4378 i = interface_h_mu%jjs1(ni,ms)
4380 i = interface_h_mu%jjs2(ni,ms)
4382 x =
rjs(ls,ms2)*ray*wwsi(ni)/2
4383 src_h(i,1) = src_h(i,1)+x*(-jsolh_anal(3)*normi(2))
4384 src_h(i,2) = src_h(i,2)+x*(-jsolh_anal(4)*normi(2))
4385 src_h(i,3) = src_h(i,3)+x*(jsolh_anal(1)*normi(2)-jsolh_anal(5)*normi(1))
4386 src_h(i,4) = src_h(i,4)+x*(jsolh_anal(2)*normi(2)-jsolh_anal(6)*normi(1))
4387 src_h(i,5) = src_h(i,5)+x*(jsolh_anal(3)*normi(1))
4388 src_h(i,6) = src_h(i,6)+x*(jsolh_anal(4)*normi(1))
4394 IF (h_mesh%np /= 0)
THEN 4396 idxn = la_h%loc_to_glob(1,:)-1
4397 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,1), add_values, ierr)
4398 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,2), add_values, ierr)
4399 idxn = la_h%loc_to_glob(2,:)-1
4400 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,4), add_values, ierr)
4401 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,3), add_values, ierr)
4402 idxn = la_h%loc_to_glob(3,:)-1
4403 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,5), add_values, ierr)
4404 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,6), add_values, ierr)
4407 CALL vecassemblybegin(vb_1,ierr)
4408 CALL vecassemblyend(vb_1,ierr)
4409 CALL vecassemblybegin(vb_2,ierr)
4410 CALL vecassemblyend(vb_2,ierr)
4414 SUBROUTINE rhs_dirichlet(H_mesh,Dirichlet_bdy_H_sides,sigma,&
4415 mu_h_field,time,mode,nl,stab, la_h, vb_1, vb_2, b_ext, j_over_sigma_bdy, sigma_curl_bdy)
4421 #include "petsc/finclude/petsc.h" 4425 INTEGER,
DIMENSION(:),
INTENT(IN) :: Dirichlet_bdy_H_sides
4426 REAL(KIND=8),
INTENT(IN) :: time
4427 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN) :: sigma
4428 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
4429 INTEGER,
INTENT(IN) :: mode
4430 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: nl
4431 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
4432 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: B_ext
4433 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: J_over_sigma_bdy
4434 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: sigma_curl_bdy
4435 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: src_H
4436 REAL(KIND=8),
DIMENSION(2) :: gaussp1
4437 REAL(KIND=8) :: x, ray, stab_colle_H_mu
4438 INTEGER :: i, ni, ms, k, ls, m1, count
4440 REAL(KIND=8),
DIMENSION(6) :: JsolH_anal, B_ext_l
4441 REAL(KIND=8) :: muhl1, hm1
4442 REAL(KIND=8),
DIMENSION(6,H_mesh%gauss%l_Gs) :: Hloc, Hlocxn
4443 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%l_Gs) :: rloc
4444 REAL(KIND=8),
DIMENSION(1) :: muloc
4452 INTEGER,
DIMENSION(H_mesh%np) :: idxn
4455 petscerrorcode :: ierr
4469 stab_colle_h_mu = stab(3)
4472 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
4473 ms = dirichlet_bdy_h_sides(count)
4475 hm1 = stab_colle_h_mu/(sum(h_mesh%gauss%rjs(:,ms))*
inputs%sigma_min)
4476 m1 = h_mesh%neighs(ms)
4477 mesh_id1 = h_mesh%i_d(m1)
4478 muloc(1) = mu_h_field(h_mesh%jj(1,m1))
4480 DO ls = 1, h_mesh%gauss%l_Gs
4481 rloc(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
4482 rloc(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
4486 hloc(k,:) = hexact(h_mesh, k, rloc, mode, muloc, time)
4489 hlocxn(1,:) = hloc(3,:)*h_mesh%gauss%rnorms(2,:,ms)
4490 hlocxn(2,:) = hloc(4,:)*h_mesh%gauss%rnorms(2,:,ms)
4491 hlocxn(3,:) = hloc(5,:)*h_mesh%gauss%rnorms(1,:,ms)-hloc(1,:)*h_mesh%gauss%rnorms(2,:,ms)
4492 hlocxn(4,:) = hloc(6,:)*h_mesh%gauss%rnorms(1,:,ms)-hloc(2,:)*h_mesh%gauss%rnorms(2,:,ms)
4493 hlocxn(5,:) = -hloc(3,:)*h_mesh%gauss%rnorms(1,:,ms)
4494 hlocxn(6,:) = -hloc(4,:)*h_mesh%gauss%rnorms(1,:,ms)
4496 DO ls = 1, h_mesh%gauss%l_Gs
4501 b_ext_l(k) = sum(b_ext(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls))
4505 muhl1=sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
4506 gaussp1(1) = rloc(1,ls)
4507 gaussp1(2) = rloc(2,ls)
4515 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 4517 jsolh_anal(k) = j_over_sigma_bdy(index,k) &
4518 + sigma_curl_bdy(index,k) &
4519 + muhl1 * sum(nl(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls)) &
4524 jsolh_anal(k) = jexact_gauss(k, gaussp1, mode, one ,sigma(m1), muhl1, &
4525 time, mesh_id1, b_ext_l)/sigma(m1) &
4526 + muhl1 * sum(nl(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls)) &
4532 DO ni = 1, h_mesh%gauss%n_ws
4533 i = h_mesh%jjs(ni,ms)
4534 x = h_mesh%gauss%rjs(ls,ms)*ray*h_mesh%gauss%wws(ni,ls)
4536 src_h(i,1) = src_h(i,1)+x*(-jsolh_anal(3)*h_mesh%gauss%rnorms(2,ls,ms))
4537 src_h(i,2) = src_h(i,2)+x*(-jsolh_anal(4)*h_mesh%gauss%rnorms(2,ls,ms))
4538 src_h(i,3) = src_h(i,3)+x*(jsolh_anal(1)*h_mesh%gauss%rnorms(2,ls,ms)&
4539 -jsolh_anal(5)*h_mesh%gauss%rnorms(1,ls,ms))
4540 src_h(i,4) = src_h(i,4)+x*(jsolh_anal(2)*h_mesh%gauss%rnorms(2,ls,ms)&
4541 -jsolh_anal(6)*h_mesh%gauss%rnorms(1,ls,ms))
4542 src_h(i,5) = src_h(i,5)+x*(jsolh_anal(3)*h_mesh%gauss%rnorms(1,ls,ms))
4543 src_h(i,6) = src_h(i,6)+x*(jsolh_anal(4)*h_mesh%gauss%rnorms(1,ls,ms))
4550 IF (h_mesh%np /= 0)
THEN 4552 idxn = la_h%loc_to_glob(1,:)-1
4553 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,1), add_values, ierr)
4554 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,2), add_values, ierr)
4555 idxn = la_h%loc_to_glob(2,:)-1
4556 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,4), add_values, ierr)
4557 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,3), add_values, ierr)
4558 idxn = la_h%loc_to_glob(3,:)-1
4559 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,5), add_values, ierr)
4560 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,6), add_values, ierr)
4563 CALL vecassemblybegin(vb_1,ierr)
4564 CALL vecassemblyend(vb_1,ierr)
4565 CALL vecassemblybegin(vb_2,ierr)
4566 CALL vecassemblyend(vb_2,ierr)
4575 #include "petsc/finclude/petsc.h" 4580 INTEGER,
POINTER,
DIMENSION(:) :: js_D
4581 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: on_proc_loc, on_proc, not_cav_loc, not_cav
4582 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: is_ok, j_tmp
4583 INTEGER,
DIMENSION(1) :: loc
4584 INTEGER :: m, ms, i, nb_dom, idx, nb_cav, ni
4586 mpi_comm,
INTENT(IN) :: communicator
4588 petscerrorcode :: ierr
4590 IF (
inputs%nb_dom_phi==0)
RETURN 4592 CALL mpi_comm_rank(communicator, rank, ierr)
4594 nb_dom =
inputs%nb_dom_phi
4595 ALLOCATE(on_proc_loc(nb_dom), on_proc(nb_dom))
4596 ALLOCATE(not_cav_loc(nb_dom), not_cav(nb_dom))
4604 IF (minval(abs(
inputs%list_dom_phi-i)) /= 0)
THEN 4605 WRITE(*,*)
'error in dirichlet cavities' 4607 loc = minloc(abs(
inputs%list_dom_phi-i))
4608 on_proc_loc(loc(1)) = rank
4610 IF (mesh%mes /= 0)
THEN 4611 ALLOCATE(is_ok(mesh%mes))
4612 is_ok = mesh%i_d(mesh%neighs)
4613 IF (interface_h_phi%mes /=0)
THEN 4614 is_ok(interface_h_phi%mesh2) = 0
4617 IF (sum(abs(mesh%rr(1,mesh%jjs(:,ms)))) .LT. 1.d-12*mesh%global_diameter)
THEN 4620 IF (
inputs%my_periodic%nb_periodic_pairs /=0)
THEN 4621 IF (minval(abs(
inputs%my_periodic%list_periodic-mesh%sides(ms))) == 0)
THEN 4628 IF (is_ok(ms) == 0) cycle
4630 IF (minval(abs(
inputs%list_dom_phi-i)) /= 0)
THEN 4631 WRITE(*,*)
'error in dirichlet cavities' 4633 loc = minloc(abs(
inputs%list_dom_phi-i))
4634 not_cav_loc(loc(1)) = rank
4637 CALL mpi_allreduce(on_proc_loc, on_proc, nb_dom, mpi_integer, mpi_max, communicator, ierr)
4638 CALL mpi_allreduce(not_cav_loc, not_cav, nb_dom, mpi_integer, mpi_max, communicator, ierr)
4640 ALLOCATE(j_tmp(
SIZE(js_d)+nb_dom))
4641 j_tmp(1:
SIZE(js_d)) = js_d
4644 IF ( (not_cav(i)==-1) .AND. (on_proc(i)==rank) )
THEN 4648 IF (mesh%i_d(m) ==
inputs%list_dom_phi(i))
THEN 4649 DO ni = 1, mesh%gauss%n_w
4650 IF (minval(abs(mesh%jjs-mesh%jj(ni,m))) /=0)
THEN 4651 j_tmp(idx) = mesh%jj(ni,m)
4657 WRITE(*,*)
'add ', j_tmp(idx),
'in dom ',
inputs%list_dom_phi(i),
' : proc ', rank
4658 WRITE(*,*)
'add ', mesh%rr(:,j_tmp(idx)), mesh%i_d(m)
4666 nb_cav = idx -
SIZE(js_d)
4667 IF (nb_cav /= 0)
THEN 4673 DEALLOCATE(on_proc_loc, on_proc, j_tmp)
4674 DEALLOCATE(not_cav_loc, not_cav)
4675 IF (
ALLOCATED(is_ok))
DEALLOCATE(is_ok)
4677 WRITE(*,
'(a,x,i2,x,a,x,i2)')
'I have detected', nb_cav,
' cavity(ies) on proc', rank
4681 SUBROUTINE smb_sigma_prod_curl(communicator, mesh, jj_v_to_H, list_mode, H_in, one_over_sigma_in, sigma_nj_m,&
4689 #include "petsc/finclude/petsc.h" 4693 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
4694 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
4695 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: H_in
4696 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_in
4697 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%me),
INTENT(IN) :: sigma_nj_m
4698 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
4699 REAL(KIND=8),
DIMENSION(:,:,:) :: V_out
4700 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: H_gauss, RotH
4701 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: one_over_sigma_gauss
4702 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: RotH_bar
4703 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
4704 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
4705 INTEGER :: m, l , i, mode, index, k
4706 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,6) :: H_in_loc
4707 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: one_over_sigma_in_loc
4708 REAL(KIND=8) :: ray, sigma_np_gauss
4709 INTEGER :: nb_procs, bloc_size, m_max_pad, code
4710 mpi_comm :: communicator
4712 DO i = 1,
SIZE(list_mode)
4716 j_loc = mesh%jj(:,m)
4718 h_in_loc(:,k) = h_in(j_loc,k,i)
4721 one_over_sigma_in_loc(:,k) = one_over_sigma_in(j_loc,k,i)
4724 DO l = 1, mesh%gauss%l_G
4726 dw_loc = mesh%gauss%dw(:,:,l,m)
4729 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
4732 h_gauss(index,1,i) = sum(h_in_loc(:,1)*mesh%gauss%ww(:,l))
4733 h_gauss(index,3,i) = sum(h_in_loc(:,3)*mesh%gauss%ww(:,l))
4734 h_gauss(index,5,i) = sum(h_in_loc(:,5)*mesh%gauss%ww(:,l))
4736 h_gauss(index,2,i) = sum(h_in_loc(:,2)*mesh%gauss%ww(:,l))
4737 h_gauss(index,4,i) = sum(h_in_loc(:,4)*mesh%gauss%ww(:,l))
4738 h_gauss(index,6,i) = sum(h_in_loc(:,6)*mesh%gauss%ww(:,l))
4741 roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
4742 -sum(h_in_loc(:,3)*dw_loc(2,:))
4743 roth(index,4,i) = sum(h_in_loc(:,2)*dw_loc(2,:)) &
4744 -sum(h_in_loc(:,6)*dw_loc(1,:))
4745 roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
4746 +sum(h_in_loc(:,3)*dw_loc(1,:)) &
4747 -mode/ray*h_gauss(index,2,i)
4749 roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
4750 -sum(h_in_loc(:,4)*dw_loc(2,:))
4751 roth(index,3,i) = sum(h_in_loc(:,1)*dw_loc(2,:)) &
4752 -sum(h_in_loc(:,5)*dw_loc(1,:))
4753 roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
4754 +sum(h_in_loc(:,4)*dw_loc(1,:))&
4755 +mode/ray*h_gauss(index,1,i)
4757 IF (jj_v_to_h(mesh%jj(1,m)) == -1)
THEN 4758 one_over_sigma_gauss(index,1,i) = 0.d0
4759 one_over_sigma_gauss(index,2,i) = 0.d0
4761 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m)
4764 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc(:,1)*mesh%gauss%ww(:,l))
4765 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc(:,2)*mesh%gauss%ww(:,l))
4768 sigma_np_gauss = sum(sigma_nj_m(:,m)*mesh%gauss%ww(:,l))
4770 roth_bar(index,k,i) = roth(index,k,i)/sigma_np_gauss
4776 CALL mpi_comm_size(communicator, nb_procs, code)
4777 bloc_size =
SIZE(one_over_sigma_gauss,1)/nb_procs+1
4778 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
4780 roth, one_over_sigma_gauss, v_out, 1, nb_procs, bloc_size, m_max_pad)
4782 v_out = roth_bar - v_out
4787 h_in, one_over_sigma_in, sigma_np, sigma, v_out)
4794 #include "petsc/finclude/petsc.h" 4798 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
4799 INTEGER,
DIMENSION(:),
INTENT(IN) :: Dirichlet_bdy_H_sides
4800 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
4801 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: H_in
4802 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_in
4803 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_np
4804 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
4805 REAL(KIND=8),
DIMENSION(:,:,:) :: V_out
4806 REAL(KIND=8),
DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),6,SIZE(list_mode)) :: H_gauss, RotH
4807 REAL(KIND=8),
DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),6,SIZE(list_mode)) :: RotH_bar
4808 REAL(KIND=8),
DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),2,SIZE(list_mode)) :: one_over_sigma_gauss
4809 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
4810 INTEGER :: ms, ls , i, mode, index, k
4811 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc
4812 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,6) :: H_in_loc
4813 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: one_over_sigma_in_loc
4815 INTEGER :: nb_procs, bloc_size, m_max_pad, code, count, m1
4816 mpi_comm :: communicator
4818 DO i = 1,
SIZE(list_mode)
4821 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
4822 ms = dirichlet_bdy_h_sides(count)
4823 m1 = mesh%neighs(ms)
4825 j_loc = mesh%jjs(:,ms)
4827 h_in_loc(:,k) = h_in(j_loc,k,i)
4830 one_over_sigma_in_loc(:,k) = one_over_sigma_in(j_loc,k,i)
4833 DO ls = 1, mesh%gauss%l_Gs
4835 dw_loc = mesh%gauss%dw_s(:,:,ls,ms)
4838 ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
4841 h_gauss(index,1,i) = sum(h_in_loc(:,1)*mesh%gauss%wws(:,ls))
4842 h_gauss(index,3,i) = sum(h_in_loc(:,3)*mesh%gauss%wws(:,ls))
4843 h_gauss(index,5,i) = sum(h_in_loc(:,5)*mesh%gauss%wws(:,ls))
4845 h_gauss(index,2,i) = sum(h_in_loc(:,2)*mesh%gauss%wws(:,ls))
4846 h_gauss(index,4,i) = sum(h_in_loc(:,4)*mesh%gauss%wws(:,ls))
4847 h_gauss(index,6,i) = sum(h_in_loc(:,6)*mesh%gauss%wws(:,ls))
4850 roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
4851 -sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(2,:))
4853 roth(index,4,i) = sum(h_in(mesh%jj(:,m1),2,i)*dw_loc(2,:)) &
4854 -sum(h_in(mesh%jj(:,m1),6,i)*dw_loc(1,:))
4856 roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
4857 +sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(1,:)) &
4858 -mode/ray*h_gauss(index,2,i)
4861 roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
4862 -sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(2,:))
4864 roth(index,3,i) = sum(h_in(mesh%jj(:,m1),1,i)*dw_loc(2,:)) &
4865 -sum(h_in(mesh%jj(:,m1),5,i)*dw_loc(1,:))
4867 roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
4868 +sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(1,:))&
4869 +mode/ray*h_gauss(index,1,i)
4871 IF (jj_v_to_h(mesh%jj(1,m1)) == -1)
THEN 4872 one_over_sigma_gauss(index,1,i) = 0.d0
4873 one_over_sigma_gauss(index,2,i) = 0.d0
4875 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m1)
4878 roth_bar(index,k,i) = roth(index,k,i)/sigma(m1)
4881 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc(:,1)*mesh%gauss%wws(:,ls))
4882 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc(:,2)*mesh%gauss%wws(:,ls))
4884 roth_bar(index,k,i) = roth(index,k,i)/sum(sigma_np(mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
4891 IF (
SIZE(dirichlet_bdy_h_sides).GE.1)
THEN 4892 CALL mpi_comm_size(communicator, nb_procs, code)
4893 bloc_size =
SIZE(one_over_sigma_gauss,1)/nb_procs+1
4894 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
4896 roth, one_over_sigma_gauss, v_out, 1, nb_procs, bloc_size, m_max_pad)
4898 v_out = roth_bar - v_out
4904 h_in, one_over_sigma_in, sigma_np, sigma, v_out)
4911 #include "petsc/finclude/petsc.h" 4915 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
4917 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
4918 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: H_in
4919 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_in
4920 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_np
4921 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
4922 REAL(KIND=8),
DIMENSION(:,:,:) :: V_out
4923 REAL(KIND=8),
DIMENSION(2*mesh%gauss%l_Gs*interface_H_mu%mes,6,SIZE(list_mode)) :: H_gauss, RotH
4924 REAL(KIND=8),
DIMENSION(2*mesh%gauss%l_Gs*interface_H_mu%mes,6,SIZE(list_mode)) :: RotH_bar
4925 REAL(KIND=8),
DIMENSION(2*mesh%gauss%l_Gs*interface_H_mu%mes,2,SIZE(list_mode)) :: one_over_sigma_gauss
4926 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
4927 INTEGER :: ms, ls , i, mode, index, k
4928 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc1, j_loc2
4929 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,6) :: H_in_loc1, H_in_loc2
4930 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: one_over_sigma_in_loc1, one_over_sigma_in_loc2
4931 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,mesh%gauss%l_Gs) :: w_cs
4932 REAL(KIND=8) :: ray, diff, ref
4933 INTEGER :: nb_procs, bloc_size, m_max_pad, code
4934 INTEGER :: ms1, ms2, m1, m2, mesh_id1, mesh_id2
4935 mpi_comm :: communicator
4937 DO ms = 1, interface_h_mu%mes
4938 ms1 = interface_h_mu%mesh1(ms)
4939 ms2 = interface_h_mu%mesh2(ms)
4940 m1 = mesh%neighs(ms1)
4941 m2 = mesh%neighs(ms2)
4942 ref = 1.d-8+sum((mesh%rr(:,mesh%jjs(1,ms1)) - mesh%rr(:,mesh%jjs(2,ms1)))**2)
4943 diff = sum((mesh%rr(:,mesh%jjs(1,ms1)) - mesh%rr(:,mesh%jjs(1,ms2)))**2)
4944 IF (diff/ref .LT. 1.d-10)
THEN 4945 w_cs = mesh%gauss%wws
4947 DO ls = 1, mesh%gauss%l_Gs
4948 w_cs(1,ls)= mesh%gauss%wws(2,ls)
4949 w_cs(2,ls)= mesh%gauss%wws(1,ls)
4950 IF (mesh%gauss%n_ws==3) w_cs(mesh%gauss%n_ws,ls) = mesh%gauss%wws(mesh%gauss%n_ws,ls)
4951 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 4956 DO i = 1,
SIZE(list_mode)
4959 DO ms = 1, interface_h_mu%mes
4960 ms2 = interface_h_mu%mesh2(ms)
4961 ms1 = interface_h_mu%mesh1(ms)
4962 m2 = mesh%neighs(ms2)
4963 m1 = mesh%neighs(ms1)
4964 mesh_id1 = mesh%i_d(m1)
4965 mesh_id2 = mesh%i_d(m2)
4966 j_loc1 = mesh%jjs(:,ms1)
4967 j_loc2 = mesh%jjs(:,ms2)
4969 h_in_loc1(:,k) = h_in(j_loc1,k,i)
4970 h_in_loc2(:,k) = h_in(j_loc2,k,i)
4973 one_over_sigma_in_loc1(:,k) = one_over_sigma_in(j_loc1,k,i)
4974 one_over_sigma_in_loc2(:,k) = one_over_sigma_in(j_loc2,k,i)
4977 DO ls = 1, mesh%gauss%l_Gs
4980 dw_loc = mesh%gauss%dw_s(:,:,ls,ms1)
4982 ray = sum(mesh%rr(1,mesh%jjs(:,ms1))*w_cs(:,ls))
4984 h_gauss(index,1,i) = sum(h_in_loc1(:,1)*w_cs(:,ls))
4985 h_gauss(index,3,i) = sum(h_in_loc1(:,3)*w_cs(:,ls))
4986 h_gauss(index,5,i) = sum(h_in_loc1(:,5)*w_cs(:,ls))
4987 h_gauss(index,2,i) = sum(h_in_loc1(:,2)*w_cs(:,ls))
4988 h_gauss(index,4,i) = sum(h_in_loc1(:,4)*w_cs(:,ls))
4989 h_gauss(index,6,i) = sum(h_in_loc1(:,6)*w_cs(:,ls))
4992 roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
4993 -sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(2,:))
4994 roth(index,4,i) = sum(h_in(mesh%jj(:,m1),2,i)*dw_loc(2,:)) &
4995 -sum(h_in(mesh%jj(:,m1),6,i)*dw_loc(1,:))
4996 roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
4997 +sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(1,:)) &
4998 -mode/ray*h_gauss(index,2,i)
5000 roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
5001 -sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(2,:))
5002 roth(index,3,i) = sum(h_in(mesh%jj(:,m1),1,i)*dw_loc(2,:)) &
5003 -sum(h_in(mesh%jj(:,m1),5,i)*dw_loc(1,:))
5004 roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
5005 +sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(1,:))&
5006 +mode/ray*h_gauss(index,1,i)
5008 IF (jj_v_to_h(mesh%jj(1,m1)) == -1)
THEN 5009 one_over_sigma_gauss(index,1,i) = 0.d0
5010 one_over_sigma_gauss(index,2,i) = 0.d0
5012 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m1)
5015 roth_bar(index,k,i) = roth(index,k,i)/sigma(m1)
5018 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc1(:,1)*w_cs(:,ls))
5019 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc1(:,2)*w_cs(:,ls))
5021 roth_bar(index,k,i) = roth(index,k,i)/sum(sigma_np(mesh%jjs(:,ms1))*w_cs(:,ls))
5027 dw_loc = mesh%gauss%dw_s(:,:,ls,ms2)
5029 ray = sum(mesh%rr(1,mesh%jjs(:,ms2))*mesh%gauss%wws(:,ls))
5031 h_gauss(index,1,i) = sum(h_in_loc2(:,1)*mesh%gauss%wws(:,ls))
5032 h_gauss(index,3,i) = sum(h_in_loc2(:,3)*mesh%gauss%wws(:,ls))
5033 h_gauss(index,5,i) = sum(h_in_loc2(:,5)*mesh%gauss%wws(:,ls))
5034 h_gauss(index,2,i) = sum(h_in_loc2(:,2)*mesh%gauss%wws(:,ls))
5035 h_gauss(index,4,i) = sum(h_in_loc2(:,4)*mesh%gauss%wws(:,ls))
5036 h_gauss(index,6,i) = sum(h_in_loc2(:,6)*mesh%gauss%wws(:,ls))
5039 roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
5040 -sum(h_in(mesh%jj(:,m2),3,i)*dw_loc(2,:))
5041 roth(index,4,i) = sum(h_in(mesh%jj(:,m2),2,i)*dw_loc(2,:)) &
5042 -sum(h_in(mesh%jj(:,m2),6,i)*dw_loc(1,:))
5043 roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
5044 +sum(h_in(mesh%jj(:,m2),3,i)*dw_loc(1,:)) &
5045 -mode/ray*h_gauss(index,2,i)
5047 roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
5048 -sum(h_in(mesh%jj(:,m2),4,i)*dw_loc(2,:))
5049 roth(index,3,i) = sum(h_in(mesh%jj(:,m2),1,i)*dw_loc(2,:)) &
5050 -sum(h_in(mesh%jj(:,m2),5,i)*dw_loc(1,:))
5051 roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
5052 +sum(h_in(mesh%jj(:,m2),4,i)*dw_loc(1,:))&
5053 +mode/ray*h_gauss(index,1,i)
5055 IF (jj_v_to_h(mesh%jj(1,m2)) == -1)
THEN 5056 one_over_sigma_gauss(index,1,i) = 0.d0
5057 one_over_sigma_gauss(index,2,i) = 0.d0
5059 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m2)
5062 roth_bar(index,k,i) = roth(index,k,i)/sigma(m2)
5065 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc2(:,1)*mesh%gauss%wws(:,ls))
5066 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc2(:,2)*mesh%gauss%wws(:,ls))
5068 roth_bar(index,k,i) = roth(index,k,i)/sum(sigma_np(mesh%jjs(:,ms2))*mesh%gauss%wws(:,ls))
5075 IF (interface_h_mu%mes.GE.1)
THEN 5076 CALL mpi_comm_size(communicator, nb_procs, code)
5077 bloc_size =
SIZE(one_over_sigma_gauss,1)/nb_procs+1
5078 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5080 roth, one_over_sigma_gauss, v_out, 1, nb_procs, bloc_size, m_max_pad)
5082 v_out = roth_bar - v_out
5088 mu_h_field, mu_phi, one_over_sigma_tot, time, sigma, j_over_sigma_gauss)
5094 #include "petsc/finclude/petsc.h" 5098 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
5099 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5100 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: B_in
5101 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_tot
5102 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
5103 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
5104 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
5105 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: J_over_sigma_gauss
5106 REAL(KIND=8),
DIMENSION(mesh%me*mesh%gauss%l_G,6,SIZE(list_mode)) :: J_exact_gauss
5107 REAL(KIND=8),
DIMENSION(mesh%me*mesh%gauss%l_G,2,SIZE(list_mode)) :: one_over_sigma_gauss
5108 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
5109 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,6) :: B_in_loc
5110 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: one_over_sigma_in_loc
5111 REAL(KIND=8),
DIMENSION(6) :: B_ext_l
5112 REAL(KIND=8) :: muhl, ray
5113 REAL(KIND=8),
DIMENSION(2) :: gaussp
5114 INTEGER :: mode, mesh_id1, k, i, index, l, m, ni
5115 INTEGER :: nb_procs, bloc_size, m_max_pad, code
5116 mpi_comm :: communicator
5118 DO i = 1,
SIZE(list_mode)
5122 mesh_id1 = mesh%i_d(m)
5123 j_loc = mesh%jj(:,m)
5125 b_in_loc(:,k) = b_in(j_loc,k,i)
5128 one_over_sigma_in_loc(:,k) = one_over_sigma_tot(j_loc,k,i)
5130 DO l = 1, mesh%gauss%l_G
5134 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
5137 b_ext_l(1) = sum(b_in_loc(:,1)*mesh%gauss%ww(:,l))
5138 b_ext_l(3) = sum(b_in_loc(:,3)*mesh%gauss%ww(:,l))
5139 b_ext_l(5) = sum(b_in_loc(:,5)*mesh%gauss%ww(:,l))
5140 b_ext_l(2) = sum(b_in_loc(:,2)*mesh%gauss%ww(:,l))
5141 b_ext_l(4) = sum(b_in_loc(:,4)*mesh%gauss%ww(:,l))
5142 b_ext_l(6) = sum(b_in_loc(:,6)*mesh%gauss%ww(:,l))
5144 DO ni = 1, mesh%gauss%n_w
5145 gaussp = gaussp + mesh%rr(:,mesh%jj(ni,m))*mesh%gauss%ww(ni,l)
5147 muhl=sum(mu_h_field(mesh%jj(:,m))*mesh%gauss%ww(:,l))
5150 j_exact_gauss(index,k,i)=jexact_gauss(k, gaussp, mode, mu_phi, sigma(m), &
5151 muhl, time, mesh_id1, b_ext_l)
5154 IF (jj_v_to_h(mesh%jj(1,m)) == -1)
THEN 5155 one_over_sigma_gauss(index,1,i) = 0.d0
5156 one_over_sigma_gauss(index,2,i) = 0.d0
5158 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m)
5161 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc(:,1)*mesh%gauss%ww(:,l))
5162 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc(:,2)*mesh%gauss%ww(:,l))
5169 CALL mpi_comm_size(communicator, nb_procs, code)
5170 bloc_size =
SIZE(j_exact_gauss,1)/nb_procs+1
5171 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5173 j_exact_gauss, one_over_sigma_gauss, j_over_sigma_gauss,&
5174 1, nb_procs, bloc_size, m_max_pad)
5179 mu_h_field, mu_phi, one_over_sigma_tot, time, sigma, j_over_sigma_gauss)
5185 #include "petsc/finclude/petsc.h" 5189 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
5190 INTEGER,
DIMENSION(:),
INTENT(IN) :: Dirichlet_bdy_H_sides
5191 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5192 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: B_in
5193 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_tot
5194 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
5195 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
5196 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
5197 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: J_over_sigma_gauss
5198 REAL(KIND=8),
DIMENSION(SIZE(Dirichlet_bdy_H_sides)*mesh%gauss%l_Gs,6,SIZE(list_mode)) :: J_exact_gauss
5199 REAL(KIND=8),
DIMENSION(SIZE(Dirichlet_bdy_H_sides)*mesh%gauss%l_Gs,2,SIZE(list_mode)) :: one_over_sigma_gauss
5200 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc
5201 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,6) :: B_in_loc
5202 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: one_over_sigma_in_loc
5203 REAL(KIND=8),
DIMENSION(6) :: B_ext_l
5204 REAL(KIND=8) :: muhl, ray
5205 REAL(KIND=8),
DIMENSION(2) :: gaussp
5206 INTEGER :: mode, mesh_id1, k, i, count, index, ls, ms, ni
5207 INTEGER :: nb_procs, bloc_size, m_max_pad, code, m1
5208 mpi_comm :: communicator
5210 DO i = 1,
SIZE(list_mode)
5213 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
5214 ms = dirichlet_bdy_h_sides(count)
5215 m1 = mesh%neighs(ms)
5216 mesh_id1 = mesh%i_d(m1)
5217 j_loc = mesh%jjs(:,ms)
5219 b_in_loc(:,k) = b_in(j_loc,k,i)
5222 one_over_sigma_in_loc(:,k) = one_over_sigma_tot(j_loc,k,i)
5225 DO ls = 1, mesh%gauss%l_Gs
5229 ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
5232 b_ext_l(1) = sum(b_in_loc(:,1)*mesh%gauss%wws(:,ls))
5233 b_ext_l(3) = sum(b_in_loc(:,3)*mesh%gauss%wws(:,ls))
5234 b_ext_l(5) = sum(b_in_loc(:,5)*mesh%gauss%wws(:,ls))
5235 b_ext_l(2) = sum(b_in_loc(:,2)*mesh%gauss%wws(:,ls))
5236 b_ext_l(4) = sum(b_in_loc(:,4)*mesh%gauss%wws(:,ls))
5237 b_ext_l(6) = sum(b_in_loc(:,6)*mesh%gauss%wws(:,ls))
5239 DO ni = 1, mesh%gauss%n_ws
5240 gaussp = gaussp + mesh%rr(:,mesh%jjs(ni,ms))*mesh%gauss%wws(ni,ls)
5242 muhl=sum(mu_h_field(mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
5245 j_exact_gauss(index,k,i)=jexact_gauss(k, gaussp, mode, mu_phi, sigma(m1),&
5246 muhl, time, mesh_id1, b_ext_l)
5249 IF (jj_v_to_h(mesh%jj(1,m1)) == -1)
THEN 5250 one_over_sigma_gauss(index,1,i) = 0.d0
5251 one_over_sigma_gauss(index,2,i) = 0.d0
5253 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m1)
5256 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc(:,1)*mesh%gauss%wws(:,ls))
5257 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc(:,2)*mesh%gauss%wws(:,ls))
5263 IF (
SIZE(dirichlet_bdy_h_sides).GE.1)
THEN 5264 CALL mpi_comm_size(communicator, nb_procs, code)
5265 bloc_size =
SIZE(j_exact_gauss,1)/nb_procs+1
5266 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5268 j_exact_gauss, one_over_sigma_gauss, j_over_sigma_gauss,&
5269 1, nb_procs, bloc_size, m_max_pad)
5275 mu_h_field, mu_phi, one_over_sigma_tot, time, sigma, j_over_sigma_gauss)
5281 #include "petsc/finclude/petsc.h" 5285 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
5287 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5288 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: B_in
5289 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_tot
5290 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
5291 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
5292 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
5293 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: J_over_sigma_gauss
5294 REAL(KIND=8),
DIMENSION(2*interface_H_mu%mes*mesh%gauss%l_Gs,6,SIZE(list_mode)) :: J_exact_gauss
5295 REAL(KIND=8),
DIMENSION(2*interface_H_mu%mes*mesh%gauss%l_Gs,2,SIZE(list_mode)) :: one_over_sigma_gauss
5296 REAL(KIND=8),
DIMENSION(6) :: B_ext_l
5297 REAL(KIND=8) :: muhl, diff, ref, ray
5298 REAL(KIND=8),
DIMENSION(2) :: gaussp
5299 INTEGER :: mode, k, i, mesh_id1, mesh_id2, ni
5300 INTEGER :: nb_procs, bloc_size, m_max_pad, code
5301 INTEGER :: ms, ms1, ms2, m1, m2, ls, index
5302 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc1, j_loc2
5303 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,6) :: B_in_loc1, B_in_loc2
5304 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: one_over_sigma_in_loc1,one_over_sigma_in_loc2
5305 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,mesh%gauss%l_Gs) :: w_cs
5306 mpi_comm :: communicator
5309 DO ms = 1, interface_h_mu%mes
5310 ms1 = interface_h_mu%mesh1(ms)
5311 ms2 = interface_h_mu%mesh2(ms)
5312 m1 = mesh%neighs(ms1)
5313 m2 = mesh%neighs(ms2)
5314 ref = 1.d-8+sum((mesh%rr(:,mesh%jjs(1,ms1)) - mesh%rr(:,mesh%jjs(2,ms1)))**2)
5315 diff = sum((mesh%rr(:,mesh%jjs(1,ms1)) - mesh%rr(:,mesh%jjs(1,ms2)))**2)
5316 IF (diff/ref .LT. 1.d-10)
THEN 5317 w_cs = mesh%gauss%wws
5319 DO ls = 1, mesh%gauss%l_Gs
5320 w_cs(1,ls)= mesh%gauss%wws(2,ls)
5321 w_cs(2,ls)= mesh%gauss%wws(1,ls)
5322 IF (mesh%gauss%n_ws==3) w_cs(mesh%gauss%n_ws,ls) = mesh%gauss%wws(mesh%gauss%n_ws,ls)
5323 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 5328 DO i = 1,
SIZE(list_mode)
5331 DO ms = 1, interface_h_mu%mes
5332 ms2 = interface_h_mu%mesh2(ms)
5333 ms1 = interface_h_mu%mesh1(ms)
5334 m2 = mesh%neighs(ms2)
5335 m1 = mesh%neighs(ms1)
5336 mesh_id1 = mesh%i_d(m1)
5337 mesh_id2 = mesh%i_d(m2)
5338 j_loc1 = mesh%jjs(:,ms1)
5339 j_loc2 = mesh%jjs(:,ms2)
5341 b_in_loc1(:,k) = b_in(j_loc1,k,i)
5342 b_in_loc2(:,k) = b_in(j_loc2,k,i)
5345 one_over_sigma_in_loc1(:,k) = one_over_sigma_tot(j_loc1,k,i)
5346 one_over_sigma_in_loc2(:,k) = one_over_sigma_tot(j_loc2,k,i)
5349 DO ls = 1, mesh%gauss%l_Gs
5353 ray = sum(mesh%rr(1,mesh%jjs(:,ms1))*w_cs(:,ls))
5355 b_ext_l(1) = sum(b_in_loc1(:,1)*w_cs(:,ls))
5356 b_ext_l(3) = sum(b_in_loc1(:,3)*w_cs(:,ls))
5357 b_ext_l(5) = sum(b_in_loc1(:,5)*w_cs(:,ls))
5358 b_ext_l(2) = sum(b_in_loc1(:,2)*w_cs(:,ls))
5359 b_ext_l(4) = sum(b_in_loc1(:,4)*w_cs(:,ls))
5360 b_ext_l(6) = sum(b_in_loc1(:,6)*w_cs(:,ls))
5362 DO ni = 1, mesh%gauss%n_ws
5363 gaussp = gaussp + mesh%rr(:,mesh%jjs(ni,ms1))*w_cs(ni,ls)
5365 muhl=sum(mu_h_field(mesh%jjs(:,ms1))*w_cs(:,ls))
5368 j_exact_gauss(index,k,i)=jexact_gauss(k, gaussp, mode, mu_phi, sigma(m1),&
5369 muhl, time, mesh_id1, b_ext_l)
5372 IF (jj_v_to_h(mesh%jj(1,m1)) == -1)
THEN 5373 one_over_sigma_gauss(index,1,i) = 0.d0
5374 one_over_sigma_gauss(index,2,i) = 0.d0
5376 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m1)
5379 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc1(:,1)*w_cs(:,ls))
5380 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc1(:,2)*w_cs(:,ls))
5386 ray = sum(mesh%rr(1,mesh%jjs(:,ms2))*mesh%gauss%wws(:,ls))
5388 b_ext_l(1) = sum(b_in_loc2(:,1)*mesh%gauss%wws(:,ls))
5389 b_ext_l(3) = sum(b_in_loc2(:,3)*mesh%gauss%wws(:,ls))
5390 b_ext_l(5) = sum(b_in_loc2(:,5)*mesh%gauss%wws(:,ls))
5391 b_ext_l(2) = sum(b_in_loc2(:,2)*mesh%gauss%wws(:,ls))
5392 b_ext_l(4) = sum(b_in_loc2(:,4)*mesh%gauss%wws(:,ls))
5393 b_ext_l(6) = sum(b_in_loc2(:,6)*mesh%gauss%wws(:,ls))
5395 DO ni = 1, mesh%gauss%n_ws
5396 gaussp = gaussp + mesh%rr(:,mesh%jjs(ni,ms2))*mesh%gauss%wws(ni,ls)
5398 muhl=sum(mu_h_field(mesh%jjs(:,ms2))*mesh%gauss%wws(:,ls))
5401 j_exact_gauss(index,k,i)=jexact_gauss(k, gaussp, mode, mu_phi, sigma(m2),&
5402 muhl, time, mesh_id2, b_ext_l)
5405 IF (jj_v_to_h(mesh%jj(1,m2)) == -1)
THEN 5406 one_over_sigma_gauss(index,1,i) = 0.d0
5407 one_over_sigma_gauss(index,2,i) = 0.d0
5409 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m2)
5412 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc2(:,1)*mesh%gauss%wws(:,ls))
5413 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc2(:,2)*mesh%gauss%wws(:,ls))
5419 IF (interface_h_mu%mes.GE.1)
THEN 5420 CALL mpi_comm_size(communicator, nb_procs, code)
5421 bloc_size =
SIZE(j_exact_gauss,1)/nb_procs+1
5422 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5424 j_exact_gauss, one_over_sigma_gauss, j_over_sigma_gauss,&
5425 1, nb_procs, bloc_size, m_max_pad)
5430 SUBROUTINE smb_sigma_neumann(communicator, mesh, Neumann_bdy_H_sides, list_mode, &
5431 one_over_sigma_tot, sigma_tot_gauss_neumann)
5437 #include "petsc/finclude/petsc.h" 5441 INTEGER,
DIMENSION(:),
INTENT(IN) :: Neumann_bdy_H_sides
5442 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5443 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_tot
5444 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: sigma_tot_gauss_Neumann
5445 REAL(KIND=8),
DIMENSION(SIZE(Neumann_bdy_H_sides)*mesh%gauss%l_Gs,2,SIZE(list_mode)):: one_over_sigma_gauss
5446 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc
5447 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: one_over_sigma_in_loc
5448 INTEGER :: mode, k, i, count, index, ls, ms
5449 INTEGER :: nb_procs, bloc_size, m_max_pad, code
5450 mpi_comm :: communicator
5452 DO i = 1,
SIZE(list_mode)
5455 DO count = 1,
SIZE(neumann_bdy_h_sides)
5456 ms = neumann_bdy_h_sides(count)
5457 j_loc = mesh%jjs(:,ms)
5459 one_over_sigma_in_loc(:,k) = one_over_sigma_tot(j_loc,k,i)
5462 DO ls = 1, mesh%gauss%l_Gs
5464 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc(:,1)*mesh%gauss%wws(:,ls))
5465 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc(:,2)*mesh%gauss%wws(:,ls))
5470 sigma_tot_gauss_neumann=one_over_sigma_gauss
5472 IF (
SIZE(neumann_bdy_h_sides).GE.1)
THEN 5473 CALL mpi_comm_size(communicator, nb_procs, code)
5474 bloc_size =
SIZE(one_over_sigma_gauss,1)/nb_procs+1
5475 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5485 vv = x/ max(x*x, 1.d-20)
5492
subroutine solver(my_ksp, b, x, reinit, verbose)
subroutine smb_sigma_prod_curl_inter_mu(communicator, mesh, jj_v_to_H, interface_H_mu, list_mode, H_in, one_over_sigma_in, sigma_np, sigma, V_out)
real(kind=8) function one_over_x(x)
subroutine, public extract(xghost, ks, ke, LA, phi)
integer, dimension(:), allocatable neumann_bdy_h_sides
subroutine, public create_my_ghost(mesh, LA, ifrom)
subroutine surf_int(H_mesh, phi_mesh, pmag_mesh, interface_H_phi, interface_H_mu, list_dirichlet_sides_H, sigma, mu_phi, mu_H_field, time, mode, LA_H, LA_phi, LA_pmag, vb_1, vb_2, sigma_tot_gauss, R_fourier, index_fourier)
subroutine smb_sigma_neumann(communicator, mesh, Neumann_bdy_H_sides, list_mode, one_over_sigma_tot, sigma_tot_gauss_Neumann)
subroutine dirichlet_nodes(jjs_in, sides_in, dir_in, js_d)
subroutine vector_without_bc_glob_js_d(vv_mesh, list_mode, vv_3_LA, vv_mode_global_js_D)
subroutine, public fft_scalar_vect_no_overshoot(communicator, scalar_bounds, V1_in, V2_in, V_out, pb, nb_procs, bloc_size, m_max_pad, temps)
subroutine, public fft_par_scal_funct(communicator, c1_inout, funct, nb_procs, bloc_size, m_max_pad, temps)
real(kind=8), dimension(:,:,:), pointer rnorms
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
subroutine mat_maxwell_mu(H_mesh, jj_v_to_H, interface_H_mu, mode, stab, mu_H_field, sigma, LA_H, H_p_phi_mat1, H_p_phi_mat2, sigma_np)
subroutine error_petsc(string)
real(kind=8) function user_time()
subroutine, public fft_par_cross_prod_dcl(communicator, V1_in, V2_in, V_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine smb_current_over_sigma_inter_mu(communicator, mesh, jj_v_to_H, interface_H_mu, list_mode, B_in, mu_H_field, mu_phi, one_over_sigma_tot, time, sigma, J_over_sigma_gauss)
subroutine scalar_with_bc_glob_js_d(pp_mesh, list_mode, pp_1_LA, pp_js_D, pp_mode_global_js_D)
subroutine smb_current_over_sigma(communicator, mesh, jj_v_to_H, list_mode, B_in, mu_H_field, mu_phi, one_over_sigma_tot, time, sigma, J_over_sigma_gauss)
real(kind=8), parameter, private alpha
subroutine mat_h_p_phi_maxwell(H_mesh, pmag_mesh, phi_mesh, interface_H_phi, mode, mu_H_field, mu_phi, c_mass, stab, R_fourier, index_fourier, LA_H, LA_pmag, LA_phi, H_p_phi_mat1, H_p_phi_mat2, sigma_nj_m, sigma)
subroutine smb_current_over_sigma_bdy(communicator, mesh, jj_v_to_H, Dirichlet_bdy_H_sides, list_mode, B_in, mu_H_field, mu_phi, one_over_sigma_tot, time, sigma, J_over_sigma_gauss)
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
subroutine dirichlet_rhs(js_D, bs_D, b)
subroutine dirichlet_cavities(communicator, interface_H_phi, mesh, js_D)
subroutine dirichlet_nodes_parallel(mesh, list_dirichlet_sides, js_d)
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
subroutine rhs_dirichlet(H_mesh, Dirichlet_bdy_H_sides, sigma, mu_H_field, time, mode, nl, stab, LA_H, vb_1, vb_2, B_ext, J_over_sigma_bdy, sigma_curl_bdy)
subroutine mat_dirichlet_maxwell(H_mesh, jj_v_to_H, Dirichlet_bdy_H_sides, mode, stab, LA_H, H_p_phi_mat1, H_p_phi_mat2, sigma_np, sigma)
subroutine dirichlet_m_parallel(matrix, glob_js_D)
type(my_verbose), public talk_to_me
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
subroutine courant_int_by_parts(H_mesh, phi_mesh, interface_H_phi, sigma, mu_phi, mu_H_field, time, mode, rhs_H, nl, LA_H, LA_phi, vb_1, vb_2, B_ext, sigma_curl_gauss, J_over_sigma_gauss)
integer, dimension(:), allocatable neumann_bdy_pmag_sides
real(kind=8), dimension(:,:), pointer rjs
real(kind=8), dimension(:,:), pointer wws
subroutine courant_mu(H_mesh, interface_H_mu, sigma, mu_H_field, time, mode, nl, LA_H, vb_1, vb_2, B_ext, J_over_sigma_gauss, sigma_curl_gauss)
subroutine smb_sigma_prod_curl(communicator, mesh, jj_v_to_H, list_mode, H_in, one_over_sigma_in, sigma_nj_m, sigma, V_out)
subroutine smb_sigma_prod_curl_bdy(communicator, mesh, jj_v_to_H, Dirichlet_bdy_H_sides, list_mode, H_in, one_over_sigma_in, sigma_np, sigma, V_out)
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)
real(kind=8), dimension(:,:,:,:), pointer dw_s
integer, dimension(:), allocatable neumann_bdy_phi_sides
subroutine, public maxwell_decouple_with_h(comm_one_d, H_mesh, pmag_mesh, phi_mesh, interface_H_phi, interface_H_mu, Hn, Bn, phin, Hn1, Bn1, phin1, vel, stab_in, sigma_in, R_fourier, index_fourier, mu_H_field, mu_phi, time, dt_in, Rem, list_mode, H_phi_per, LA_H, LA_pmag, LA_phi, LA_mhd, one_over_sigma_ns_in, jj_v_to_H)