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
64 REAL(KIND=8),
SAVE :: dt
65 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE,
SAVE :: sigma_ns_bar
67 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE,
SAVE :: sigma_np
68 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE,
SAVE :: sigma
69 REAL(KIND=8),
SAVE :: sigma_min
70 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: H_mode_global_js_D
71 TYPE(
dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: H_global_D
73 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: pmag_mode_global_js_D
74 TYPE(
dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: pmag_global_D
76 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: phi_mode_global_js_D
77 TYPE(
dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: phi_global_D
78 INTEGER,
DIMENSION(:),
POINTER,
SAVE :: pmag_js_D, phi_js_D
79 LOGICAL,
SAVE :: once=.true.
80 INTEGER,
SAVE :: m_max_c
82 REAL(KIND=8),
DIMENSION(3),
SAVE :: stab
83 INTEGER,
SAVE :: my_petscworld_rank
84 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: sigma_curl_gauss_bdy
85 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: J_over_sigma_gauss_bdy
86 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: sigma_curl_gauss_inter_mu
87 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: J_over_sigma_gauss_inter_mu
88 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: sigma_tot_gauss_Neumann
89 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:),
SAVE :: sigma_nj_m
90 REAL(KIND=8),
DIMENSION(H_mesh%gauss%l_G*H_mesh%me,6,SIZE(list_mode)) :: sigma_curl_gauss
91 REAL(KIND=8),
DIMENSION(H_mesh%gauss%l_G*H_mesh%me,6,SIZE(list_mode)) :: J_over_sigma_gauss
92 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),6,SIZE(Hn,3)) :: H_ns
93 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),2,SIZE(Hn,3)) :: one_over_sigma_tot
94 LOGICAL,
ALLOCATABLE,
DIMENSION(:) :: Dir_pmag
95 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: rhs_H
96 REAL(KIND=8),
DIMENSION(phi_mesh%np,2) :: rhs_phi
99 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),6,SIZE(Hn,3)) :: H_pert
101 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),6,SIZE(Hn,3)) :: NL, B_ext
102 REAL(KIND=8),
DIMENSION(3) :: temps_par
103 INTEGER,
POINTER,
DIMENSION(:) :: H_ifrom, pmag_ifrom, phi_ifrom, H_p_phi_ifrom
104 REAL(KIND=8),
DIMENSION(phi_mesh%np, 2) :: phin_p1
105 REAL(KIND=8),
DIMENSION(H_mesh%np, 6) :: Hn_p1
106 LOGICAL,
DIMENSION(H_mesh%mes) :: virgin1
107 LOGICAL,
DIMENSION(phi_mesh%mes) :: virgin2
108 INTEGER :: mode, k, i, n, m, ms, code, nj, j, count
109 INTEGER :: nb_procs, bloc_size, m_max_pad
110 REAL(KIND=8) :: tps, nr_vel, tps_tot, tps_cumul, norm
112 REAL(KIND=8) :: one_and_half
113 DATA one_and_half/1.5d0/
115 petscerrorcode :: ierr
116 mpi_comm,
DIMENSION(:),
POINTER :: comm_one_d
117 mat,
DIMENSION(:),
POINTER,
SAVE :: h_p_phi_mat1, h_p_phi_mat2
118 mat :: tampon1, tampon2, precond1, precond2
119 ksp,
DIMENSION(:),
POINTER,
SAVE :: h_p_phi_ksp1, h_p_phi_ksp2
120 vec,
SAVE :: vx_1, vb_1, vx_1_ghost, vx_2, vb_2, vx_2_ghost
133 CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
140 IF (
inputs%if_quasi_static_approx)
THEN 142 ELSE IF (
inputs%type_pb==
"fhd")
THEN 150 n =
SIZE(h_ifrom)+
SIZE(pmag_ifrom)+
SIZE(phi_ifrom)
151 ALLOCATE(h_p_phi_ifrom(n))
152 IF (
SIZE(h_ifrom)/=0)
THEN 153 h_p_phi_ifrom(1:
SIZE(h_ifrom)) = h_ifrom
155 IF (
SIZE(pmag_ifrom)/=0)
THEN 156 h_p_phi_ifrom(
SIZE(h_ifrom)+1:
SIZE(h_ifrom)+
SIZE(pmag_ifrom)) = pmag_ifrom
158 IF (
SIZE(phi_ifrom)/=0)
THEN 159 h_p_phi_ifrom(
SIZE(h_ifrom)+
SIZE(pmag_ifrom)+1:)=phi_ifrom
162 n = 3*h_mesh%dom_np + pmag_mesh%dom_np + phi_mesh%dom_np
163 CALL veccreateghost(comm_one_d(1), n, &
164 petsc_determine,
SIZE(h_p_phi_ifrom), h_p_phi_ifrom, vx_1, ierr)
165 CALL vecghostgetlocalform(vx_1, vx_1_ghost, ierr)
166 CALL vecduplicate(vx_1, vb_1, ierr)
167 CALL veccreateghost(comm_one_d(1), n, &
168 petsc_determine,
SIZE(h_p_phi_ifrom), h_p_phi_ifrom, vx_2, ierr)
169 CALL vecghostgetlocalform(vx_2, vx_2_ghost, ierr)
170 CALL vecduplicate(vx_2, vb_2, ierr)
174 ALLOCATE(sigma(
SIZE(sigma_in)))
175 sigma = sigma_in * rem
178 CALL mpi_allreduce(minval(sigma),sigma_min,1,mpi_double_precision, mpi_min,comm_one_d(1), ierr)
222 m_max_c =
SIZE(list_mode)
226 ALLOCATE(sigma_nj_m(h_mesh%gauss%n_w,h_mesh%me))
227 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 228 ALLOCATE(sigma_ns_bar(
SIZE(hn,1)))
229 sigma_ns_bar = sigma_bar_in_fourier_space(h_mesh)*rem
233 DO nj = 1, h_mesh%gauss%n_w
235 IF (jj_v_to_h(j) == -1)
THEN 236 sigma_nj_m(nj,m) = sigma(m)
238 sigma_nj_m(nj,m) = sigma_ns_bar(j)
244 sigma_nj_m(:,m) = sigma(m)
247 ALLOCATE(sigma_np(
SIZE(hn,1)))
250 DO nj = 1, h_mesh%gauss%n_w
251 sigma_np(h_mesh%jj(nj,m)) = sigma_nj_m(nj,m)
260 ALLOCATE (dir_pmag(maxval(pmag_mesh%sides)))
262 DO ms = 1,
SIZE(dir_pmag)
263 IF (minval(abs(
inputs%list_dirichlet_sides_H-ms)) == 0)
THEN 264 dir_pmag(ms) = .true.
266 IF (minval(abs(
inputs%list_inter_H_phi-ms)) == 0)
THEN 267 dir_pmag(ms) = .true.
271 CALL dirichlet_nodes(pmag_mesh%jjs, pmag_mesh%sides, dir_pmag, pmag_js_d)
276 ALLOCATE(pmag_global_d(m_max_c))
278 ALLOCATE(pmag_global_d(i)%DRL(
SIZE(pmag_mode_global_js_d(i)%DIL)))
279 pmag_global_d(i)%DRL = 0.d0
287 IF (interface_h_phi%mes/=0)
THEN 288 virgin1(interface_h_phi%mesh1) = .false.
289 virgin2(interface_h_phi%mesh2) = .false.
291 IF (interface_h_mu%mes/=0)
THEN 292 virgin1(interface_h_mu%mesh1) = .false.
293 virgin1(interface_h_mu%mesh2) = .false.
297 DO ms = 1, h_mesh%mes
298 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))).LT.1d-12*h_mesh%global_diameter) cycle
299 IF (.NOT.virgin1(ms)) cycle
300 IF(minval(abs(h_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))==0) cycle
302 IF (
inputs%my_periodic%nb_periodic_pairs /=0)
THEN 303 IF (minval(abs(
inputs%my_periodic%list_periodic-h_mesh%sides(ms))) == 0) cycle
310 DO ms = 1, h_mesh%mes
311 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))).LT.1d-12*h_mesh%global_diameter) cycle
312 IF (.NOT.virgin1(ms)) cycle
313 IF(minval(abs(h_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))==0) cycle
315 IF (
inputs%my_periodic%nb_periodic_pairs /=0)
THEN 316 IF (minval(abs(
inputs%my_periodic%list_periodic-h_mesh%sides(ms))) == 0) cycle
324 DO ms = 1, pmag_mesh%mes
325 IF (maxval(abs(pmag_mesh%rr(1,pmag_mesh%jjs(:,ms)))).LT.1d-12*pmag_mesh%global_diameter) cycle
326 IF(minval(abs(pmag_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))==0) cycle
327 IF(minval(abs(pmag_mesh%sides(ms)-
inputs%list_inter_H_phi))==0) cycle
329 IF (
inputs%my_periodic%nb_periodic_pairs /=0)
THEN 330 IF (minval(abs(
inputs%my_periodic%list_periodic-pmag_mesh%sides(ms))) == 0) cycle
337 DO ms = 1, pmag_mesh%mes
338 IF (maxval(abs(pmag_mesh%rr(1,pmag_mesh%jjs(:,ms)))).LT.1d-12*pmag_mesh%global_diameter) cycle
339 IF(minval(abs(pmag_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))==0) cycle
340 IF(minval(abs(pmag_mesh%sides(ms)-
inputs%list_inter_H_phi))==0) cycle
342 IF (
inputs%my_periodic%nb_periodic_pairs /=0)
THEN 343 IF (minval(abs(
inputs%my_periodic%list_periodic-pmag_mesh%sides(ms))) == 0) cycle
351 DO ms = 1, phi_mesh%mes
353 IF (phi_mesh%sides(ms)==index_fourier) cycle
355 IF (.NOT.virgin2(ms)) cycle
356 IF (maxval(abs(phi_mesh%rr(1,phi_mesh%jjs(:,ms)))).LT.1d-12*phi_mesh%global_diameter) cycle
357 IF (minval(abs(phi_mesh%sides(ms)-
inputs%phi_list_dirichlet_sides))==0) cycle
362 DO ms = 1, phi_mesh%mes
364 IF (phi_mesh%sides(ms)==index_fourier) cycle
366 IF (.NOT.virgin2(ms)) cycle
367 IF (maxval(abs(phi_mesh%rr(1,phi_mesh%jjs(:,ms)))).LT.1d-12*phi_mesh%global_diameter) cycle
368 IF (minval(abs(phi_mesh%sides(ms)-
inputs%phi_list_dirichlet_sides))==0) cycle
377 DO ms = 1, h_mesh%mes
378 IF (minval(abs(h_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))/=0) cycle
379 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))) .LT.1d-12*h_mesh%global_diameter) cycle
384 DO ms = 1, h_mesh%mes
385 IF (minval(abs(h_mesh%sides(ms)-
inputs%list_dirichlet_sides_H))/=0) cycle
386 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))) .LT.1d-12*h_mesh%global_diameter) cycle
392 ALLOCATE(h_global_d(m_max_c))
394 ALLOCATE(h_global_d(i)%DRL(
SIZE(h_mode_global_js_d(i)%DIL)))
403 ALLOCATE(phi_global_d(m_max_c))
405 ALLOCATE(phi_global_d(i)%DRL(
SIZE(phi_mode_global_js_d(i)%DIL)))
406 phi_global_d(i)%DRL = 0.d0
411 ALLOCATE(h_p_phi_mat1(m_max_c),h_p_phi_ksp1(m_max_c))
412 ALLOCATE(h_p_phi_mat2(m_max_c),h_p_phi_ksp2(m_max_c))
418 ALLOCATE(sigma_curl_gauss_bdy(0,0,0))
419 ALLOCATE(j_over_sigma_gauss_bdy(0,0,0))
420 sigma_curl_gauss_bdy = 0.d0
421 j_over_sigma_gauss_bdy = 0.d0
424 IF (interface_h_mu%mes.GE.1)
THEN 425 ALLOCATE(sigma_curl_gauss_inter_mu(2*h_mesh%gauss%l_Gs*interface_h_mu%mes,6,
SIZE(list_mode)))
426 ALLOCATE(j_over_sigma_gauss_inter_mu(2*h_mesh%gauss%l_Gs*interface_h_mu%mes,6,
SIZE(list_mode)))
428 ALLOCATE(sigma_curl_gauss_inter_mu(0,0,0))
429 ALLOCATE(j_over_sigma_gauss_inter_mu(0,0,0))
430 sigma_curl_gauss_inter_mu = 0.d0
431 j_over_sigma_gauss_inter_mu = 0.d0
435 IF(my_petscworld_rank==0)
THEN 436 WRITE(*,*)
"WARNING, Neumann BC: either sigma or CURL(H)_Neumann is axisymmetric." 438 ALLOCATE(sigma_tot_gauss_neumann(h_mesh%gauss%l_Gs*
SIZE(
neumann_bdy_h_sides),2,
SIZE(list_mode)))
440 ALLOCATE(sigma_tot_gauss_neumann(0,0,0))
441 sigma_tot_gauss_neumann = 0.d0
463 mode,mu_h_field, mu_phi, one_and_half/dt, stab, r_fourier, index_fourier, &
464 la_h, la_pmag, la_phi, h_p_phi_mat1(i), h_p_phi_mat2(i), sigma_nj_m, sigma)
469 CALL mat_maxwell_mu(h_mesh, jj_v_to_h, interface_h_mu, mode, stab, &
470 mu_h_field, sigma, la_h, h_p_phi_mat1(i), h_p_phi_mat2(i), sigma_np)
476 stab, mu_h_field, la_h, h_p_phi_mat1(i), h_p_phi_mat2(i), sigma_np, sigma)
478 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 480 h_phi_per%perlist, h_p_phi_mat1(i), la_mhd)
482 h_phi_per%perlist, h_p_phi_mat2(i), la_mhd)
499 CALL init_solver(
inputs%my_par_H_p_phi,h_p_phi_ksp1(i),h_p_phi_mat1(i),comm_one_d(1),&
501 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 IF (h_mesh%me/=0)
THEN 520 IF (
inputs%if_permeability_variable_in_theta)
THEN 521 CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
522 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
523 bloc_size =
SIZE(b_ext,1)/nb_procs+1
525 h_mesh, hn, bn, nb_procs, bloc_size, m_max_pad, time,temps_par)
527 h_mesh, hn1, bn1, nb_procs, bloc_size, m_max_pad, time,temps_par)
531 bn(:,k,i) = mu_h_field*hn(:,k,i)
532 bn1(:,k,i) = mu_h_field*hn1(:,k,i)
542 nr_vel =
norm_sf(comm_one_d,
'L2', h_mesh, list_mode, vel)
544 IF(
inputs%if_quasi_static_approx)
THEN 547 b_ext(:,:,i) = h_b_quasi_static(
'B', h_mesh%rr, mode)
553 IF (nr_vel .LE. 1.d-10)
THEN 555 ELSE IF (
inputs%type_pb==
"fhd")
THEN 558 CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
559 bloc_size =
SIZE(vel,1)/nb_procs+1
560 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
567 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 569 one_over_sigma_tot = 0.d0
571 DO nj = 1, h_mesh%gauss%n_w
574 IF (jj_v_to_h(j) /= -1)
THEN 575 h_ns(j,:,:) = 2*hn(j,:,:)- hn1(j,:,:)
576 one_over_sigma_tot(j,:,:) = one_over_sigma_ns_in(jj_v_to_h(j),:,:)/rem
578 DO i = 1,
SIZE(list_mode)
581 one_over_sigma_tot(j,1,i) = 1.d0/sigma(m)
590 one_over_sigma_tot, sigma_nj_m, sigma, sigma_curl_gauss)
593 one_over_sigma_tot, sigma_np, sigma, sigma_curl_gauss_bdy)
595 sigma_curl_gauss_bdy = 0.d0
597 IF (interface_h_mu%mes.GE.1)
THEN 599 one_over_sigma_tot, sigma_np, sigma, sigma_curl_gauss_inter_mu)
601 sigma_curl_gauss_inter_mu = 0.d0
606 mu_h_field, mu_phi, one_over_sigma_tot, time, sigma, j_over_sigma_gauss)
609 list_mode, b_ext, mu_h_field, mu_phi, one_over_sigma_tot, time, sigma,&
610 j_over_sigma_gauss_bdy)
612 j_over_sigma_gauss_bdy = 0.d0
614 IF (interface_h_mu%mes.GE.1)
THEN 616 list_mode, b_ext, mu_h_field, mu_phi, one_over_sigma_tot, time, sigma,&
617 j_over_sigma_gauss_inter_mu)
619 j_over_sigma_gauss_inter_mu=0.d0
625 list_mode, one_over_sigma_tot, sigma_tot_gauss_neumann)
627 sigma_tot_gauss_neumann = 0.d0
630 sigma_curl_gauss = 0.d0
631 sigma_curl_gauss_bdy = 0.d0
632 sigma_curl_gauss_inter_mu = 0.d0
633 j_over_sigma_gauss = 0.d0
634 j_over_sigma_gauss_bdy = 0.d0
635 j_over_sigma_gauss_inter_mu= 0.d0
636 sigma_tot_gauss_neumann = 0.d0
639 IF (
inputs%if_permeability_variable_in_theta)
THEN 641 CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
642 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
643 bloc_size =
SIZE(b_ext,1)/nb_procs+1
646 h_mesh, b_ext, h_pert, nb_procs, bloc_size, m_max_pad, time,temps_par)
649 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
661 rhs_h(:,k) = (4*bn(:,k,i)-bn1(:,k,i))/(2*dt)
664 rhs_phi(:,k) = mu_phi*(4*phin(:,k,i)-phin1(:,k,i))/(2*dt)
668 rhs_h, nl(:,:,i), la_h, la_phi, vb_1, vb_2, b_ext(:,:,i), h_pert(:,:,i), &
669 sigma_curl_gauss(:,:,i), j_over_sigma_gauss(:,:,i))
676 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
680 CALL courant_mu(h_mesh, interface_h_mu, sigma, mu_h_field, time, mode, nl(:,:,i), &
681 la_h, vb_1, vb_2, b_ext(:,:,i), j_over_sigma_gauss_inter_mu(:,:,i), &
682 sigma_curl_gauss_inter_mu(:,:,i))
684 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
690 sigma, mu_h_field, time, mode, nl(:,:,i), stab, la_h, vb_1,vb_2, b_ext(:,:,i), &
691 j_over_sigma_gauss_bdy(:,:,i), sigma_curl_gauss_bdy(:,:,i))
700 CALL surf_int(h_mesh, phi_mesh, pmag_mesh, interface_h_phi, interface_h_mu,
inputs%list_dirichlet_sides_H, &
701 sigma, mu_phi, mu_h_field, time, mode, la_h, la_phi, la_pmag, vb_1, vb_2, &
702 sigma_tot_gauss_neumann(:,:,i), r_fourier, index_fourier)
704 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
710 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 711 CALL periodic_rhs_petsc(h_phi_per%n_bord, h_phi_per%list, h_phi_per%perlist, vb_1, la_mhd)
712 CALL periodic_rhs_petsc(h_phi_per%n_bord, h_phi_per%list, h_phi_per%perlist, vb_2, la_mhd)
719 pmag_global_d(i)%DRL = 0.d0
720 CALL dirichlet_rhs(pmag_mode_global_js_d(i)%DIL-1,pmag_global_d(i)%DRL,vb_1)
721 CALL dirichlet_rhs(pmag_mode_global_js_d(i)%DIL-1,pmag_global_d(i)%DRL,vb_2)
723 IF (
SIZE(phi_js_d)>0)
THEN 728 phi_global_d(i)%DRL(1:n) = phiexact(1,phi_mesh%rr(1:2,phi_js_d), mode, mu_phi, time)
729 IF (
SIZE(phi_global_d(i)%DRL)>n) phi_global_d(i)%DRL(n+1:)=0.d0
730 CALL dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_1)
731 phi_global_d(i)%DRL(1:n) = phiexact(2,phi_mesh%rr(1:2,phi_js_d), mode, mu_phi, time)
732 IF (
SIZE(phi_global_d(i)%DRL)>n) phi_global_d(i)%DRL(n+1:)=0.d0
733 CALL dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_2)
737 phi_global_d(i)%DRL=0.d0
738 CALL dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_1)
739 CALL dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_2)
743 h_global_d(i)%DRL = 0.d0
744 CALL dirichlet_rhs(h_mode_global_js_d(i)%DIL-1,h_global_d(i)%DRL,vb_1)
745 CALL dirichlet_rhs(h_mode_global_js_d(i)%DIL-1,h_global_d(i)%DRL,vb_2)
747 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
752 IF (
inputs%my_par_H_p_phi%verbose .AND. (i==1))
WRITE(*,*)
'start solving' 755 CALL solver(h_p_phi_ksp1(i),vb_1,vx_1,reinit=.false.,
verbose=
inputs%my_par_H_p_phi%verbose)
757 CALL vecghostupdatebegin(vx_1,insert_values,scatter_forward,ierr)
758 CALL vecghostupdateend(vx_1,insert_values,scatter_forward,ierr)
759 IF (h_mesh%me/=0)
THEN 760 CALL extract(vx_1_ghost,1,1,la_mhd,hn_p1(:,1))
761 CALL extract(vx_1_ghost,2,2,la_mhd,hn_p1(:,4))
762 CALL extract(vx_1_ghost,3,3,la_mhd,hn_p1(:,5))
764 IF (phi_mesh%me/=0)
THEN 765 CALL extract(vx_1_ghost,5,5,la_mhd,phin_p1(:,1))
768 CALL solver(h_p_phi_ksp2(i),vb_2,vx_2,reinit=.false.,
verbose=
inputs%my_par_H_p_phi%verbose)
769 CALL vecghostupdatebegin(vx_2,insert_values,scatter_forward,ierr)
770 CALL vecghostupdateend(vx_2,insert_values,scatter_forward,ierr)
771 IF (h_mesh%me/=0)
THEN 772 CALL extract(vx_2_ghost,1,1,la_mhd,hn_p1(:,2))
773 CALL extract(vx_2_ghost,2,2,la_mhd,hn_p1(:,3))
774 CALL extract(vx_2_ghost,3,3,la_mhd,hn_p1(:,6))
776 IF (phi_mesh%me/=0)
THEN 777 CALL extract(vx_2_ghost,5,5,la_mhd,phin_p1(:,2))
780 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
788 IF (h_mesh%me /=0)
THEN 793 IF (phi_mesh%me /=0 )
THEN 806 IF (h_mesh%me /=0)
THEN 807 hn1(:,:,i) = hn(:,:,i)
817 bn1(:,:,i) = bn(:,:,i)
819 bn(:,1,i) = hn_p1(:,1)
820 bn(:,4,i) = hn_p1(:,4)
821 bn(:,5,i) = hn_p1(:,5)
823 bn(:,2,i) = hn_p1(:,2)
824 bn(:,3,i) = hn_p1(:,3)
825 bn(:,6,i) = hn_p1(:,6)
829 IF (phi_mesh%me /= 0)
THEN 830 phin1(:,:,i) = phin(:,:,i)
832 phin(:,1,i) = phin_p1(:,1)
834 phin(:,2,i) = phin_p1(:,2)
836 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
842 IF (h_mesh%me /=0)
THEN 843 IF (
inputs%if_permeability_variable_in_theta)
THEN 844 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
845 bloc_size =
SIZE(bn,1)/nb_procs+1
847 h_mesh, bn, hn, nb_procs, bloc_size, m_max_pad, time,temps_par)
851 hn(:,k,i) = bn(:,k,i)/mu_h_field
858 IF (
inputs%verbose_divergence)
THEN 859 norm =
norm_sf(comm_one_d,
'L2', h_mesh, list_mode, bn)
871 mode, mu_h_field, mu_phi, c_mass, stab, r_fourier, index_fourier, &
872 la_h, la_pmag, la_phi, h_p_phi_mat1, h_p_phi_mat2, sigma_nj_m, sigma)
878 #include "petsc/finclude/petsc.h" 885 INTEGER,
INTENT(IN) :: mode
886 REAL(KIND=8),
INTENT(IN) :: mu_phi, c_mass
887 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
888 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
889 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_w,H_mesh%me),
INTENT(IN) :: sigma_nj_m
890 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN):: sigma
891 REAL(KIND=8),
OPTIONAL :: R_fourier
892 INTEGER,
OPTIONAL :: index_fourier
893 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_ws,phi_mesh%gauss%l_Gs) :: w_cs
894 REAL(KIND=8),
DIMENSION(2, H_mesh%gauss%n_w, phi_mesh%gauss%l_Gs, H_mesh%mes) :: dw_cs
895 INTEGER :: m, l, ms, ls, ni, nj, k, i, j, &
896 n_ws1, n_ws2, n_w2, n_w1, m1, m2, ki, kj,ib,jb, ms1, ms2
897 REAL(KIND=8) :: x, y, hm1, stab_div, stab_colle_H_phi
898 REAL(KIND=8) :: ray, error
899 LOGICAL :: mark=.false.
900 REAL(KIND=8),
DIMENSION(3,H_mesh%gauss%n_w,pmag_mesh%gauss%n_w) :: THpmag
901 REAL(KIND=8),
DIMENSION(pmag_mesh%gauss%n_w,pmag_mesh%gauss%n_w) :: Tpmag
902 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: TH
903 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: HtoB
904 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_w,phi_mesh%gauss%n_w):: TPhi
932 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: Hsij
933 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_w,phi_mesh%gauss%n_w) :: Phisij
934 REAL(KIND=8),
DIMENSION(6,phi_mesh%gauss%n_w,phi_mesh%gauss%n_w) :: Sij
935 REAL(KIND=8),
DIMENSION(6,phi_mesh%gauss%n_w,phi_mesh%gauss%n_w) :: Smuij
956 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: ww_H
958 REAL(KIND=8),
DIMENSION(:,:,:,:),
POINTER :: dw_H
960 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: rj_H
962 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: ww_phi
964 REAL(KIND=8),
DIMENSION(:,:,:,:),
POINTER :: dw_phi
967 REAL(KIND=8),
DIMENSION(2,pmag_mesh%gauss%n_w,H_mesh%gauss%l_G) :: dwp
968 REAL(KIND=8),
DIMENSION(pmag_mesh%gauss%n_w,H_mesh%gauss%l_G) :: wwp
970 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: rj_phi
971 REAL(KIND=8),
DIMENSION(2,phi_mesh%gauss%l_Gs) :: gauss1, gauss2
973 REAL(KIND=8) :: ref, diff, mu_H, c_mu_phi, muhl, &
974 dzmuhl, drmuhl, c_div, hloc, viscolm, xij, eps
976 REAL(KIND=8) :: c_sym=.0d0
979 REAL(KIND=8) :: c_lap
984 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
985 INTEGER ,
DIMENSION(3*H_mesh%gauss%n_w+pmag_mesh%gauss%n_w+ &
phi_mesh%gauss%n_w) :: idxn, jdxn
988 INTEGER :: n_wpmag, n_wH, n_wphi, ix, jx
990 REAL(KIND=8),
DIMENSION(2) :: Dmu_field
991 REAL(KIND=8),
DIMENSION(2) :: gauss_pt
992 INTEGER,
DIMENSION(1) :: gauss_pt_id
993 REAL(KIND=8),
DIMENSION(1) :: dummy_mu_bar
996 REAL(KIND=8) :: sigma_np_gauss
998 mat :: h_p_phi_mat1, h_p_phi_mat2
999 petscerrorcode :: ierr
1000 CALL matzeroentries (h_p_phi_mat1, ierr)
1001 CALL matzeroentries (h_p_phi_mat2, ierr)
1002 CALL matsetoption (h_p_phi_mat1, mat_row_oriented, petsc_false, ierr)
1003 CALL matsetoption (h_p_phi_mat2, mat_row_oriented, petsc_false, ierr)
1007 stab_colle_h_phi = stab(2)
1011 c_mu_phi = c_mass*mu_phi
1013 ww_h => h_mesh%gauss%ww
1014 dw_h => h_mesh%gauss%dw
1015 rj_h => h_mesh%gauss%rj
1016 ww_phi => phi_mesh%gauss%ww
1017 dw_phi => phi_mesh%gauss%dw
1018 rj_phi => phi_mesh%gauss%rj
1020 n_wh = h_mesh%gauss%n_w
1021 n_wpmag = pmag_mesh%gauss%n_w
1022 n_wphi = phi_mesh%gauss%n_w
1029 mesh_id = h_mesh%i_d(m)
1034 DO l = 1, h_mesh%gauss%l_G
1037 hloc = (sqrt(sum(h_mesh%gauss%rj(:,m)))/h_mesh%global_diameter)**(2*
alpha)
1042 IF(
inputs%if_use_fem_integration_for_mu_bar)
THEN 1043 muhl = sum(mu_h_field(h_mesh%jj(:,m))*ww_h(:,l))
1044 drmuhl = sum(mu_h_field(h_mesh%jj(:,m))*dw_h(1,:,l,m))
1045 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m))*dw_h(2,:,l,m))
1047 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jj(:,m))*ww_h(:,l))
1048 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jj(:,m))*ww_h(:,l))
1050 dummy_mu_bar(:) = mu_bar_in_fourier_space(h_mesh,1,1,gauss_pt,gauss_pt_id)
1051 muhl=dummy_mu_bar(1)
1052 dmu_field = grad_mu_bar_in_fourier_space(gauss_pt,gauss_pt_id)
1053 drmuhl =dmu_field(1)
1054 dzmuhl =dmu_field(2)
1056 sigma_np_gauss = sum(sigma_nj_m(:,m)*ww_h(:,l))
1069 c_div = stab_div*hloc/(
inputs%mu_min**2*
inputs%sigma_min)
1074 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
1075 ray = ray + h_mesh%rr(1,i)*ww_h(ni,l)
1078 DO ni = 1, h_mesh%gauss%n_w
1079 DO nj = 1, h_mesh%gauss%n_w
1082 th(1,ni,nj) = th(1,ni,nj) + rj_h(l,m) * ray* ( &
1083 c_mass*ww_h(ni,l)*ww_h(nj,l) &
1084 + (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*muhl) &
1087 + c_div*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) + ww_h(ni,l)*0) &
1088 *(ww_h(nj,l)/ray+dw_h(1,nj,l,m))) &
1089 + rj_h(l,m) * ray* ( &
1090 (dw_h(2,ni,l,m)*dzmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
1093 th(2,ni,nj) = th(2,ni,nj)+ rj_h(l,m) * ray* ( &
1094 mode/ray**2 * ww_h(ni,l)*(ww_h(nj,l)+ray*dw_h(1,nj,l,m))/(sigma_np_gauss*muhl) &
1097 + c_div*mode/ray*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) + ww_h(ni,l)*0)*ww_h(nj,l)) &
1098 + rj_h(l,m) * ray* ( &
1099 ((mode/ray)*ww_h(ni,l)*drmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
1101 htob(2,ni,nj) = htob(2,ni,nj)+ rj_h(l,m) * ray* ( &
1102 - mode/ray**2 * ww_h(ni,l)*(ww_h(nj,l)+ray*dw_h(1,nj,l,m))/(sigma_np_gauss*muhl) &
1105 - c_div*mode/ray*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) + ww_h(ni,l)*0)*ww_h(nj,l)) &
1106 - rj_h(l,m) * ray* ( &
1107 ((mode/ray)*ww_h(ni,l)*drmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
1109 th(3,ni,nj) = th(3,ni,nj)+ rj_h(l,m) * ray* ( &
1110 - dw_h(2,ni,l,m)*dw_h(1,nj,l,m)/(sigma_np_gauss*muhl) &
1113 + c_div*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) + ww_h(ni,l)*0)*(dw_h(2,nj,l,m))) &
1114 + rj_h(l,m) * ray* ( &
1115 - dw_h(2,ni,l,m)*drmuhl*ww_h(nj,l)*(-1/(sigma_np_gauss*muhl**2)))
1117 th(4,ni,nj) = th(4,ni,nj)+ rj_h(l,m) * ray* ( &
1118 mode/ray**2 * ww_h(nj,l)*(ww_h(ni,l)+ray*dw_h(1,ni,l,m))/(sigma_np_gauss*muhl) &
1121 + c_div*(mode/ray)*ww_h(ni,l)*muhl*(ww_h(nj,l)/ray+dw_h(1,nj,l,m)))
1123 htob(4,ni,nj) = htob(4,ni,nj)+ rj_h(l,m) * ray* ( &
1124 - mode/ray**2 * ww_h(nj,l)*(ww_h(ni,l)+ray*dw_h(1,ni,l,m))/(sigma_np_gauss*muhl) &
1127 - c_div*(mode/ray)*ww_h(ni,l)*muhl*(ww_h(nj,l)/ray+dw_h(1,nj,l,m)))
1129 th(5,ni,nj) = th(5,ni,nj) + rj_h(l,m) * ray* ( &
1130 c_mass*ww_h(ni,l)*ww_h(nj,l) &
1131 + (dw_h(2,ni,l,m)*dw_h(2,nj,l,m) &
1132 + 1/ray**2*(ww_h(ni,l)+ray*dw_h(1,ni,l,m))*(ww_h(nj,l)+ray*dw_h(1,nj,l,m)))/(sigma_np_gauss*muhl) &
1135 +c_div*muhl*mode**2/ray**2*ww_h(ni,l)*ww_h(nj,l)) &
1136 + rj_h(l,m)*ray*((dw_h(2,ni,l,m)*dzmuhl*ww_h(nj,l) &
1137 + (ww_h(ni,l)/ray+dw_h(1,ni,l,m))*drmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
1139 th(6,ni,nj) = th(6,ni,nj) + rj_h(l,m) * ray* (&
1140 + mode/ray*dw_h(2,ni,l,m)*ww_h(nj,l)/(sigma_np_gauss*muhl) &
1143 +c_div*mode/ray*muhl*ww_h(ni,l)*(dw_h(2,nj,l,m)))
1145 htob(6,ni,nj) = htob(6,ni,nj) + rj_h(l,m) * ray* (&
1146 - mode/ray*dw_h(2,ni,l,m)*ww_h(nj,l)/(sigma_np_gauss*muhl) &
1149 - c_div*mode/ray*muhl*ww_h(ni,l)*(dw_h(2,nj,l,m)))
1151 th(7,ni,nj) = th(7,ni,nj)+ rj_h(l,m) * ray* ( &
1152 - dw_h(1,ni,l,m)*dw_h(2,nj,l,m)/(sigma_np_gauss*muhl) &
1155 + c_div*(muhl*dw_h(2,ni,l,m)+0*ww_h(ni,l))*(ww_h(nj,l)/ray+dw_h(1,nj,l,m))) &
1156 + rj_h(l,m)*ray*((-dw_h(1,ni,l,m)*dzmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
1158 th(8,ni,nj) = th(8,ni,nj) + rj_h(l,m) * ray* (&
1159 + (mode/ray)*ww_h(ni,l)*dw_h(2,nj,l,m)/(sigma_np_gauss*muhl) &
1162 + c_div*mode/ray*muhl*ww_h(nj,l)*(dw_h(2,ni,l,m))) &
1163 + rj_h(l,m)*ray*(((mode/ray)*ww_h(ni,l)*dzmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
1167 htob(8,ni,nj) = htob(8,ni,nj) + rj_h(l,m) * ray* (&
1168 - mode/ray*dw_h(2,nj,l,m)*ww_h(ni,l)/(sigma_np_gauss*muhl) &
1171 - c_div*mode/ray*muhl*ww_h(nj,l)*(dw_h(2,ni,l,m))) &
1172 - rj_h(l,m)*ray*(((mode/ray)*ww_h(ni,l)*dzmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
1176 th(9,ni,nj) = th(9,ni,nj) + rj_h(l,m) * ray* ( &
1177 c_mass*ww_h(ni,l)*ww_h(nj,l) &
1178 + (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*muhl) &
1181 + c_div*(muhl*dw_h(2,ni,l,m) + ww_h(ni,l)*0)*(dw_h(2,nj,l,m))) &
1182 + rj_h(l,m)*ray*((dw_h(1,ni,l,m)*drmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
1290 i = h_mesh%jj(ni, m)
1291 ib = la_h%loc_to_glob(ki,i)
1296 j = h_mesh%jj(nj, m)
1297 jb = la_h%loc_to_glob(kj,j)
1301 IF ((ki == 1) .AND. (kj == 1))
THEN 1302 mat_loc1(ix,jx) = th(1,ni,nj)
1303 mat_loc2(ix,jx) = th(1,ni,nj)
1304 ELSEIF ((ki == 1) .AND. (kj == 2))
THEN 1305 mat_loc1(ix,jx) = th(2,ni,nj)
1306 mat_loc2(ix,jx) = htob(2,ni,nj)
1307 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 1308 mat_loc1(ix,jx) = th(3,ni,nj)
1309 mat_loc2(ix,jx) = th(3,ni,nj)
1310 ELSEIF ((ki == 2) .AND. (kj == 1))
THEN 1311 mat_loc1(ix,jx) = th(4,ni,nj)
1312 mat_loc2(ix,jx) = htob(4,ni,nj)
1313 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 1314 mat_loc1(ix,jx) = th(5,ni,nj)
1315 mat_loc2(ix,jx) = th(5,ni,nj)
1316 ELSEIF ((ki == 2) .AND. (kj == 3))
THEN 1317 mat_loc1(ix,jx) = th(6,ni,nj)
1318 mat_loc2(ix,jx) = htob(6,ni,nj)
1319 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 1320 mat_loc1(ix,jx) = th(7,ni,nj)
1321 mat_loc2(ix,jx) = th(7,ni,nj)
1322 ELSEIF ((ki == 3) .AND. (kj == 2))
THEN 1323 mat_loc1(ix,jx) = th(8,ni,nj)
1324 mat_loc2(ix,jx) = htob(8,ni,nj)
1325 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 1326 mat_loc1(ix,jx) = th(9,ni,nj)
1327 mat_loc2(ix,jx) = th(9,ni,nj)
1334 CALL matsetvalues(h_p_phi_mat1, 3*n_wh, idxn(1:3*n_wh), 3*n_wh, jdxn(1:3*n_wh), &
1335 mat_loc1(1:3*n_wh,1:3*n_wh), add_values, ierr)
1336 CALL matsetvalues(h_p_phi_mat2, 3*n_wh, idxn(1:3*n_wh), 3*n_wh, jdxn(1:3*n_wh), &
1337 mat_loc2(1:3*n_wh,1:3*n_wh), add_values, ierr)
1341 DO m = 1, pmag_mesh%me
1346 hloc = stab_div*(sqrt(sum(pmag_mesh%gauss%rj(:,m)))/pmag_mesh%global_diameter)**(2*(1-
alpha))
1347 mesh_id = pmag_mesh%i_d(m)
1350 DO l = 1, pmag_mesh%gauss%l_G
1353 IF(
inputs%if_use_fem_integration_for_mu_bar)
THEN 1354 muhl = sum(mu_h_field(h_mesh%jj(1:3,m))*pmag_mesh%gauss%ww(:,l))
1356 gauss_pt(1)=sum(pmag_mesh%rr(1,pmag_mesh%jj(:,m))*pmag_mesh%gauss%ww(:,l))
1357 gauss_pt(2)=sum(pmag_mesh%rr(2,pmag_mesh%jj(:,m))*pmag_mesh%gauss%ww(:,l))
1359 dummy_mu_bar(:) = mu_bar_in_fourier_space(pmag_mesh,1,1,gauss_pt,gauss_pt_id)
1360 muhl=dummy_mu_bar(1)
1362 sigma_np_gauss = sum(sigma_nj_m(1:3,m)*pmag_mesh%gauss%ww(:,l))
1373 DO ni = 1, pmag_mesh%gauss%n_w
1374 i = pmag_mesh%jj(ni,m)
1375 ray = ray + pmag_mesh%rr(1,i)*pmag_mesh%gauss%ww(ni,l)
1382 viscolm = (pmag_mesh%global_diameter)**2*
inputs%mu_min**2*
inputs%sigma_min*hloc*pmag_mesh%gauss%rj(l,m)
1384 DO nj = 1, pmag_mesh%gauss%n_w
1385 j = pmag_mesh%jj(nj, m)
1386 DO ni = 1, pmag_mesh%gauss%n_w
1387 i = pmag_mesh%jj(ni, m)
1391 xij = xij + pmag_mesh%gauss%dw(k,nj,l,m) * pmag_mesh%gauss%dw(k,ni,l,m)
1394 tpmag(ni,nj) = tpmag(ni,nj) + ray * viscolm* xij &
1395 + viscolm*mode**2*pmag_mesh%gauss%ww(ni,l)*pmag_mesh%gauss%ww(nj,l)/ray
1400 DO ni = 1, pmag_mesh%gauss%n_w
1401 i = pmag_mesh%jj(ni, m)
1402 ib = la_pmag%loc_to_glob(1,i)
1404 DO nj = 1, pmag_mesh%gauss%n_w
1405 j = pmag_mesh%jj(nj, m)
1406 jb = la_pmag%loc_to_glob(1,j)
1410 CALL matsetvalues(h_p_phi_mat1, n_wpmag, idxn(1:n_wpmag), n_wpmag, jdxn(1:n_wpmag), &
1411 tpmag(1:n_wpmag,1:n_wpmag), add_values, ierr)
1412 CALL matsetvalues(h_p_phi_mat2, n_wpmag, idxn(1:n_wpmag), n_wpmag, jdxn(1:n_wpmag), &
1413 tpmag(1:n_wpmag,1:n_wpmag), add_values, ierr)
1418 DO m = 1, pmag_mesh%me
1419 mesh_id = h_mesh%i_d(m)
1420 IF (h_mesh%gauss%n_w==3)
THEN 1421 dwp=h_mesh%gauss%dw(:,:,:,m)
1424 dwp(:,1,:) = h_mesh%gauss%dw(:,1,:,m) + 0.5d0*(h_mesh%gauss%dw(:,5,:,m)+h_mesh%gauss%dw(:,6,:,m))
1425 dwp(:,2,:) = h_mesh%gauss%dw(:,2,:,m) + 0.5d0*(h_mesh%gauss%dw(:,6,:,m)+h_mesh%gauss%dw(:,4,:,m))
1426 dwp(:,3,:) = h_mesh%gauss%dw(:,3,:,m) + 0.5d0*(h_mesh%gauss%dw(:,4,:,m)+h_mesh%gauss%dw(:,5,:,m))
1427 wwp(1,:) = h_mesh%gauss%ww(1,:) + 0.5d0*(h_mesh%gauss%ww(5,:)+h_mesh%gauss%ww(6,:))
1428 wwp(2,:) = h_mesh%gauss%ww(2,:) + 0.5d0*(h_mesh%gauss%ww(6,:)+h_mesh%gauss%ww(4,:))
1429 wwp(3,:) = h_mesh%gauss%ww(3,:) + 0.5d0*(h_mesh%gauss%ww(4,:)+h_mesh%gauss%ww(5,:))
1433 DO l = 1, h_mesh%gauss%l_G
1435 DO ni = 1, h_mesh%gauss%n_w
1437 ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
1441 IF(
inputs%if_use_fem_integration_for_mu_bar)
THEN 1442 muhl = stab_div*ray*h_mesh%gauss%rj(l,m)*sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
1445 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jj(:,m))*ww_h(:,l))
1446 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jj(:,m))*ww_h(:,l))
1447 gauss_pt_id(1)=mesh_id
1448 dummy_mu_bar(:) = mu_bar_in_fourier_space(h_mesh,1,1,gauss_pt,gauss_pt_id)
1449 muhl=dummy_mu_bar(1)
1450 muhl = stab_div*ray*h_mesh%gauss%rj(l,m)*muhl
1452 DO nj = 1, pmag_mesh%gauss%n_w
1453 j = pmag_mesh%jj(nj, m)
1454 DO ni = 1, h_mesh%gauss%n_w
1455 i = h_mesh%jj(ni, m)
1456 thpmag(1,ni,nj) = thpmag(1,ni,nj) + muhl*dwp(1,nj,l)*h_mesh%gauss%ww(ni,l)
1457 thpmag(2,ni,nj) = thpmag(2,ni,nj) - muhl*mode*wwp(nj,l)*h_mesh%gauss%ww(ni,l)/ray
1458 thpmag(3,ni,nj) = thpmag(3,ni,nj) + muhl*dwp(2,nj,l)*h_mesh%gauss%ww(ni,l)
1468 i = h_mesh%jj(ni, m)
1475 ib = la_h%loc_to_glob(k,i)
1476 ix = (k-1)*n_wh + ni
1479 j = pmag_mesh%jj(nj, m)
1480 jb = la_pmag%loc_to_glob(1,j)
1483 mat_loc1(ix,jx) = thpmag(k,ni,nj)
1484 mat_loc2(ix,jx) = eps*thpmag(k,ni,nj)
1489 CALL matsetvalues(h_p_phi_mat1, 3*n_wh, idxn(1:3*n_wh), n_wpmag, jdxn(1:n_wpmag), &
1490 mat_loc1(1:3*n_wh,1:n_wpmag), add_values, ierr)
1491 CALL matsetvalues(h_p_phi_mat2, 3*n_wh, idxn(1:3*n_wh), n_wpmag, jdxn(1:n_wpmag), &
1492 mat_loc2(1:3*n_wh,1:n_wpmag), add_values, ierr)
1496 DO l = 1, h_mesh%gauss%l_G
1498 DO ni = 1, h_mesh%gauss%n_w
1500 ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
1502 x = stab_div*ray*h_mesh%gauss%rj(l,m)
1503 DO nj = 1, h_mesh%gauss%n_w
1504 j = h_mesh%jj(nj, m)
1505 DO ni = 1, pmag_mesh%gauss%n_w
1506 i = pmag_mesh%jj(ni, m)
1507 thpmag(1,nj,ni) = thpmag(1,nj,ni) - x*dwp(1,ni,l)*h_mesh%gauss%ww(nj,l)
1508 thpmag(2,nj,ni) = thpmag(2,nj,ni) + x*mode*wwp(ni,l)*h_mesh%gauss%ww(nj,l)/ray
1509 thpmag(3,nj,ni) = thpmag(3,nj,ni) - x*dwp(2,ni,l)*h_mesh%gauss%ww(nj,l)
1517 i = pmag_mesh%jj(ni, m)
1518 ib = la_pmag%loc_to_glob(1,i)
1528 j = h_mesh%jj(nj, m)
1529 jb = la_h%loc_to_glob(k,j)
1530 jx = (k-1)*n_wh + nj
1532 mat_loc1(ix,jx) = thpmag(k,nj,ni)
1533 mat_loc2(ix,jx) = eps*thpmag(k,nj,ni)
1537 CALL matsetvalues(h_p_phi_mat1, n_wpmag, idxn(1:n_wpmag), 3*n_wh, jdxn(1:3*n_wh), &
1538 mat_loc1(1:n_wpmag,1:3*n_wh), add_values, ierr)
1539 CALL matsetvalues(h_p_phi_mat2, n_wpmag, idxn(1:n_wpmag), 3*n_wh, jdxn(1:3*n_wh), &
1540 mat_loc2(1:n_wpmag,1:3*n_wh), add_values, ierr)
1547 DO m = 1,phi_mesh%me
1551 DO l = 1, phi_mesh%gauss%l_G
1555 DO ni = 1, phi_mesh%gauss%n_w; i = phi_mesh%jj(ni,m)
1556 ray = ray + phi_mesh%rr(1,i)*ww_phi(ni,l)
1559 DO ni = 1, phi_mesh%gauss%n_w
1560 DO nj = 1, phi_mesh%gauss%n_w
1568 tphi(ni,nj) = tphi(ni,nj) + rj_phi(l,m) * ray* (c_mass+c_lap)*mu_phi &
1569 *(dw_phi(1,ni,l,m)*dw_phi(1,nj,l,m)+dw_phi(2,ni,l,m)*dw_phi(2,nj,l,m) &
1570 +mode**2/ray**2*ww_phi(ni,l)*ww_phi(nj,l))
1577 DO ni = 1, phi_mesh%gauss%n_w
1578 i = phi_mesh%jj(ni, m)
1579 ib = la_phi%loc_to_glob(1,i)
1581 DO nj = 1, phi_mesh%gauss%n_w
1582 j = phi_mesh%jj(nj, m)
1583 jb = la_phi%loc_to_glob(1,j)
1587 CALL matsetvalues(h_p_phi_mat1, n_wphi, idxn(1:n_wphi), n_wphi, jdxn(1:n_wphi), &
1588 tphi(1:n_wphi,1:n_wphi), add_values, ierr)
1589 CALL matsetvalues(h_p_phi_mat2, n_wphi, idxn(1:n_wphi), n_wphi, jdxn(1:n_wphi), &
1590 tphi(1:n_wphi,1:n_wphi), add_values, ierr)
1598 CALL gauss(phi_mesh)
1599 n_ws1 = h_mesh%gauss%n_ws
1600 n_ws2 = phi_mesh%gauss%n_ws
1601 n_w1 = h_mesh%gauss%n_w
1602 n_w2 = phi_mesh%gauss%n_w
1604 IF (h_mesh%gauss%n_ws ==
n_ws)
THEN 1606 DO ms = 1, interface_h_phi%mes
1608 ms2 = interface_h_phi%mesh2(ms)
1609 m2 = phi_mesh%neighs(ms2)
1610 ms1 = interface_h_phi%mesh1(ms)
1611 m1 = h_mesh%neighs(ms1)
1613 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
1614 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - phi_mesh%rr(:,phi_mesh%jjs(1,ms2)))**2)
1615 IF (diff/ref .LT. 1.d-10)
THEN 1619 w_cs(1,ls)=
wws(2,ls)
1620 w_cs(2,ls)=
wws(1,ls)
1621 w_cs(3,ls)=
wws(3,ls)
1622 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 1627 gauss2(1,ls) = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))*phi_mesh%gauss%wws(:,ls))
1628 gauss2(2,ls) = sum(phi_mesh%rr(2,phi_mesh%jjs(:,ms2))*phi_mesh%gauss%wws(:,ls))
1629 gauss1(1,ls) = sum( h_mesh%rr(1, h_mesh%jjs(:,ms1))* h_mesh%gauss%wws(:,ls))
1630 gauss1(2,ls) = sum( h_mesh%rr(2, h_mesh%jjs(:,ms1))* h_mesh%gauss%wws(:,ls))
1634 ref = sqrt(1.d-8+sum(gauss2(:,ls2)**2))
1637 diff = sqrt(sum((gauss2(:,ls2)-gauss1(:,ls1))**2))
1638 IF (diff .LT. 1.d-10)
THEN 1639 dw_cs(:,:,ls2,ms1) = h_mesh%gauss%dw_s(:,:,ls1,ms1)
1644 IF (.NOT.mark)
WRITE(*,*)
' BUG ' 1650 DO ms = 1, interface_h_phi%mes
1652 ms2 = interface_h_phi%mesh2(ms)
1653 m2 = phi_mesh%neighs(ms2)
1654 ms1 = interface_h_phi%mesh1(ms)
1655 m1 = h_mesh%neighs(ms1)
1657 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
1658 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - phi_mesh%rr(:,phi_mesh%jjs(1,ms2)))**2)
1659 IF (diff/ref .LT. 1.d-10)
THEN 1661 w_cs(1,ls)=
wws(1,ls)+0.5*
wws(3,ls)
1662 w_cs(2,ls)=
wws(2,ls)+0.5*
wws(3,ls)
1667 w_cs(1,ls)=
wws(2,ls)+0.5*
wws(3,ls)
1668 w_cs(2,ls)=
wws(1,ls)+0.5*
wws(3,ls)
1670 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 1675 dw_cs(1,:,ls,ms1) = h_mesh%gauss%dw(1,:,1,m1)
1676 dw_cs(2,:,ls,ms1) = h_mesh%gauss%dw(2,:,1,m1)
1683 DO ms = 1, interface_h_phi%mes
1685 ms2 = interface_h_phi%mesh2(ms)
1686 ms1 = interface_h_phi%mesh1(ms)
1687 m2 = phi_mesh%neighs(ms2)
1688 m1 = h_mesh%neighs(ms1)
1689 mu_h = sum(mu_h_field(h_mesh%jj(:,m1)))/h_mesh%gauss%n_w
1694 hm1 = stab_colle_h_phi/(sum(
rjs(:,ms2))*
inputs%sigma_min)
1707 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1708 x = hm1*
rjs(ls,ms2)*ray
1710 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1716 y = x * w_cs(ni,ls)*w_cs(nj,ls)/muhl
1718 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(
rnorms(2,ls,ms2)**2)
1719 hsij(4,ni,nj) = hsij(4,ni,nj) - y*
rnorms(1,ls,ms2)*
rnorms(2,ls,ms2)
1720 hsij(5,ni,nj) = hsij(5,ni,nj) + y
1721 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(
rnorms(1,ls,ms2)**2)
1732 i = interface_h_phi%jjs1(ni,ms)
1733 ib = la_h%loc_to_glob(ki,i)
1734 ix = (ki-1)*n_ws1+ni
1738 j = interface_h_phi%jjs1(nj,ms)
1739 jb = la_h%loc_to_glob(kj,j)
1740 jx = (kj-1)*n_ws1+nj
1742 IF ((ki == 1) .AND. (kj == 1))
THEN 1743 mat_loc1(ix,jx) = hsij(1,ni,nj)
1744 mat_loc2(ix,jx) = hsij(1,ni,nj)
1745 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 1746 mat_loc1(ix,jx) = hsij(4,ni,nj)
1747 mat_loc2(ix,jx) = hsij(4,ni,nj)
1748 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 1749 mat_loc1(ix,jx) = hsij(4,nj,ni)
1750 mat_loc2(ix,jx) = hsij(4,nj,ni)
1751 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 1752 mat_loc1(ix,jx) = hsij(5,ni,nj)
1753 mat_loc2(ix,jx) = hsij(5,ni,nj)
1754 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 1755 mat_loc1(ix,jx) = hsij(6,ni,nj)
1756 mat_loc2(ix,jx) = hsij(6,ni,nj)
1762 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1763 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1764 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1765 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1771 DO ls = 1, phi_mesh%gauss%l_Gs
1774 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1775 x =
rjs(ls,ms2)*ray/sigma(m1)
1780 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1781 drmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(1,:,ls,ms1))
1782 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(2,:,ls,ms1))
1788 y = x*w_cs(ni,ls)*w_cs(nj,ls)
1795 hsij(2,ni,nj) = hsij(2,ni,nj) + y * (-mode/ray)*(-
rnorms(1,ls,ms2))/muhl
1796 hsij(3,ni,nj) = hsij(3,ni,nj) + y * mode/ray *(-
rnorms(1,ls,ms2))/muhl
1797 hsij(5,ni,nj) = hsij(5,ni,nj) + y * (-1/ray) *(-
rnorms(1,ls,ms2))/muhl
1798 hsij(8,ni,nj) = hsij(8,ni,nj) + y * (-mode/ray)*(-
rnorms(2,ls,ms2))/muhl
1799 hsij(9,ni,nj) = hsij(9,ni,nj) + y * mode/ray *(-
rnorms(2,ls,ms2))/muhl
1801 hsij(1,ni,nj) = hsij(1,ni,nj) - y*(-
rnorms(2,ls,ms2))*(-(dzmuhl/muhl**2))
1802 hsij(4,ni,nj) = hsij(4,ni,nj) - y*(-
rnorms(2,ls,ms2))*( (drmuhl/muhl**2))
1803 hsij(5,ni,nj) = hsij(5,ni,nj) &
1804 + y*(-
rnorms(1,ls,ms2)*drmuhl-
rnorms(2,ls,ms2)*dzmuhl)/muhl**2
1805 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(-
rnorms(1,ls,ms2))*( (drmuhl/muhl**2))
1806 hsij(7,ni,nj) = hsij(7,ni,nj) + y*(-
rnorms(1,ls,ms2))*(-(dzmuhl/muhl**2))
1821 i = interface_h_phi%jjs1(ni,ms)
1822 ib = la_h%loc_to_glob(ki,i)
1823 ix = (ki-1)*n_ws1 + ni
1827 j = interface_h_phi%jjs1(nj,ms)
1828 jb = la_h%loc_to_glob(kj,j)
1829 jx = (kj-1)*n_ws1 + nj
1832 IF ( (ki == 1) .AND. (kj == 1))
THEN 1833 mat_loc1(ix,jx) = hsij(1,ni,nj)
1834 mat_loc2(ix,jx) = hsij(1,ni,nj)
1835 ELSE IF ( (ki == 1) .AND. (kj == 3))
THEN 1836 mat_loc1(ix,jx) = hsij(4,ni,nj)
1837 mat_loc2(ix,jx) = hsij(4,ni,nj)
1839 ELSE IF ( (ki == 2) .AND. (kj == 1))
THEN 1840 mat_loc1(ix,jx) = hsij(2,ni,nj)
1841 mat_loc2(ix,jx) = hsij(3,ni,nj)
1842 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 1843 mat_loc1(ix,jx) = hsij(5,ni,nj)
1844 mat_loc2(ix,jx) = hsij(5,ni,nj)
1845 ELSEIF ( (ki == 2) .AND. (kj == 3))
THEN 1846 mat_loc1(ix,jx) = hsij(8,ni,nj)
1847 mat_loc2(ix,jx) = hsij(9,ni,nj)
1849 ELSE IF ( (ki == 3) .AND. (kj == 1))
THEN 1850 mat_loc1(ix,jx) = hsij(7,ni,nj)
1851 mat_loc2(ix,jx) = hsij(7,ni,nj)
1852 ELSE IF ( (ki == 3) .AND. (kj == 3))
THEN 1853 mat_loc1(ix,jx) = hsij(6,ni,nj)
1854 mat_loc2(ix,jx) = hsij(6,ni,nj)
1861 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1862 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1863 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1864 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1873 i = interface_h_phi%jjs1(ni,ms)
1874 ib = la_h%loc_to_glob(ki,i)
1875 ix = (ki-1)*n_ws1 + ni
1879 j = interface_h_phi%jjs1(nj,ms)
1880 jb = la_h%loc_to_glob(kj,j)
1881 jx = (kj-1)*n_ws1 + nj
1883 IF ( (kj == 2) .AND. (ki == 1))
THEN 1884 mat_loc1(ix,jx) = hsij(2,nj,ni)
1885 mat_loc2(ix,jx) = hsij(3,nj,ni)
1886 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN 1887 mat_loc1(ix,jx) = hsij(5,nj,ni)
1888 mat_loc2(ix,jx) = hsij(5,nj,ni)
1889 ELSEIF ( (kj == 2) .AND. (ki == 3))
THEN 1890 mat_loc1(ix,jx) = hsij(8,nj,ni)
1891 mat_loc2(ix,jx) = hsij(9,nj,ni)
1898 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1899 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1900 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1901 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1904 DO ls = 1, phi_mesh%gauss%l_Gs
1907 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1908 x =
rjs(ls,ms2)*ray /sigma(m1)
1913 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1926 hsij(1,ni,nj) = hsij(1,ni,nj) - y*(-
rnorms(2,ls,ms2))*(dw_cs(2,nj,ls,ms1)/muhl)
1927 hsij(4,ni,nj) = hsij(4,ni,nj) - y*(-
rnorms(2,ls,ms2))*(-dw_cs(1,nj,ls,ms1)/muhl)
1928 hsij(5,ni,nj) = hsij(5,ni,nj) &
1929 + y*(-
rnorms(2,ls,ms2))*(-dw_cs(2,nj,ls,ms1)/muhl) &
1930 - y*(-
rnorms(1,ls,ms2))*( dw_cs(1,nj,ls,ms1)/muhl)
1931 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(-
rnorms(1,ls,ms2))*(-dw_cs(1,nj,ls,ms1)/muhl)
1932 hsij(7,ni,nj) = hsij(7,ni,nj) + y*(-
rnorms(1,ls,ms2))*(dw_cs(2,nj,ls,ms1)/muhl)
1946 i = interface_h_phi%jjs1(ni,ms)
1947 ib = la_h%loc_to_glob(ki,i)
1948 ix = (ki-1)*n_ws1 + ni
1952 j = h_mesh%jj(nj,m1)
1953 jb = la_h%loc_to_glob(kj,j)
1954 jx = (kj-1)*n_w1 + nj
1956 IF ((ki == 1) .AND. (kj == 1))
THEN 1957 mat_loc1(ix,jx) = hsij(1,ni,nj)
1958 mat_loc2(ix,jx) = hsij(1,ni,nj)
1959 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 1960 mat_loc1(ix,jx) = hsij(4,ni,nj)
1961 mat_loc2(ix,jx) = hsij(4,ni,nj)
1962 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 1963 mat_loc1(ix,jx) = hsij(5,ni,nj)
1964 mat_loc2(ix,jx) = hsij(5,ni,nj)
1965 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 1966 mat_loc1(ix,jx) = hsij(6,ni,nj)
1967 mat_loc2(ix,jx) = hsij(6,ni,nj)
1968 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 1969 mat_loc1(ix,jx) = hsij(7,ni,nj)
1970 mat_loc2(ix,jx) = hsij(7,ni,nj)
1977 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
1978 mat_loc1(1:3*n_ws1,1:3*n_w1), add_values, ierr)
1979 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
1980 mat_loc2(1:3*n_ws1,1:3*n_w1), add_values, ierr)
1988 i = h_mesh%jj(ni,m1)
1989 ib = la_h%loc_to_glob(ki,i)
1990 ix = (ki-1)*n_w1 + ni
1994 j = interface_h_phi%jjs1(nj,ms)
1995 jb = la_h%loc_to_glob(kj,j)
1996 jx = (kj-1)*n_ws1 + nj
1998 IF ((kj == 1) .AND. (ki == 1))
THEN 1999 mat_loc1(ix,jx) = hsij(1,nj,ni)
2000 mat_loc2(ix,jx) = hsij(1,nj,ni)
2001 ELSEIF ((kj == 1) .AND. (ki == 3))
THEN 2002 mat_loc1(ix,jx) = hsij(4,nj,ni)
2003 mat_loc2(ix,jx) = hsij(4,nj,ni)
2004 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN 2005 mat_loc1(ix,jx) = hsij(5,nj,ni)
2006 mat_loc2(ix,jx) = hsij(5,nj,ni)
2007 ELSEIF ((kj == 3) .AND. (ki == 3))
THEN 2008 mat_loc1(ix,jx) = hsij(6,nj,ni)
2009 mat_loc2(ix,jx) = hsij(6,nj,ni)
2010 ELSEIF ((kj == 3) .AND. (ki == 1))
THEN 2011 mat_loc1(ix,jx) = hsij(7,nj,ni)
2012 mat_loc2(ix,jx) = hsij(7,nj,ni)
2018 CALL matsetvalues(h_p_phi_mat1, 3*n_w1, idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
2019 mat_loc1(1:3*n_w1,1:3*n_ws1), add_values, ierr)
2020 CALL matsetvalues(h_p_phi_mat2, 3*n_w1, idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
2021 mat_loc2(1:3*n_w1,1:3*n_ws1), add_values, ierr)
2034 DO ls = 1, phi_mesh%gauss%l_Gs
2037 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2038 x = hm1*
rjs(ls,ms2)*ray
2043 phisij(ni,nj) = phisij(ni,nj) + x*mode**2/ray**2*
wws(ni,ls)*
wws(nj,ls)
2054 i = interface_h_phi%jjs2(ni,ms)
2055 ib = la_phi%loc_to_glob(1,i)
2058 j = interface_h_phi%jjs2(nj,ms)
2059 jb = la_phi%loc_to_glob(1,j)
2063 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
2064 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
2065 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
2066 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
2072 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2073 x = hm1*
rjs(ls,ms2)*ray
2078 phisij(ni,nj) = phisij(ni,nj) + x*( &
2079 (
dw_s(2,ni,ls,ms2)*
rnorms(1,ls,ms2) -
dw_s(1,ni,ls,ms2)*
rnorms(2,ls,ms2))* &
2080 (
dw_s(2,nj,ls,ms2)*
rnorms(1,ls,ms2) -
dw_s(1,nj,ls,ms2)*
rnorms(2,ls,ms2)))
2090 i = phi_mesh%jj(ni, m2)
2091 ib = la_phi%loc_to_glob(1,i)
2094 j = phi_mesh%jj(nj, m2)
2095 jb = la_phi%loc_to_glob(1,j)
2099 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), n_w2, jdxn(1:n_w2), &
2100 phisij(1:n_w2,1:n_w2), add_values, ierr)
2101 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), n_w2, jdxn(1:n_w2), &
2102 phisij(1:n_w2,1:n_w2), add_values, ierr)
2116 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2117 x = hm1*
rjs(ls,ms2)*ray
2119 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2124 sij(3,ni,nj) = sij(3,ni,nj) + x*(mode/ray)*w_cs(ni,ls)*
wws(nj,ls)
2125 smuij(3,ni,nj) = smuij(3,ni,nj) + x*(mode/ray)*w_cs(ni,ls)*
wws(nj,ls)/muhl
2129 sij(4,:,:) = -sij(3,:,:)
2131 smuij(4,:,:) = -smuij(3,:,:)
2141 i = interface_h_phi%jjs1(ni,ms)
2142 ib = la_h%loc_to_glob(ki,i)
2145 j = interface_h_phi%jjs2(nj,ms)
2146 jb = la_phi%loc_to_glob(1,j)
2150 CALL matsetvalues(h_p_phi_mat1, n_ws1, idxn(1:n_ws1), n_ws2, jdxn(1:n_ws2), &
2151 sij(3,1:n_ws1,1:n_ws2), add_values, ierr)
2152 CALL matsetvalues(h_p_phi_mat2, n_ws1, idxn(1:n_ws1), n_ws2, jdxn(1:n_ws2), &
2153 sij(4,1:n_ws1,1:n_ws2), add_values, ierr)
2162 i = interface_h_phi%jjs2(ni,ms)
2163 ib = la_phi%loc_to_glob(1,i)
2166 j = interface_h_phi%jjs1(nj,ms)
2167 jb = la_h%loc_to_glob(kj,j)
2172 mat_loc1(ni,nj) = smuij(3,nj,ni)
2173 mat_loc2(ni,nj) = smuij(4,nj,ni)
2178 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws1, jdxn(1:n_ws1), &
2179 mat_loc1(1:n_ws2,1:n_ws1), add_values, ierr)
2180 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws1, jdxn(1:n_ws1), &
2181 mat_loc2(1:n_ws2,1:n_ws1), add_values, ierr)
2191 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2192 x = hm1*
rjs(ls,ms2)*ray
2195 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2201 sij(1,ni,nj) = sij(1,ni,nj) + &
2202 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))
2204 smuij(1,ni,nj) = smuij(1,ni,nj) + &
2205 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))/muhl
2207 sij(5,ni,nj) = sij(5,ni,nj) + &
2208 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))
2210 smuij(5,ni,nj) = smuij(5,ni,nj) + &
2211 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))/muhl
2225 i = interface_h_phi%jjs1(ni,ms)
2226 ib = la_h%loc_to_glob(ki,i)
2227 ix = (ki-1)*n_ws1 + ni
2230 j = phi_mesh%jj(nj,m2)
2231 jb = la_phi%loc_to_glob(1,j)
2235 mat_loc1(ix,jx) = sij(1,ni,nj)
2236 mat_loc2(ix,jx) = sij(1,ni,nj)
2237 ELSEIF (ki == 3)
THEN 2238 mat_loc1(ix,jx) = sij(5,ni,nj)
2239 mat_loc2(ix,jx) = sij(5,ni,nj)
2244 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), n_w2, jdxn(1:n_w2), &
2245 mat_loc1(1:3*n_ws1,1:n_w2), add_values, ierr)
2246 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), n_w2, jdxn(1:n_w2), &
2247 mat_loc2(1:3*n_ws1,1:n_w2), add_values, ierr)
2255 i = phi_mesh%jj(ni,m2)
2256 ib = la_phi%loc_to_glob(1,i)
2261 j = interface_h_phi%jjs1(nj,ms)
2262 jb = la_h%loc_to_glob(kj,j)
2263 jx = (kj-1)*n_ws1 + nj
2269 mat_loc1(ix,jx) = smuij(1,nj,ni)
2270 mat_loc2(ix,jx) = smuij(1,nj,ni)
2272 ELSEIF (kj == 3)
THEN 2276 mat_loc1(ix,jx) = smuij(5,nj,ni)
2277 mat_loc2(ix,jx) = smuij(5,nj,ni)
2284 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
2285 mat_loc1(1:n_w2,1:3*n_ws1), add_values, ierr)
2286 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
2287 mat_loc2(1:n_w2,1:3*n_ws1), add_values, ierr)
2302 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2303 x =
rjs(ls,ms2)*ray/sigma(m1)
2308 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2309 drmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(1,:,ls,ms1))
2310 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(2,:,ls,ms1))
2317 y = x *
wws(ni,ls)*w_cs(nj,ls)
2323 sij(1,ni,nj) = sij(1,ni,nj) + y*( mode/ray)**2*
rnorms(1,ls,ms2)/muhl
2324 sij(3,ni,nj) = sij(3,ni,nj) + y*( mode/ray**2)*
rnorms(1,ls,ms2)/muhl
2325 sij(4,ni,nj) = sij(4,ni,nj) + y*(-mode/ray**2)*
rnorms(1,ls,ms2)/muhl
2326 sij(5,ni,nj) = sij(5,ni,nj) + y*( mode/ray)**2*
rnorms(2,ls,ms2)/muhl
2328 sij(3,ni,nj) = sij(3,ni,nj) &
2329 + y*(mode/ray)*(drmuhl*
rnorms(1,ls,ms2)+dzmuhl*
rnorms(2,ls,ms2))*(-1/muhl**2)
2330 sij(4,ni,nj) = sij(4,ni,nj) &
2331 - y*(mode/ray)*(drmuhl*
rnorms(1,ls,ms2)+dzmuhl*
rnorms(2,ls,ms2))*(-1/muhl**2)
2345 i = interface_h_phi%jjs2(ni,ms)
2346 ib = la_phi%loc_to_glob(1,i)
2351 j = interface_h_phi%jjs1(nj,ms)
2352 jb = la_h%loc_to_glob(kj,j)
2353 jx = (kj-1)*n_ws1 + nj
2356 mat_loc1(ix,jx) = sij(1,ni,nj)
2357 mat_loc2(ix,jx) = sij(1,ni,nj)
2358 ELSEIF (kj == 2)
THEN 2359 mat_loc1(ix,jx) = sij(3,ni,nj)
2360 mat_loc2(ix,jx) = sij(4,ni,nj)
2361 ELSEIF (kj == 3)
THEN 2362 mat_loc1(ix,jx) = sij(5,ni,nj)
2363 mat_loc2(ix,jx) = sij(5,ni,nj)
2368 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), 3*n_ws1, jdxn(1:3*n_ws1), &
2369 mat_loc1(1:n_ws2,1:3*n_ws1), add_values, ierr)
2370 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), 3*n_ws1, jdxn(1:3*n_ws1), &
2371 mat_loc2(1:n_ws2,1:3*n_ws1), add_values, ierr)
2379 i = interface_h_phi%jjs1(ni,ms)
2380 ib = la_h%loc_to_glob(ki,i)
2381 ix = (ki-1)*n_ws1 + ni
2384 j = interface_h_phi%jjs2(nj,ms)
2385 jb = la_phi%loc_to_glob(1,j)
2389 mat_loc1(ix,jx) = sij(1,nj,ni)
2390 mat_loc2(ix,jx) = sij(1,nj,ni)
2391 ELSEIF (ki == 2)
THEN 2392 mat_loc1(ix,jx) = sij(3,nj,ni)
2393 mat_loc2(ix,jx) = sij(4,nj,ni)
2394 ELSEIF (ki == 3)
THEN 2395 mat_loc1(ix,jx) = sij(5,nj,ni)
2396 mat_loc2(ix,jx) = sij(5,nj,ni)
2401 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), n_ws2, jdxn(1:n_ws2), &
2402 mat_loc1(1:3*n_ws1,1:n_ws2), add_values, ierr)
2403 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), n_ws2, jdxn(1:n_ws2), &
2404 mat_loc2(1:3*n_ws1,1:n_ws2), add_values, ierr)
2411 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2412 x =
rjs(ls,ms2)*ray/sigma(m1)
2417 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2421 y = x*
wws(ni,ls)*mode/ray
2423 sij(3,ni,nj) = sij(3,ni,nj) + &
2426 y*(dw_cs(2,nj,ls,ms1)*
rnorms(2,ls,ms2)/muhl + dw_cs(1,nj,ls,ms1)*
rnorms(1,ls,ms2)/muhl)
2431 sij(4,:,:) = -sij(3,:,:)
2437 i = interface_h_phi%jjs2(ni,ms)
2438 ib = la_phi%loc_to_glob(1,i)
2441 j = h_mesh%jj(nj,m1)
2442 jb = la_h%loc_to_glob(kj,j)
2446 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_w1, jdxn(1:n_w1), &
2447 sij(3,1:n_ws2,1:n_w1), add_values, ierr)
2448 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_w1, jdxn(1:n_w1), &
2449 sij(4,1:n_ws2,1:n_w1), add_values, ierr)
2458 i = h_mesh%jj(ni,m1)
2459 ib = la_h%loc_to_glob(ki,i)
2462 j = interface_h_phi%jjs2(nj,ms)
2463 jb = la_phi%loc_to_glob(1,j)
2465 mat_loc1(ix,jx) = sij(3,nj,ni)
2466 mat_loc2(ix,jx) = sij(4,nj,ni)
2469 CALL matsetvalues(h_p_phi_mat1, n_w1, idxn(1:n_w1), n_ws2, jdxn(1:n_ws2), &
2470 mat_loc1(1:n_w1,1:n_ws2), add_values, ierr)
2471 CALL matsetvalues(h_p_phi_mat2, n_w1, idxn(1:n_w1), n_ws2, jdxn(1:n_ws2), &
2472 mat_loc2(1:n_w1,1:n_ws2), add_values, ierr)
2479 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2480 x =
rjs(ls,ms2)*ray/sigma(m1)
2485 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2489 y = x*(
dw_s(2,ni,ls,ms2)*
rnorms(1,ls,ms2) -
dw_s(1,ni,ls,ms2)*
rnorms(2,ls,ms2))
2494 sij(1,ni,nj) = sij(1,ni,nj) + y * (dw_cs(2,nj,ls,ms1)/muhl)
2495 sij(5,ni,nj) = sij(5,ni,nj) + (-y)*(dw_cs(1,nj,ls,ms1)/muhl)
2509 i = phi_mesh%jj(ni,m2)
2510 ib = la_phi%loc_to_glob(1,i)
2514 j = h_mesh%jj(nj,m1)
2516 jb = la_h%loc_to_glob(kj,j)
2517 jx = (kj-1)*n_w1 + nj
2520 mat_loc1(ix,jx) = sij(1,ni,nj)
2521 mat_loc2(ix,jx) = sij(1,ni,nj)
2522 ELSEIF (kj == 3)
THEN 2523 mat_loc1(ix,jx) = sij(5,ni,nj)
2524 mat_loc2(ix,jx) = sij(5,ni,nj)
2529 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 3*n_w1, jdxn(1:3*n_w1), &
2530 mat_loc1(1:n_w2,1:3*n_w1), add_values, ierr)
2531 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 3*n_w1, jdxn(1:3*n_w1), &
2532 mat_loc2(1:n_w2,1:3*n_w1), add_values, ierr)
2538 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2539 x =
rjs(ls,ms2)*ray/sigma(m1)
2543 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2544 drmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(1,:,ls,ms1))
2545 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(2,:,ls,ms1))
2548 y = x*(
dw_s(2,ni,ls,ms2)*
rnorms(1,ls,ms2) -
dw_s(1,ni,ls,ms2)*
rnorms(2,ls,ms2))
2550 sij(1,ni,nj) = sij(1,ni,nj) + y *( (-1/muhl**2)*dzmuhl)*w_cs(nj,ls)
2551 sij(5,ni,nj) = sij(5,ni,nj) + y *(-(-1/muhl**2)*drmuhl)*w_cs(nj,ls)
2563 i = phi_mesh%jj(ni,m2)
2564 ib = la_phi%loc_to_glob(1,i)
2569 j = interface_h_phi%jjs1(nj,ms)
2570 jb = la_h%loc_to_glob(kj,j)
2571 jx = (kj-1)*n_ws1 + nj
2574 mat_loc1(ix,jx) = sij(1,ni,nj)
2575 mat_loc2(ix,jx) = sij(1,ni,nj)
2576 ELSEIF (kj == 3)
THEN 2577 mat_loc1(ix,jx) = sij(5,ni,nj)
2578 mat_loc2(ix,jx) = sij(5,ni,nj)
2583 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
2584 mat_loc1(1:n_w2,1:3*n_ws1), add_values, ierr)
2585 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
2586 mat_loc2(1:n_w2,1:3*n_ws1), add_values, ierr)
2595 i = h_mesh%jj(ni,m1)
2596 ib = la_h%loc_to_glob(ki,i)
2597 ix = (ki-1)*n_w1 + ni
2600 j = phi_mesh%jj(nj,m2)
2601 jb = la_phi%loc_to_glob(1,j)
2605 mat_loc1(ix,jx) = sij(1,nj,ni)
2606 mat_loc2(ix,jx) = sij(1,nj,ni)
2607 ELSEIF (ki == 3)
THEN 2608 mat_loc1(ix,jx) = sij(5,nj,ni)
2609 mat_loc2(ix,jx) = sij(5,nj,ni)
2614 CALL matsetvalues(h_p_phi_mat1, 3*n_w1, idxn(1:3*n_w1), n_w2, jdxn(1:n_w2), &
2615 mat_loc1(1:3*n_w1,1:n_w2), add_values, ierr)
2616 CALL matsetvalues(h_p_phi_mat2, 3*n_w1, idxn(1:3*n_w1), n_w2, jdxn(1:n_w2), &
2617 mat_loc2(1:3*n_w1,1:n_w2), add_values, ierr)
2624 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2625 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2628 x = c_lap*
rjs(ls,ms2)*ray
2632 sij(1,ni,nj) = sij(1,ni,nj) - x*w_cs(nj,ls)*
wws(ni,ls)*
rnorms(1,ls,ms2)
2633 sij(5,ni,nj) = sij(5,ni,nj) - x*w_cs(nj,ls)*
wws(ni,ls)*
rnorms(2,ls,ms2)
2644 i = interface_h_phi%jjs2(ni,ms)
2645 ib = la_phi%loc_to_glob(1,i)
2649 j = interface_h_phi%jjs1(nj,ms)
2650 jb = la_h%loc_to_glob(1,j)
2653 mat_loc1(ix,jx) = sij(1,ni,nj)
2654 mat_loc2(ix,jx) = sij(1,ni,nj)
2656 jb = la_h%loc_to_glob(3,j)
2659 mat_loc1(ix,jx) = sij(5,ni,nj)
2660 mat_loc2(ix,jx) = sij(5,ni,nj)
2663 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), 2*n_ws1, jdxn(1:2*n_ws1), &
2664 mat_loc1(1:n_ws2,1:2*n_ws1), add_values, ierr)
2665 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), 2*n_ws1, jdxn(1:2*n_ws1), &
2666 mat_loc2(1:n_ws2,1:2*n_ws1), add_values, ierr)
2677 IF (stab(2) > 1.d-12)
THEN 2686 ms2 = interface_h_phi%mesh2(ms)
2689 hm1 =(sum(
rjs(:,ms2))/h_mesh%global_diameter)**(2*
alpha-1)/(
inputs%sigma_min*
inputs%mu_min**2*h_mesh%global_diameter)
2695 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2698 DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,ms2)
2699 ray = ray + phi_mesh%rr(1,i)* phi_mesh%gauss%wws(ni,ls)
2704 ray = stab_div*ray*hm1*
rjs(ls,ms2)
2712 x = muhl*w_cs(ni,ls)*w_cs(nj,ls)*ray
2714 hsij(1,ni,nj) = hsij(1,ni,nj) + x*
rnorms(1,ls,ms2)**2
2715 hsij(4,ni,nj) = hsij(4,ni,nj) + x*
rnorms(1,ls,ms2)*
rnorms(2,ls,ms2)
2716 hsij(6,ni,nj) = hsij(6,ni,nj) + x*
rnorms(2,ls,ms2)**2
2720 x = muhl*mu_phi*w_cs(ni,ls)*(
dw_s(1,nj,ls,ms2)*
rnorms(1,ls,ms2) +&
2722 sij(1,ni,nj) = sij(1,ni,nj) - x*
rnorms(1,ls,ms2)
2723 sij(5,ni,nj) = sij(5,ni,nj) - x*
rnorms(2,ls,ms2)
2730 x = mu_phi*w_cs(nj,ls)*(
dw_s(1,ni,ls,ms2)*
rnorms(1,ls,ms2) +&
2732 smuij(1,ni,nj) = smuij(1,ni,nj) - x*
rnorms(1,ls,ms2)
2733 smuij(5,ni,nj) = smuij(5,ni,nj) - x*
rnorms(2,ls,ms2)
2741 x = mu_phi**2*(
dw_s(1,ni,ls,ms2)*
rnorms(1,ls,ms2) +
dw_s(2,ni,ls,ms2)*
rnorms(2,ls,ms2))* &
2742 (
dw_s(1,nj,ls,ms2)*
rnorms(1,ls,ms2) +
dw_s(2,nj,ls,ms2)*
rnorms(2,ls,ms2))*ray
2743 phisij(ni,nj) = phisij(ni,nj) + x
2748 sij(2,:,:) = sij(1,:,:)
2749 sij(6,:,:) = sij(5,:,:)
2755 i = h_mesh%jjs(ni,ms1)
2757 ib = la_h%loc_to_glob(ki,i)
2758 ix = (ki/2)*n_ws1 + ni
2761 j = h_mesh%jjs(nj,ms1)
2763 jb = la_h%loc_to_glob(kj,j)
2764 jx = (kj/2)*n_ws1 + nj
2767 mat_loc1(ix,jx) = hsij(1,ni,nj)
2768 mat_loc2(ix,jx) = hsij(1,ni,nj)
2769 ELSE IF (ki*kj==9)
THEN 2770 mat_loc1(ix,jx) = hsij(6,ni,nj)
2771 mat_loc2(ix,jx) = hsij(6,ni,nj)
2772 ELSE IF (ki*kj==3)
THEN 2773 mat_loc1(ix,jx) = hsij(4,ni,nj)
2774 mat_loc2(ix,jx) = hsij(4,ni,nj)
2780 j = phi_mesh%jj(nj,m2)
2781 jb = la_phi%loc_to_glob(1,j)
2784 mat_loc1(ix,jx) = sij(2*ki-1,ni,nj)
2785 mat_loc2(ix,jx) = sij(2*ki-1,ni,nj)
2789 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), &
2790 mat_loc1(1:2*n_ws1,1:2*n_ws1+n_w2), add_values, ierr)
2791 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), &
2792 mat_loc2(1:2*n_ws1,1:2*n_ws1+n_w2), add_values, ierr)
2797 i = phi_mesh%jj(ni,m2)
2798 ib = la_phi%loc_to_glob(1,i)
2802 j = h_mesh%jjs(nj,ms1)
2804 jb = la_h%loc_to_glob(kj,j)
2805 jx = (kj/2)*n_ws1 + nj
2810 mat_loc1(ix,jx) = smuij(2*kj-1,ni,nj)
2811 mat_loc2(ix,jx) = smuij(2*kj-1,ni,nj)
2817 j = phi_mesh%jj(nj,m2)
2818 jb = la_phi%loc_to_glob(1,j)
2821 mat_loc1(ix,jx) = phisij(ni,nj)
2822 mat_loc2(ix,jx) = phisij(ni,nj)
2825 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2826 mat_loc1(1:n_w2,1:2*n_ws1+n_w2), add_values, ierr)
2827 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2828 mat_loc2(1:n_w2,1:2*n_ws1+n_w2), add_values, ierr)
2839 IF (.NOT.
PRESENT(index_fourier) .OR. .NOT.
PRESENT(r_fourier))
RETURN 2840 IF (r_fourier.GT.0.d0)
THEN 2842 DO ms = 1, phi_mesh%mes
2843 IF (phi_mesh%sides(ms) /= index_fourier) cycle
2847 DO ls = 1, phi_mesh%gauss%l_Gs
2850 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms))* phi_mesh%gauss%wws(:,ls))
2852 x = c_mu_phi*
rjs(ls,ms)*ray/r_fourier
2854 DO ni=1, phi_mesh%gauss%n_ws
2855 DO nj=1, phi_mesh%gauss%n_ws
2856 phisij(ni,nj) = phisij(ni,nj) + x*
wws(ni,ls)*
wws(nj,ls)
2863 DO ni = 1, phi_mesh%gauss%n_ws
2864 i = phi_mesh%jjs(ni,ms)
2865 ib = la_phi%loc_to_glob(1,i)
2867 DO nj = 1, phi_mesh%gauss%n_ws
2868 j = phi_mesh%jjs(nj,ms)
2869 jb = la_phi%loc_to_glob(1,j)
2873 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
2874 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
2875 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
2876 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
2880 CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
2881 CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
2882 CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
2883 CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
2890 mode, stab, mu_h_field, la_h, h_p_phi_mat1, h_p_phi_mat2, sigma_np, sigma)
2897 #include "petsc/finclude/petsc.h" 2901 INTEGER,
DIMENSION(:),
INTENT(IN) :: jj_v_to_H
2902 INTEGER,
DIMENSION(:),
INTENT(IN) :: Dirichlet_bdy_H_sides
2903 INTEGER,
INTENT(IN) :: mode
2904 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
2905 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field, sigma_np
2906 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN) :: sigma
2907 INTEGER :: ms, ls, ni, nj, i, j, &
2908 n_ws1, n_w1, m1, ki, kj, ib, jb
2909 REAL(KIND=8) :: x, y, hm1
2910 REAL(KIND=8) :: ray, error, stab_colle_h_mu
2911 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: Hsij
2931 REAL(KIND=8) :: muhl, dzmuhl, drmuhl, norm
2933 REAL(KIND=8) :: c_sym=.0d0
2938 REAL(KIND=8),
DIMENSION(3*H_mesh%gauss%n_w,3*H_mesh%gauss%n_w) :: mat_loc1, mat_loc2
2939 INTEGER ,
DIMENSION(3*H_mesh%gauss%n_w) :: idxn, jdxn
2945 REAL(KIND=8),
DIMENSION(2) :: Dmu_field
2946 REAL(KIND=8),
DIMENSION(2) :: gauss_pt
2947 INTEGER,
DIMENSION(1) :: gauss_pt_id
2948 REAL(KIND=8),
DIMENSION(1) :: dummy_mu_bar
2950 petscerrorcode :: ierr
2951 mat :: h_p_phi_mat1, h_p_phi_mat2
2954 stab_colle_h_mu = stab(3)
2961 n_ws1 = h_mesh%gauss%n_ws
2962 n_w1 = h_mesh%gauss%n_w
2970 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
2971 ms = dirichlet_bdy_h_sides(count)
2974 hm1 = stab_colle_h_mu/(sum(h_mesh%gauss%rjs(:,ms))*
inputs%sigma_min)
2976 m1 = h_mesh%neighs(ms)
2977 mesh_id = h_mesh%i_d(m1)
2986 DO ls = 1, h_mesh%gauss%l_Gs
2988 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2991 IF(
inputs%if_use_fem_integration_for_mu_bar)
THEN 2992 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2994 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2995 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2996 gauss_pt_id(1)=mesh_id
2997 dummy_mu_bar(:) = mu_bar_in_fourier_space(h_mesh,1,1,gauss_pt,gauss_pt_id)
2998 muhl=dummy_mu_bar(1)
3002 x = hm1*h_mesh%gauss%rjs(ls,ms)*ray /muhl
3004 DO ni = 1, h_mesh%gauss%n_ws
3005 DO nj = 1, h_mesh%gauss%n_ws
3006 y = x * h_mesh%gauss%wws(ni,ls)*h_mesh%gauss%wws(nj,ls)
3008 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(h_mesh%gauss%rnorms(2,ls,ms)**2)
3009 hsij(4,ni,nj) = hsij(4,ni,nj) - y*h_mesh%gauss%rnorms(1,ls,ms)*h_mesh%gauss%rnorms(2,ls,ms)
3010 hsij(5,ni,nj) = hsij(5,ni,nj) + y
3011 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(h_mesh%gauss%rnorms(1,ls,ms)**2)
3026 i = h_mesh%jjs(ni,ms)
3027 ib = la_h%loc_to_glob(ki,i)
3028 ix = ni + (ki-1)*n_ws1
3032 j = h_mesh%jjs(nj,ms)
3033 jb = la_h%loc_to_glob(kj,j)
3034 jx = nj + (kj-1)*n_ws1
3036 IF ((ki == 1) .AND. (kj == 1))
THEN 3037 mat_loc1(ix,jx) = hsij(1,ni,nj)
3038 mat_loc2(ix,jx) = hsij(1,ni,nj)
3039 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 3040 mat_loc1(ix,jx) = hsij(4,ni,nj)
3041 mat_loc2(ix,jx) = hsij(4,ni,nj)
3042 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 3043 mat_loc1(ix,jx) = hsij(4,nj,ni)
3044 mat_loc2(ix,jx) = hsij(4,nj,ni)
3045 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 3046 mat_loc1(ix,jx) = hsij(5,ni,nj)
3047 mat_loc2(ix,jx) = hsij(5,ni,nj)
3048 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 3049 mat_loc1(ix,jx) = hsij(6,ni,nj)
3050 mat_loc2(ix,jx) = hsij(6,ni,nj)
3057 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
3058 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
3059 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
3060 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
3071 DO ls = 1, h_mesh%gauss%l_Gs
3073 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3076 IF(
inputs%if_use_fem_integration_for_mu_bar)
THEN 3077 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3078 drmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*
dw_s(1,:,ls,ms))
3079 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*
dw_s(2,:,ls,ms))
3081 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3082 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3083 gauss_pt_id(1)=mesh_id
3084 dummy_mu_bar(:) = mu_bar_in_fourier_space(h_mesh,1,1,gauss_pt,gauss_pt_id)
3085 muhl=dummy_mu_bar(1)
3086 dmu_field = grad_mu_bar_in_fourier_space(gauss_pt,gauss_pt_id)
3087 drmuhl =dmu_field(1)
3088 dzmuhl =dmu_field(2)
3092 IF (jj_v_to_h(h_mesh%jj(1,m1)) == -1)
THEN 3093 x = h_mesh%gauss%rjs(ls,ms)*ray/sigma(m1)
3095 x = h_mesh%gauss%rjs(ls,ms)*ray/sum(sigma_np(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3103 y = x*h_mesh%gauss%wws(ni,ls)*h_mesh%gauss%wws(nj,ls)
3105 hsij(2,ni,nj) = hsij(2,ni,nj) + y * (-mode/ray)*(
rnorms(1,ls,ms))/muhl
3106 hsij(3,ni,nj) = hsij(3,ni,nj) + y * mode/ray *(
rnorms(1,ls,ms))/muhl
3107 hsij(5,ni,nj) = hsij(5,ni,nj) + y * (-1/ray) *(
rnorms(1,ls,ms))/muhl
3108 hsij(8,ni,nj) = hsij(8,ni,nj) + y * (-mode/ray)*(
rnorms(2,ls,ms))/muhl
3109 hsij(9,ni,nj) = hsij(9,ni,nj) + y * mode/ray *(
rnorms(2,ls,ms))/muhl
3111 hsij(1,ni,nj) = hsij(1,ni,nj) - y*(
rnorms(2,ls,ms))*(-(dzmuhl/muhl**2))
3112 hsij(4,ni,nj) = hsij(4,ni,nj) - y*(
rnorms(2,ls,ms))*( (drmuhl/muhl**2))
3113 hsij(5,ni,nj) = hsij(5,ni,nj) &
3114 + y*(
rnorms(1,ls,ms)*drmuhl+
rnorms(2,ls,ms)*dzmuhl)/muhl**2
3115 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(
rnorms(1,ls,ms))*( (drmuhl/muhl**2))
3116 hsij(7,ni,nj) = hsij(7,ni,nj) + y*(
rnorms(1,ls,ms))*(-(dzmuhl/muhl**2))
3133 i = h_mesh%jjs(ni,ms)
3134 ib = la_h%loc_to_glob(ki,i)
3135 ix = ni + (ki-1)*n_ws1
3139 j = h_mesh%jjs(nj,ms)
3140 jb = la_h%loc_to_glob(kj,j)
3141 jx = nj + (kj-1)*n_ws1
3143 IF ((ki == 1) .AND. (kj == 1))
THEN 3144 mat_loc1(ix,jx) = hsij(1,ni,nj)
3145 mat_loc2(ix,jx) = hsij(1,ni,nj)
3146 ELSE IF ( (ki == 1) .AND. (kj == 3))
THEN 3147 mat_loc1(ix,jx) = hsij(4,ni,nj)
3148 mat_loc2(ix,jx) = hsij(4,ni,nj)
3149 ELSE IF ( (ki == 2) .AND. (kj == 1))
THEN 3150 mat_loc1(ix,jx) = hsij(2,ni,nj)
3151 mat_loc2(ix,jx) = hsij(3,ni,nj)
3152 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 3153 mat_loc1(ix,jx) = hsij(5,ni,nj)
3154 mat_loc2(ix,jx) = hsij(5,ni,nj)
3155 ELSEIF ( (ki == 2) .AND. (kj == 3))
THEN 3156 mat_loc1(ix,jx) = hsij(8,ni,nj)
3157 mat_loc2(ix,jx) = hsij(9,ni,nj)
3158 ELSEIF ( (ki == 3) .AND. (kj == 1))
THEN 3159 mat_loc1(ix,jx) = hsij(7,ni,nj)
3160 mat_loc2(ix,jx) = hsij(7,ni,nj)
3161 ELSEIF ( (ki == 3) .AND. (kj == 3))
THEN 3162 mat_loc1(ix,jx) = hsij(6,ni,nj)
3163 mat_loc2(ix,jx) = hsij(6,ni,nj)
3170 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
3171 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
3172 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
3173 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
3181 i = h_mesh%jjs(ni,ms)
3182 ib = la_h%loc_to_glob(ki,i)
3183 ix = ni + (ki-1)*n_ws1
3187 j = h_mesh%jjs(nj,ms)
3188 jb = la_h%loc_to_glob(kj,j)
3189 jx = nj + (kj-1)*n_ws1
3191 IF ( (kj == 2) .AND. (ki == 1))
THEN 3192 mat_loc1(ix,jx) = hsij(2,nj,ni)
3193 mat_loc2(ix,jx) = hsij(3,nj,ni)
3194 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN 3195 mat_loc1(ix,jx) = hsij(5,nj,ni)
3196 mat_loc2(ix,jx) = hsij(5,nj,ni)
3197 ELSEIF ( (kj == 2) .AND. (ki == 3))
THEN 3198 mat_loc1(ix,jx) = hsij(8,nj,ni)
3199 mat_loc2(ix,jx) = hsij(9,nj,ni)
3205 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
3206 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
3207 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
3208 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
3214 DO ls = 1, h_mesh%gauss%l_Gs
3216 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3219 IF(
inputs%if_use_fem_integration_for_mu_bar)
THEN 3220 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3222 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3223 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3224 gauss_pt_id(1)=mesh_id
3225 dummy_mu_bar(:) = mu_bar_in_fourier_space(h_mesh,1,1,gauss_pt,gauss_pt_id)
3226 muhl=dummy_mu_bar(1)
3230 IF (jj_v_to_h(h_mesh%jj(1,m1)) == -1)
THEN 3231 x = h_mesh%gauss%rjs(ls,ms)*ray/(sigma(m1)*muhl)
3233 x = h_mesh%gauss%rjs(ls,ms)*ray/(sum(sigma_np(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))*muhl)
3240 y = x*h_mesh%gauss%wws(ni,ls)
3242 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(-h_mesh%gauss%dw_s(2,nj,ls,ms))*(
rnorms(2,ls,ms))
3243 hsij(4,ni,nj) = hsij(4,ni,nj) + y* h_mesh%gauss%dw_s(1,nj,ls,ms) *(
rnorms(2,ls,ms))
3244 hsij(5,ni,nj) = hsij(5,ni,nj) + &
3245 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)))
3246 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(-h_mesh%gauss%dw_s(1,nj,ls,ms))*(
rnorms(1,ls,ms))
3247 hsij(7,ni,nj) = hsij(7,ni,nj) + y* h_mesh%gauss%dw_s(2,nj,ls,ms) *(
rnorms(1,ls,ms))
3261 i = h_mesh%jjs(ni,ms)
3262 ib = la_h%loc_to_glob(ki,i)
3263 ix = ni + (ki-1)*n_ws1
3267 j = h_mesh%jj(nj,m1)
3268 jb = la_h%loc_to_glob(kj,j)
3269 jx = nj + (kj-1)*n_w1
3271 IF ((ki == 1) .AND. (kj == 1))
THEN 3272 mat_loc1(ix,jx) = hsij(1,ni,nj)
3273 mat_loc2(ix,jx) = hsij(1,ni,nj)
3274 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 3275 mat_loc1(ix,jx) = hsij(4,ni,nj)
3276 mat_loc2(ix,jx) = hsij(4,ni,nj)
3277 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 3278 mat_loc1(ix,jx) = hsij(5,ni,nj)
3279 mat_loc2(ix,jx) = hsij(5,ni,nj)
3280 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 3281 mat_loc1(ix,jx) = hsij(6,ni,nj)
3282 mat_loc2(ix,jx) = hsij(6,ni,nj)
3283 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 3284 mat_loc1(ix,jx) = hsij(7,ni,nj)
3285 mat_loc2(ix,jx) = hsij(7,ni,nj)
3292 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
3293 mat_loc1(1:3*n_ws1,1:3*n_w1), add_values, ierr)
3294 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
3295 mat_loc2(1:3*n_ws1,1:3*n_w1), add_values, ierr)
3303 i = h_mesh%jj(ni,m1)
3304 ib = la_h%loc_to_glob(ki,i)
3305 ix = ni + (ki-1)*n_w1
3309 j = h_mesh%jjs(nj,ms)
3310 jb = la_h%loc_to_glob(kj,j)
3311 jx = nj + (kj-1)*n_ws1
3313 IF ((kj == 1) .AND. (ki == 1))
THEN 3314 mat_loc1(ix,jx) = hsij(1,nj,ni)
3315 mat_loc2(ix,jx) = hsij(1,nj,ni)
3316 ELSEIF ((kj == 1) .AND. (ki == 3))
THEN 3317 mat_loc1(ix,jx) = hsij(4,nj,ni)
3318 mat_loc2(ix,jx) = hsij(4,nj,ni)
3319 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN 3320 mat_loc1(ix,jx) = hsij(5,nj,ni)
3321 mat_loc2(ix,jx) = hsij(5,nj,ni)
3322 ELSEIF ((kj == 3) .AND. (ki == 3))
THEN 3323 mat_loc1(ix,jx) = hsij(6,nj,ni)
3324 mat_loc2(ix,jx) = hsij(6,nj,ni)
3325 ELSEIF ((kj == 3) .AND. (ki == 1))
THEN 3326 mat_loc1(ix,jx) = hsij(7,nj,ni)
3327 mat_loc2(ix,jx) = hsij(7,nj,ni)
3334 CALL matsetvalues(h_p_phi_mat1, 3*n_w1 , idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
3335 mat_loc1(1:3*n_w1,1:3*n_ws1), add_values, ierr)
3336 CALL matsetvalues(h_p_phi_mat2, 3*n_w1 , idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
3337 mat_loc2(1:3*n_w1,1:3*n_ws1), add_values, ierr)
3345 /(h_mesh%global_diameter*
inputs%sigma_min*
inputs%mu_min**2)
3350 DO ls = 1, h_mesh%gauss%l_Gs
3352 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3355 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3356 x = norm*muhl*h_mesh%gauss%rjs(ls,ms)*ray
3360 DO ni = 1, h_mesh%gauss%n_ws
3361 DO nj = 1, h_mesh%gauss%n_ws
3362 y = x*h_mesh%gauss%wws(ni,ls)*h_mesh%gauss%wws(nj,ls)
3363 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(h_mesh%gauss%rnorms(1,ls,ms))**2
3364 hsij(4,ni,nj) = hsij(4,ni,nj) + y*h_mesh%gauss%rnorms(1,ls,ms)*h_mesh%gauss%rnorms(2,ls,ms)
3365 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(h_mesh%gauss%rnorms(2,ls,ms))**2
3366 hsij(7,ni,nj) = hsij(7,ni,nj) + y*h_mesh%gauss%rnorms(1,ls,ms)*h_mesh%gauss%rnorms(2,ls,ms)
3375 i = h_mesh%jjs(ni,ms)
3376 ib = la_h%loc_to_glob(ki,i)
3377 ix = ni + (ki-1)*n_ws1
3381 j = h_mesh%jjs(nj,ms)
3382 jb = la_h%loc_to_glob(kj,j)
3383 jx = nj + (kj-1)*n_ws1
3385 IF ((ki == 1) .AND. (kj == 1))
THEN 3386 mat_loc1(ix,jx) = hsij(1,ni,nj)
3387 mat_loc2(ix,jx) = hsij(1,ni,nj)
3388 ELSE IF ((ki == 1) .AND. (kj == 3))
THEN 3389 mat_loc1(ix,jx) = hsij(4,ni,nj)
3390 mat_loc2(ix,jx) = hsij(4,ni,nj)
3391 ELSE IF ((ki == 3) .AND. (kj == 1))
THEN 3392 mat_loc1(ix,jx) = hsij(7,ni,nj)
3393 mat_loc2(ix,jx) = hsij(7,ni,nj)
3394 ELSE IF ((ki == 3) .AND. (kj == 3))
THEN 3395 mat_loc1(ix,jx) = hsij(6,ni,nj)
3396 mat_loc2(ix,jx) = hsij(6,ni,nj)
3403 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
3404 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
3405 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
3406 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
3411 CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
3412 CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
3413 CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
3414 CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
3424 SUBROUTINE courant_int_by_parts(H_mesh,phi_mesh,interface_H_phi,sigma,mu_phi,mu_H_field,time,mode,&
3425 rhs_h,nl, la_h, la_phi, vb_1, vb_2, b_ext, h_pert, sigma_curl_gauss, j_over_sigma_gauss)
3435 #include "petsc/finclude/petsc.h" 3438 TYPE(
mesh_type),
INTENT(IN) :: H_mesh, phi_mesh
3440 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
3441 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma
3442 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
3443 INTEGER,
INTENT(IN) :: mode
3444 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: nl
3445 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: B_ext
3446 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rhs_H
3447 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: H_pert
3448 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: J_over_sigma_gauss
3449 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: sigma_curl_gauss
3451 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: src_H
3452 REAL(KIND=8),
DIMENSION(phi_mesh%np,2) :: src_phi
3453 REAL(KIND=8),
DIMENSION(6) :: c
3454 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_ws,phi_mesh%gauss%l_Gs) :: w_cs
3455 REAL(KIND=8),
DIMENSION(2) :: gaussp
3457 INTEGER :: m, l, i, ni, k, ms, ls, n_ws1, n_ws2, ms1, ms2, H_bloc_size, n_w2, m1
3459 REAL(KIND=8),
DIMENSION(6) :: JsolH_anal, rhs_Hl
3460 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%n_w) :: dwH
3463 REAL(KIND=8) :: ray_rjl, muhl, drmuhl, dzmuhl
3468 INTEGER,
DIMENSION(H_mesh%np) :: idxn_H
3469 INTEGER,
DIMENSION(phi_mesh%np) :: idxn_phi
3471 REAL(KIND=8),
DIMENSION(2) :: Dmu_field
3472 REAL(KIND=8),
DIMENSION(2) :: gauss_pt
3473 INTEGER,
DIMENSION(1) :: gauss_pt_id
3474 REAL(KIND=8),
DIMENSION(1) :: dummy_mu_bar
3477 REAL(KIND=8),
DIMENSION(6) :: B_ext_l
3478 petscerrorcode :: ierr
3483 CALL veczeroentries(vb_1, ierr)
3484 CALL veczeroentries(vb_2, ierr)
3496 mesh_id1 = h_mesh%i_d(m)
3497 DO l = 1, h_mesh%gauss%l_G
3502 muhl=sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
3504 dwh = h_mesh%gauss%dw(:,:,l,m)
3507 b_ext_l(k) = sum(b_ext(h_mesh%jj(:,m),k)*h_mesh%gauss%ww(:,l))
3513 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
3514 gaussp = gaussp + h_mesh%rr(:,i)*h_mesh%gauss%ww(ni,l)
3515 jsolh_anal(:) = jsolh_anal(:) + nl(i,:)*h_mesh%gauss%ww(ni,l)
3516 rhs_hl(:) = rhs_hl(:) + rhs_h(i,:)*h_mesh%gauss%ww(ni,l)
3519 ray_rjl = h_mesh%gauss%rj(l,m)*ray
3521 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 3523 jsolh_anal(k) = jsolh_anal(k) + j_over_sigma_gauss(index,k) + sigma_curl_gauss(index,k)
3528 jsolh_anal(k) = jsolh_anal(k) + &
3529 jexact_gauss(k, gaussp, mode, mu_phi, sigma(m), muhl, time, mesh_id1, b_ext_l)/sigma(m)
3534 IF (
inputs%if_permeability_variable_in_theta)
THEN 3536 IF(
inputs%if_use_fem_integration_for_mu_bar)
THEN 3537 muhl = sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
3538 drmuhl = sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%dw(1,:,l,m))
3539 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%dw(2,:,l,m))
3541 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
3542 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
3543 gauss_pt_id(1)=mesh_id1
3544 dummy_mu_bar(:) = mu_bar_in_fourier_space(h_mesh,1,1,gauss_pt,gauss_pt_id)
3545 muhl=dummy_mu_bar(1)
3546 dmu_field = grad_mu_bar_in_fourier_space(gauss_pt,gauss_pt_id)
3547 drmuhl =dmu_field(1)
3548 dzmuhl =dmu_field(2)
3552 DO ni = 1,h_mesh%gauss%n_w
3558 c(1) = c(1) + ( mode/ray*h_pert(i,6)*h_mesh%gauss%ww(ni,l) &
3559 - h_pert(i,3)*h_mesh%gauss%dw(2,ni,l,m)) &
3560 + (1/muhl)*( mode/ray*b_ext(i,6)*h_mesh%gauss%ww(ni,l) &
3561 - b_ext(i,3)*h_mesh%gauss%dw(2,ni,l,m)) &
3562 -(1/muhl**2)*(-dzmuhl*b_ext(i,3)*h_mesh%gauss%ww(ni,l))
3564 c(2) = c(2) + (-mode/ray*h_pert(i,5)*h_mesh%gauss%ww(ni,l) &
3565 - h_pert(i,4)*h_mesh%gauss%dw(2,ni,l,m)) &
3566 + (1/muhl)*(-mode/ray*b_ext(i,5)*h_mesh%gauss%ww(ni,l) &
3567 - b_ext(i,4)*h_mesh%gauss%dw(2,ni,l,m)) &
3568 -(1/muhl**2)*(-dzmuhl*b_ext(i,4)*h_mesh%gauss%ww(ni,l))
3571 c(3) = c(3) + (h_pert(i,1)*h_mesh%gauss%dw(2,ni,l,m) &
3572 - h_pert(i,5)*h_mesh%gauss%dw(1,ni,l,m)) &
3573 + (1/muhl)*(b_ext(i,1)*h_mesh%gauss%dw(2,ni,l,m) &
3574 - b_ext(i,5)*h_mesh%gauss%dw(1,ni,l,m)) &
3575 -(1/muhl**2)*(dzmuhl*b_ext(i,1)-drmuhl*b_ext(i,5))*h_mesh%gauss%ww(ni,l)
3577 c(4) = c(4) + (h_pert(i,2)*h_mesh%gauss%dw(2,ni,l,m) &
3578 - h_pert(i,6)*h_mesh%gauss%dw(1,ni,l,m)) &
3579 + (1/muhl)*(b_ext(i,2)*h_mesh%gauss%dw(2,ni,l,m) &
3580 - b_ext(i,6)*h_mesh%gauss%dw(1,ni,l,m)) &
3581 -(1/muhl**2)*(dzmuhl*b_ext(i,2)-drmuhl*b_ext(i,6))*h_mesh%gauss%ww(ni,l)
3584 c(5) = c(5) + (h_pert(i,3)*h_mesh%gauss%dw(1,ni,l,m) &
3585 + h_pert(i,3)*h_mesh%gauss%ww(ni,l)/ray &
3586 - mode/ray*h_pert(i,2)*h_mesh%gauss%ww(ni,l)) &
3587 + (1/muhl)*(b_ext(i,3)*h_mesh%gauss%dw(1,ni,l,m) &
3588 + b_ext(i,3)*h_mesh%gauss%ww(ni,l)/ray &
3589 - mode/ray*b_ext(i,2)*h_mesh%gauss%ww(ni,l)) &
3590 -(1/muhl**2)*(drmuhl*b_ext(i,3)*h_mesh%gauss%ww(ni,l))
3592 c(6) = c(6) + (h_pert(i,4)*h_mesh%gauss%dw(1,ni,l,m) &
3593 + h_pert(i,4)*h_mesh%gauss%ww(ni,l)/ray &
3594 + mode/ray*h_pert(i,1)*h_mesh%gauss%ww(ni,l)) &
3595 + (1/muhl)*(b_ext(i,4)*h_mesh%gauss%dw(1,ni,l,m) &
3596 + b_ext(i,4)*h_mesh%gauss%ww(ni,l)/ray &
3597 + mode/ray*b_ext(i,1)*h_mesh%gauss%ww(ni,l)) &
3598 -(1/muhl**2)*(drmuhl*b_ext(i,4)*h_mesh%gauss%ww(ni,l))
3602 jsolh_anal= jsolh_anal + c
3606 DO ni = 1,h_mesh%gauss%n_w
3611 src_h(i,1) = src_h(i,1)+ ray_rjl &
3612 *(jsolh_anal(3)*dwh(2,ni) &
3613 + mode/ray*jsolh_anal(6)*h_mesh%gauss%ww(ni,l) &
3614 + rhs_hl(1)*h_mesh%gauss%ww(ni,l))
3616 src_h(i,2) = src_h(i,2)+ ray_rjl &
3617 *(jsolh_anal(4)*dwh(2,ni) &
3618 - mode/ray*jsolh_anal(5)*h_mesh%gauss%ww(ni,l) &
3619 + rhs_hl(2)*h_mesh%gauss%ww(ni,l))
3622 src_h(i,3) = src_h(i,3)+ ray_rjl &
3623 * (-jsolh_anal(1)*dwh(2,ni) &
3624 + 1/ray*jsolh_anal(5)*(ray*dwh(1,ni) + h_mesh%gauss%ww(ni,l)) &
3625 + rhs_hl(3)*h_mesh%gauss%ww(ni,l))
3627 src_h(i,4) = src_h(i,4)+ ray_rjl &
3628 * (-jsolh_anal(2)*dwh(2,ni) &
3629 + 1/ray*jsolh_anal(6)*(ray*dwh(1,ni) + h_mesh%gauss%ww(ni,l)) &
3630 + rhs_hl(4)*h_mesh%gauss%ww(ni,l))
3633 src_h(i,5) = src_h(i,5)+ ray_rjl* &
3634 (-mode/ray*jsolh_anal(2)*h_mesh%gauss%ww(ni,l) &
3635 - jsolh_anal(3)*dwh(1,ni) &
3636 + rhs_hl(5)*h_mesh%gauss%ww(ni,l))
3638 src_h(i,6) = src_h(i,6)+ ray_rjl* &
3639 (mode/ray*jsolh_anal(1)*h_mesh%gauss%ww(ni,l) &
3640 - jsolh_anal(4)*dwh(1,ni) &
3641 + rhs_hl(6)*h_mesh%gauss%ww(ni,l))
3698 CALL gauss(phi_mesh)
3700 n_ws1 = h_mesh%gauss%n_ws
3701 n_ws2 = phi_mesh%gauss%n_ws
3702 n_w2 = phi_mesh%gauss%n_w
3704 h_bloc_size = h_mesh%np
3706 IF (interface_h_phi%mes /=0)
THEN 3707 IF (h_mesh%gauss%n_ws ==
n_ws)
THEN 3711 w_cs(1,ls)=
wws(1,ls)+0.5*
wws(3,ls)
3712 w_cs(2,ls)=
wws(2,ls)+0.5*
wws(3,ls)
3720 DO ms = 1, interface_h_phi%mes
3722 ms2 = interface_h_phi%mesh2(ms)
3723 ms1 = interface_h_phi%mesh1(ms)
3724 m = phi_mesh%neighs(ms2)
3725 m1 = h_mesh%neighs(ms1)
3726 mesh_id1 = h_mesh%i_d(m1)
3730 muhl=sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
3733 b_ext_l(k) = sum(b_ext(interface_h_phi%jjs1(1:n_ws1,ms),k)*w_cs(1:n_ws1,ls))
3738 DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,interface_h_phi%mesh2(ms))
3739 ray = ray + phi_mesh%rr(1,i)*
wws(ni,ls)
3744 i=phi_mesh%jjs(ni,ms2)
3745 gaussp = gaussp + phi_mesh%rr(:,i)*phi_mesh%gauss%wws(ni,ls)
3749 jsolh_anal(k) = jexact_gauss(k, gaussp, mode, mu_phi ,sigma(m1), &
3750 muhl, time, mesh_id1, b_ext_l)/sigma(m1) &
3751 + sum(nl(h_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
3771 i = interface_h_phi%jjs1(ni,ms)
3772 src_h(i,1) = src_h(i,1)+
rjs(ls,ms2)*ray*( &
3773 -jsolh_anal(3)*w_cs(ni,ls)*(-
rnorms(2,ls,ms2)))
3775 src_h(i,2) = src_h(i,2)+
rjs(ls,ms2)*ray*( &
3776 -jsolh_anal(4)*w_cs(ni,ls)*(-
rnorms(2,ls,ms2)))
3778 src_h(i,3) = src_h(i,3)+
rjs(ls,ms2)*ray*( &
3779 jsolh_anal(1)*w_cs(ni,ls)*(-
rnorms(2,ls,ms2)) &
3780 -jsolh_anal(5)*w_cs(ni,ls)*(-
rnorms(1,ls,ms2)))
3782 src_h(i,4) = src_h(i,4)+
rjs(ls,ms2)*ray*( &
3783 jsolh_anal(2)*w_cs(ni,ls)*(-
rnorms(2,ls,ms2)) &
3784 -jsolh_anal(6)*w_cs(ni,ls)*(-
rnorms(1,ls,ms2)))
3786 src_h(i,5) = src_h(i,5)+
rjs(ls,ms2)*ray*( &
3787 jsolh_anal(3)*w_cs(ni,ls)*(-
rnorms(1,ls,ms2)))
3789 src_h(i,6) = src_h(i,6)+
rjs(ls,ms2)*ray*( &
3790 jsolh_anal(4)*w_cs(ni,ls)*(-
rnorms(1,ls,ms2)))
3796 i = interface_h_phi%jjs2(ni,ms)
3799 src_phi(i,1) = src_phi(i,1)+
rjs(ls,ms2)*( &
3800 - mode*jsolh_anal(2)*
wws(ni,ls) *
rnorms(2,ls,ms2) &
3801 + mode*jsolh_anal(6)*
wws(ni,ls) *
rnorms(1,ls,ms2))
3803 src_phi(i,2) = src_phi(i,2)+
rjs(ls,ms2)*( &
3804 + mode*jsolh_anal(1)*
wws(ni,ls) *
rnorms(2,ls,ms2) &
3805 - mode*jsolh_anal(5)*
wws(ni,ls) *
rnorms(1,ls,ms2))
3811 i = phi_mesh%jj(ni,m)
3812 src_phi(i,1) = src_phi(i,1)+
rjs(ls,ms2)*ray*( &
3813 + jsolh_anal(3) *(
dw_s(2,ni,ls,ms2) *
rnorms(1,ls,ms2)&
3816 src_phi(i,2) = src_phi(i,2)+
rjs(ls,ms2)*ray*( &
3817 + jsolh_anal(4)*(
dw_s(2,ni,ls,ms2) *
rnorms(1,ls,ms2)&
3825 i = interface_h_phi%jjs1(ni,ms)
3826 rhs_hl(:) = rhs_hl(:) + rhs_h(i,:)*w_cs(ni,ls)
3830 i = interface_h_phi%jjs2(ni,ms)
3831 src_phi(i,1) = src_phi(i,1)+
rjs(ls,ms2)*ray*
wws(ni,ls)*( &
3832 rhs_hl(1)*
rnorms(1,ls,ms2) + rhs_hl(5)*
rnorms(2,ls,ms2))
3833 src_phi(i,2) = src_phi(i,2)+
rjs(ls,ms2)*ray*
wws(ni,ls)*( &
3834 rhs_hl(2)*
rnorms(1,ls,ms2) + rhs_hl(6)*
rnorms(2,ls,ms2))
3842 IF (h_mesh%np /= 0)
THEN 3844 idxn_h = la_h%loc_to_glob(1,:)-1
3845 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,1), add_values, ierr)
3846 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,2), add_values, ierr)
3847 idxn_h = la_h%loc_to_glob(2,:)-1
3848 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,4), add_values, ierr)
3849 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,3), add_values, ierr)
3850 idxn_h = la_h%loc_to_glob(3,:)-1
3851 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,5), add_values, ierr)
3852 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,6), add_values, ierr)
3855 IF (phi_mesh%np /=0)
THEN 3857 idxn_phi = la_phi%loc_to_glob(1,:)-1
3858 CALL vecsetvalues(vb_1, phi_mesh%np, idxn_phi, src_phi(:,1), add_values, ierr)
3859 CALL vecsetvalues(vb_2, phi_mesh%np, idxn_phi, src_phi(:,2), add_values, ierr)
3863 CALL vecassemblybegin(vb_1,ierr)
3864 CALL vecassemblyend(vb_1,ierr)
3865 CALL vecassemblybegin(vb_2,ierr)
3866 CALL vecassemblyend(vb_2,ierr)
3877 SUBROUTINE surf_int(H_mesh, phi_mesh, pmag_mesh, interface_H_phi, interface_H_mu, &
3878 list_dirichlet_sides_h, sigma,mu_phi, mu_h_field, time, mode, &
3879 la_h, la_phi, la_pmag, vb_1, vb_2, sigma_tot_gauss, r_fourier, index_fourier)
3885 #include "petsc/finclude/petsc.h" 3888 TYPE(
mesh_type),
INTENT(IN) :: H_mesh, phi_mesh, pmag_mesh
3889 TYPE(
interface_type),
INTENT(IN) :: interface_H_phi, interface_H_mu
3890 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_dirichlet_sides_H
3891 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
3892 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN):: sigma
3893 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
3894 INTEGER,
INTENT(IN) :: mode
3895 REAL(KIND=8),
DIMENSION(:,:) :: sigma_tot_gauss
3896 REAL(KIND=8),
OPTIONAL :: R_fourier
3897 INTEGER,
OPTIONAL :: index_fourier
3898 REAL(KIND=8),
DIMENSION(H_mesh%np, 6) :: src_H
3899 REAL(KIND=8),
DIMENSION(phi_mesh%np, 2) :: src_phi
3900 REAL(KIND=8),
DIMENSION(pmag_mesh%np, 2) :: src_pmag
3901 REAL(KIND=8),
DIMENSION(4,H_mesh%gauss%l_Gs):: Banal
3902 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%l_Gs):: rloc
3903 REAL(KIND=8),
DIMENSION(H_mesh%gauss%l_Gs) :: B_dot_n_cos, B_dot_n_sin, muloc
3904 REAL(KIND=8),
DIMENSION(4,pmag_mesh%gauss%l_Gs):: Banal_pmesh
3905 REAL(KIND=8),
DIMENSION(2,pmag_mesh%gauss%l_Gs):: rloc_pmesh
3906 REAL(KIND=8),
DIMENSION(pmag_mesh%gauss%l_Gs) :: B_dot_n_cos_pmesh, B_dot_n_sin_pmesh, muloc_pmesh
3907 REAL(KIND=8) :: ray, muhl, norm, y
3908 INTEGER :: ms, ls, ns, i, k, m, n, ni, count, index
3909 REAL(KIND=8),
DIMENSION(2) :: gaussp
3910 REAL(KIND=8),
DIMENSION(6) :: EsolH_anal, Esolphi_anal
3911 INTEGER,
DIMENSION(H_mesh%np) :: idxn_H
3912 INTEGER,
DIMENSION(phi_mesh%np) :: idxn_phi
3913 INTEGER,
DIMENSION(pmag_mesh%np) :: idxn_pmag
3915 petscerrorcode :: ierr
3938 m = h_mesh%neighs(ms)
3940 DO ls = 1, h_mesh%gauss%l_Gs
3941 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3943 rloc(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3944 rloc(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3946 banal(1,:) = muhl*hexact(h_mesh, 1, rloc, mode, muloc, time)
3947 banal(2,:) = muhl*hexact(h_mesh, 2, rloc, mode, muloc, time)
3948 banal(3,:) = muhl*hexact(h_mesh, 5, rloc, mode, muloc, time)
3949 banal(4,:) = muhl*hexact(h_mesh, 6, rloc, mode, muloc, time)
3950 b_dot_n_cos = banal(1,:)*h_mesh%gauss%rnorms(1,:,ms) + banal(3,:)*h_mesh%gauss%rnorms(2,:,ms)
3951 b_dot_n_sin = banal(2,:)*h_mesh%gauss%rnorms(1,:,ms) + banal(4,:)*h_mesh%gauss%rnorms(2,:,ms)
3954 DO ls = 1, h_mesh%gauss%l_Gs
3957 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3962 DO ni = 1, h_mesh%gauss%n_ws; i = h_mesh%jjs(ni,ms)
3963 ray = ray + h_mesh%rr(1,i)* h_mesh%gauss%wws(ni,ls)
3966 IF (ray.LT.1.d-12*h_mesh%global_diameter) cycle
3969 DO ns=1, h_mesh%gauss%n_ws
3971 gaussp = gaussp + h_mesh%rr(:,i)*h_mesh%gauss%wws(ns,ls)
3979 esolh_anal(k) = eexact_gauss(k,gaussp,mode,mu_phi,sigma_tot_gauss(index,mod(k+1,2)+1),muhl, time)
3997 DO ns=1, h_mesh%gauss%n_ws
3998 i = h_mesh%jjs(ns,ms)
3999 src_h(i,1) = src_h(i,1)-h_mesh%gauss%rjs(ls,ms)*ray*( &
4000 -esolh_anal(3)*h_mesh%gauss%wws(ns,ls)* &
4001 (h_mesh%gauss%rnorms(2,ls,ms)))
4003 src_h(i,2) = src_h(i,2)-h_mesh%gauss%rjs(ls,ms)*ray*( &
4004 -esolh_anal(4)*h_mesh%gauss%wws(ns,ls)* &
4005 (h_mesh%gauss%rnorms(2,ls,ms)))
4007 src_h(i,3) = src_h(i,3)-h_mesh%gauss%rjs(ls,ms)*ray*( &
4008 esolh_anal(1)*h_mesh%gauss%wws(ns,ls)* &
4009 (h_mesh%gauss%rnorms(2,ls,ms)) - &
4010 esolh_anal(5)*h_mesh%gauss%wws(ns,ls) * &
4011 (h_mesh%gauss%rnorms(1,ls,ms)))
4013 src_h(i,4) = src_h(i,4)-h_mesh%gauss%rjs(ls,ms)*ray*( &
4014 esolh_anal(2)*h_mesh%gauss%wws(ns,ls) * &
4015 (h_mesh%gauss%rnorms(2,ls,ms)) - &
4016 esolh_anal(6)*h_mesh%gauss%wws(ns,ls) * &
4017 (h_mesh%gauss%rnorms(1,ls,ms)))
4019 src_h(i,5) = src_h(i,5)-h_mesh%gauss%rjs(ls,ms)*ray*( &
4020 esolh_anal(3)*h_mesh%gauss%wws(ns,ls)* &
4021 (h_mesh%gauss%rnorms(1,ls,ms)))
4023 src_h(i,6) = src_h(i,6)-h_mesh%gauss%rjs(ls,ms)*ray*( &
4024 esolh_anal(4)*h_mesh%gauss%wws(ns,ls) * &
4025 (h_mesh%gauss%rnorms(1,ls,ms)))
4030 norm = (
inputs%stab(1)/
inputs%Rem)*(sum(h_mesh%gauss%rjs(:,ms))/h_mesh%global_diameter)**(2*
alpha-1)&
4031 /(h_mesh%global_diameter*
inputs%sigma_min*
inputs%mu_min**2)
4035 DO ns = 1, h_mesh%gauss%n_ws
4036 i = h_mesh%jjs(ns,ms)
4037 y = norm*muhl*h_mesh%gauss%rjs(ls,ms)*ray
4038 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)
4039 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)
4040 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)
4041 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)
4051 m = pmag_mesh%neighs(ms)
4053 DO ls = 1, pmag_mesh%gauss%l_Gs
4054 muhl = sum(mu_h_field(pmag_mesh%jjs(:,ms))*pmag_mesh%gauss%wws(:,ls))
4055 muloc_pmesh(ls) = muhl
4056 rloc_pmesh(1,ls) = sum(pmag_mesh%rr(1,pmag_mesh%jjs(:,ms))*pmag_mesh%gauss%wws(:,ls))
4057 rloc_pmesh(2,ls) = sum(pmag_mesh%rr(2,pmag_mesh%jjs(:,ms))*pmag_mesh%gauss%wws(:,ls))
4059 banal_pmesh(1,:) = muhl*hexact(pmag_mesh, 1, rloc_pmesh, mode, muloc_pmesh, time)
4060 banal_pmesh(2,:) = muhl*hexact(pmag_mesh, 2, rloc_pmesh, mode, muloc_pmesh, time)
4061 banal_pmesh(3,:) = muhl*hexact(pmag_mesh, 5, rloc_pmesh, mode, muloc_pmesh, time)
4062 banal_pmesh(4,:) = muhl*hexact(pmag_mesh, 6, rloc_pmesh, mode, muloc_pmesh, time)
4063 b_dot_n_cos_pmesh = banal_pmesh(1,:)*pmag_mesh%gauss%rnorms(1,:,ms) &
4064 + banal_pmesh(3,:)*pmag_mesh%gauss%rnorms(2,:,ms)
4065 b_dot_n_sin_pmesh = banal_pmesh(2,:)*pmag_mesh%gauss%rnorms(1,:,ms) &
4066 + banal_pmesh(4,:)*pmag_mesh%gauss%rnorms(2,:,ms)
4069 DO ls = 1, pmag_mesh%gauss%l_Gs
4071 muhl = sum(mu_h_field(pmag_mesh%jjs(:,ms))*pmag_mesh%gauss%wws(:,ls))
4076 DO ni = 1, pmag_mesh%gauss%n_ws; i = pmag_mesh%jjs(ni,ms)
4077 ray = ray + pmag_mesh%rr(1,i)* pmag_mesh%gauss%wws(ni,ls)
4080 IF (ray.LT.1.d-12*h_mesh%global_diameter) cycle
4084 DO ns=1, pmag_mesh%gauss%n_ws
4085 i = pmag_mesh%jjs(ns,ms)
4086 src_pmag(i,1) = src_pmag(i,1) -
inputs%stab(1)*b_dot_n_cos_pmesh(ls)* &
4087 pmag_mesh%gauss%wws(ns,ls)*pmag_mesh%gauss%rjs(ls,ms)*ray
4088 src_pmag(i,2) = src_pmag(i,2) -
inputs%stab(1)*b_dot_n_sin_pmesh(ls)* &
4089 pmag_mesh%gauss%wws(ns,ls)*pmag_mesh%gauss%rjs(ls,ms)*ray
4099 m = phi_mesh%neighs(ms)
4100 DO ls = 1, phi_mesh%gauss%l_Gs
4103 DO ni = 1, phi_mesh%gauss%n_ws; i = phi_mesh%jjs(ni,ms)
4104 ray = ray + phi_mesh%rr(1,i)* phi_mesh%gauss%wws(ni,ls)
4108 DO ns=1, phi_mesh%gauss%n_ws
4109 i=phi_mesh%jjs(ns,ms)
4110 gaussp = gaussp + phi_mesh%rr(:,i)*phi_mesh%gauss%wws(ns,ls)
4116 esolphi_anal(k) = eexact_gauss(k,gaussp,mode,mu_phi,sigma(1),mu_h_field(1),time)
4123 DO ns=1, phi_mesh%gauss%n_ws
4124 i = phi_mesh%jjs(ns,ms)
4125 DO n = 1, phi_mesh%gauss%n_w
4126 IF (phi_mesh%jj(n,m) == i)
EXIT 4129 src_phi(i,1) = src_phi(i,1)-phi_mesh%gauss%rjs(ls,ms)*ray*( &
4130 +esolphi_anal(3)*(phi_mesh%gauss%dw_s(2,n,ls,ms)*phi_mesh%gauss%rnorms(1,ls,ms) &
4131 -phi_mesh%gauss%dw_s(1,n,ls,ms)*phi_mesh%gauss%rnorms(2,ls,ms))) &
4132 -phi_mesh%gauss%rjs(ls,ms)*(&
4133 -mode*esolphi_anal(2)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(2,ls,ms) &
4134 +mode*esolphi_anal(6)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(1,ls,ms))
4136 src_phi(i,2) = src_phi(i,2)-phi_mesh%gauss%rjs(ls,ms)*ray*( &
4137 +esolphi_anal(4)*(phi_mesh%gauss%dw_s(2,n,ls,ms)*phi_mesh%gauss%rnorms(1,ls,ms) &
4138 -phi_mesh%gauss%dw_s(1,n,ls,ms)*phi_mesh%gauss%rnorms(2,ls,ms))) &
4139 -phi_mesh%gauss%rjs(ls,ms)*(&
4140 mode*esolphi_anal(1)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(2,ls,ms) &
4141 -mode*esolphi_anal(5)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(1,ls,ms))
4158 IF (
PRESENT(r_fourier))
THEN 4159 IF (r_fourier.GE.0.d0)
CALL error_petsc(
'maxwell_update_time_with_B: R_fourier should be -1')
4182 IF (h_mesh%mes/=0)
THEN 4184 idxn_h = la_h%loc_to_glob(1,:)-1
4185 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,1), add_values, ierr)
4186 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,2), add_values, ierr)
4187 idxn_h = la_h%loc_to_glob(2,:)-1
4188 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,4), add_values, ierr)
4189 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,3), add_values, ierr)
4190 idxn_h = la_h%loc_to_glob(3,:)-1
4191 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,5), add_values, ierr)
4192 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,6), add_values, ierr)
4195 idxn_pmag = la_pmag%loc_to_glob(1,:)-1
4196 CALL vecsetvalues(vb_1, pmag_mesh%np, idxn_pmag, src_pmag(:,1), add_values, ierr)
4197 CALL vecsetvalues(vb_2, pmag_mesh%np, idxn_pmag, src_pmag(:,2), add_values, ierr)
4202 IF (phi_mesh%mes/=0)
THEN 4204 idxn_phi = la_phi%loc_to_glob(1,:)-1
4205 CALL vecsetvalues(vb_1, phi_mesh%np, idxn_phi, src_phi(:,1), add_values, ierr)
4206 CALL vecsetvalues(vb_2, phi_mesh%np, idxn_phi, src_phi(:,2), add_values, ierr)
4210 CALL vecassemblybegin(vb_1,ierr)
4211 CALL vecassemblyend(vb_1,ierr)
4212 CALL vecassemblybegin(vb_2,ierr)
4213 CALL vecassemblyend(vb_2,ierr)
4218 count=index_fourier; count=interface_h_mu%mes; count=interface_h_phi%mes
4219 count=
SIZE(list_dirichlet_sides_h)
4223 SUBROUTINE mat_maxwell_mu(H_mesh, jj_v_to_H, interface_H_mu, mode, stab, &
4224 mu_h_field, sigma, la_h, h_p_phi_mat1, h_p_phi_mat2, sigma_np)
4229 #include "petsc/finclude/petsc.h" 4233 INTEGER,
DIMENSION(:),
INTENT(IN) :: jj_v_to_H
4235 INTEGER,
INTENT(IN) :: mode
4236 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
4237 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma, mu_H_field
4238 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_np
4239 INTEGER :: ms, ls, ni, nj, i, j, &
4240 n_ws1, n_ws2, n_w2, n_w1, m1, m2, ki, kj,ib,jb, ms1, ms2
4241 REAL(KIND=8) :: x, y, z, norm, hm1
4242 REAL(KIND=8) :: ray, stab_colle_H_mu
4243 LOGICAL :: mark=.false.
4244 REAL(KIND=8),
DIMENSION(9,SIZE(H_mesh%jj,1),SIZE(H_mesh%jj,1),2,2) :: Hsij, Gsij
4259 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_ws,H_mesh%gauss%l_Gs) :: w_cs
4260 REAL(KIND=8),
DIMENSION(2, H_mesh%gauss%n_w, H_mesh%gauss%l_Gs, H_mesh%mes) :: dw_cs
4261 REAL(KIND=8),
DIMENSION(2, H_mesh%gauss%n_w) :: dwsi,dwsj
4262 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%l_Gs) :: gauss1, gauss2
4268 REAL(KIND=8),
DIMENSION(2) :: normi, normj
4269 REAL(KIND=8),
DIMENSION(SIZE(H_mesh%jjs,1)) :: wwsi, wwsj
4270 INTEGER :: n_wsi, n_wsj, ci, cj, n_wi, n_wj
4273 REAL(KIND=8) :: ref, diff, mu_H, muhl1, muhl2, muhi, muhj, sigmai, sigmaj
4275 REAL(KIND=8) :: c_sym =.0d0
4276 REAL(KIND=8) :: wwiwwj, normt, stab_div
4278 REAL(KIND=8) :: drmuhl1, drmuhl2, drmuhj, dzmuhl1, dzmuhl2, dzmuhj
4279 REAL(KIND=8) :: sigmal1, sigmal2
4285 REAL(KIND=8),
DIMENSION(6*H_mesh%gauss%n_w,6*H_mesh%gauss%n_w) :: mat_loc1, mat_loc2
4286 INTEGER ,
DIMENSION(6*H_mesh%gauss%n_w) :: idxn, jdxn
4290 mat :: h_p_phi_mat1, h_p_phi_mat2
4291 petscerrorcode :: ierr
4294 stab_colle_h_mu = stab(3)
4304 n_ws1 = h_mesh%gauss%n_ws
4305 n_ws2 = h_mesh%gauss%n_ws
4306 n_w1 = h_mesh%gauss%n_w
4307 n_w2 = h_mesh%gauss%n_w
4319 DO ms = 1, interface_h_mu%mes
4321 ms2 = interface_h_mu%mesh2(ms)
4322 m2 = h_mesh%neighs(ms2)
4323 ms1 = interface_h_mu%mesh1(ms)
4324 m1 = h_mesh%neighs(ms1)
4326 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
4327 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(1,ms2)))**2)
4328 IF (diff/ref .LT. 1.d-10)
THEN 4332 w_cs(1,ls)=
wws(2,ls)
4333 w_cs(2,ls)=
wws(1,ls)
4334 IF (n_ws1==3) w_cs(n_ws1,ls) =
wws(n_ws1,ls)
4335 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 4340 gauss2(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))*h_mesh%gauss%wws(:,ls))
4341 gauss2(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms2))*h_mesh%gauss%wws(:,ls))
4342 gauss1(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms1))*h_mesh%gauss%wws(:,ls))
4343 gauss1(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms1))*h_mesh%gauss%wws(:,ls))
4347 ref = sqrt(1.d-8+sum(gauss2(:,ls2)**2))
4350 diff = sqrt(sum((gauss2(:,ls2)-gauss1(:,ls1))**2))
4351 IF (diff .LT. 1.d-10)
THEN 4352 dw_cs(:,:,ls2,ms1) = h_mesh%gauss%dw_s(:,:,ls1,ms1)
4357 IF (.NOT.mark)
WRITE(*,*)
' BUG ' 4362 DO ms = 1, interface_h_mu%mes
4364 ms2 = interface_h_mu%mesh2(ms)
4365 ms1 = interface_h_mu%mesh1(ms)
4366 m2 = h_mesh%neighs(ms2)
4367 m1 = h_mesh%neighs(ms1)
4368 mu_h = sum(mu_h_field(h_mesh%jj(:,m1)))/h_mesh%gauss%n_w
4371 hm1 = 1/sum(
rjs(:,ms2))
4383 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
4384 x = hm1*
rjs(ls,ms2)*ray
4387 muhl1 = sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4388 muhl2 = sum(mu_h_field(h_mesh%jjs(:,ms2))*
wws(:,ls))
4392 normt =stab_colle_h_mu/
inputs%sigma_min
4399 norm = stab_div*(sum(
rjs(:,ms2))/h_mesh%global_diameter)**(2*
alpha)/(
inputs%sigma_min*
inputs%mu_min**2)
4435 wwiwwj = x * wwsi(ni)*wwsj(nj)
4441 y = normt * wwiwwj/muhj
4442 z = norm * muhi * wwiwwj
4444 hsij(1,ni,nj,ci,cj) = hsij(1,ni,nj,ci,cj) + y*normi(2)*normj(2) &
4445 + z*normi(1)*normj(1)
4446 hsij(4,ni,nj,ci,cj) = hsij(4,ni,nj,ci,cj) - y*normj(1)*normi(2) &
4447 + z*normi(1)*normj(2)
4449 hsij(7,ni,nj,ci,cj) = hsij(7,ni,nj,ci,cj) - y*normj(2)*normi(1) &
4450 + z*normi(2)*normj(1)
4452 hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) + y*(normi(1)*normj(1) + normi(2)*normj(2))
4453 hsij(6,ni,nj,ci,cj) = hsij(6,ni,nj,ci,cj) + y*normi(1)*normj(1) &
4454 + z*normi(2)*normj(2)
4467 i = interface_h_mu%jjs1(ni,ms)
4469 i = interface_h_mu%jjs2(ni,ms)
4471 ib = la_h%loc_to_glob(ki,i)
4472 ix = ni + n_ws1*((ki-1) + 3*(ci-1))
4479 j = interface_h_mu%jjs1(nj,ms)
4481 j = interface_h_mu%jjs2(nj,ms)
4483 jb = la_h%loc_to_glob(kj,j)
4484 jx = nj + n_ws2*((kj-1) + 3*(cj-1))
4486 IF ((ki == 1) .AND. (kj == 1))
THEN 4487 mat_loc1(ix,jx) = hsij(1,ni,nj,ci,cj)
4488 mat_loc2(ix,jx) = hsij(1,ni,nj,ci,cj)
4489 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 4490 mat_loc1(ix,jx) = hsij(4,ni,nj,ci,cj)
4491 mat_loc2(ix,jx) = hsij(4,ni,nj,ci,cj)
4492 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 4497 mat_loc1(ix,jx) = hsij(7,ni,nj,ci,cj)
4498 mat_loc2(ix,jx) = hsij(7,ni,nj,ci,cj)
4499 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 4500 mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)
4501 mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)
4502 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 4503 mat_loc1(ix,jx) = hsij(6,ni,nj,ci,cj)
4504 mat_loc2(ix,jx) = hsij(6,ni,nj,ci,cj)
4513 CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
4514 mat_loc1(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
4515 CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
4516 mat_loc2(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
4525 DO ls = 1, h_mesh%gauss%l_Gs
4528 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
4531 muhl1 = sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4532 muhl2 = sum(mu_h_field(h_mesh%jjs(:,ms2))*
wws(:,ls))
4533 drmuhl1 = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(1,:,ls,ms1))
4534 drmuhl2 = sum(mu_h_field(h_mesh%jj(:,m2))*
dw_s(1,:,ls,ms2))
4535 dzmuhl1 = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(2,:,ls,ms1))
4536 dzmuhl2 = sum(mu_h_field(h_mesh%jj(:,m2))*
dw_s(2,:,ls,ms2))
4538 IF (jj_v_to_h(h_mesh%jj(1,m1)) == -1)
THEN 4541 sigmal1 = sum(sigma_np(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4543 IF (jj_v_to_h(h_mesh%jj(1,m2)) == -1)
THEN 4546 sigmal2 = sum(sigma_np(h_mesh%jjs(:,ms2))*
wws(:,ls))
4601 y = x*wwsi(ni)*wwsj(nj)/(2*sigmaj*muhj)
4603 hsij(2,ni,nj,ci,cj) = hsij(2,ni,nj,ci,cj) + y * (-mode/ray)*normi(1)
4604 hsij(3,ni,nj,ci,cj) = hsij(3,ni,nj,ci,cj) + y * mode/ray *normi(1)
4607 hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) - y * normi(1) *((+1/ray)+(-1/muhj)*drmuhj) &
4608 + y * normi(2)*(1/muhj)*dzmuhj
4610 hsij(8,ni,nj,ci,cj) = hsij(8,ni,nj,ci,cj) + y * (-mode/ray)*normi(2)
4611 hsij(9,ni,nj,ci,cj) = hsij(9,ni,nj,ci,cj) + y * mode/ray *normi(2)
4613 hsij(1,ni,nj,ci,cj) = hsij(1,ni,nj,ci,cj) - y * normi(2)*(-1/muhj)*dzmuhj
4614 hsij(4,ni,nj,ci,cj) = hsij(4,ni,nj,ci,cj) - y * normi(2)* (1/muhj)*drmuhj
4615 hsij(7,ni,nj,ci,cj) = hsij(7,ni,nj,ci,cj) + y * normi(1)*(-1/muhj)*dzmuhj
4616 hsij(6,ni,nj,ci,cj) = hsij(6,ni,nj,ci,cj) + y * normi(1)* (1/muhj)*drmuhj
4619 y = x*wwsi(ni)*wwsj(nj)/(2*sigmai)
4620 gsij(2,ni,nj,ci,cj) = gsij(2,ni,nj,ci,cj) + y * (-mode/ray)*normj(1)
4621 gsij(3,ni,nj,ci,cj) = gsij(3,ni,nj,ci,cj) + y * ( mode/ray)*normj(1)
4622 gsij(5,ni,nj,ci,cj) = gsij(5,ni,nj,ci,cj) + y * (-1/ray) *normj(1)
4623 gsij(8,ni,nj,ci,cj) = gsij(8,ni,nj,ci,cj) + y * (-mode/ray)*normj(2)
4624 gsij(9,ni,nj,ci,cj) = gsij(9,ni,nj,ci,cj) + y * mode/ray *normj(2)
4642 i = interface_h_mu%jjs1(ni,ms)
4644 i = interface_h_mu%jjs2(ni,ms)
4646 ib = la_h%loc_to_glob(ki,i)
4647 ix = ni + n_wsi*((ki-1) + 3*(ci-1))
4653 j = interface_h_mu%jjs1(nj,ms)
4655 j = interface_h_mu%jjs2(nj,ms)
4657 jb = la_h%loc_to_glob(kj,j)
4658 jx = nj + n_wsj*((kj-1) + 3*(cj-1))
4660 IF ((ki == 1) .AND. (kj == 1))
THEN 4661 mat_loc1(ix,jx) = hsij(1,ni,nj,ci,cj)
4662 mat_loc2(ix,jx) = hsij(1,ni,nj,ci,cj)
4663 ELSE IF((ki == 1) .AND. (kj == 2))
THEN 4664 mat_loc1(ix,jx) = gsij(2,ni,nj,ci,cj)
4665 mat_loc2(ix,jx) = gsij(3,ni,nj,ci,cj)
4666 ELSE IF ((ki == 1) .AND. (kj == 3))
THEN 4667 mat_loc1(ix,jx) = hsij(4,ni,nj,ci,cj)
4668 mat_loc2(ix,jx) = hsij(4,ni,nj,ci,cj)
4669 ELSE IF ((ki == 2) .AND. (kj == 1))
THEN 4670 mat_loc1(ix,jx) = hsij(2,ni,nj,ci,cj)
4671 mat_loc2(ix,jx) = hsij(3,ni,nj,ci,cj)
4672 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 4673 mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)+gsij(5,ni,nj,ci,cj)
4674 mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)+gsij(5,ni,nj,ci,cj)
4675 ELSEIF ((ki == 2) .AND. (kj == 3))
THEN 4676 mat_loc1(ix,jx) = hsij(8,ni,nj,ci,cj)
4677 mat_loc2(ix,jx) = hsij(9,ni,nj,ci,cj)
4678 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 4679 mat_loc1(ix,jx) = hsij(7,ni,nj,ci,cj)
4680 mat_loc2(ix,jx) = hsij(7,ni,nj,ci,cj)
4681 ELSEIF ((ki == 3) .AND. (kj == 2))
THEN 4682 mat_loc1(ix,jx) = gsij(8,ni,nj,ci,cj)
4683 mat_loc2(ix,jx) = gsij(9,ni,nj,ci,cj)
4684 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 4685 mat_loc1(ix,jx) = hsij(6,ni,nj,ci,cj)
4686 mat_loc2(ix,jx) = hsij(6,ni,nj,ci,cj)
4695 CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
4696 mat_loc1(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
4697 CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
4698 mat_loc2(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
4703 DO ls = 1, h_mesh%gauss%l_Gs
4706 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
4709 muhl1 = sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4710 muhl2 = sum(mu_h_field(h_mesh%jjs(:,ms2))*
wws(:,ls))
4713 IF (jj_v_to_h(h_mesh%jj(1,m1)) == -1)
THEN 4716 sigmal1 = sum(sigma_np(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4718 IF (jj_v_to_h(h_mesh%jj(1,m2)) == -1)
THEN 4721 sigmal2 = sum(sigma_np(h_mesh%jjs(:,ms2))*
wws(:,ls))
4728 dwsi = dw_cs(:,:,ls,ms1)
4738 dwsi =
dw_s(:,:,ls,ms2)
4750 dwsj = dw_cs(:,:,ls,ms1)
4761 dwsj =
dw_s(:,:,ls,ms2)
4776 y = x*wwsi(ni)/(2*sigmaj*muhj)
4778 hsij(1,ni,nj,ci,cj) = hsij(1,ni,nj,ci,cj) + y*(-dwsj(2,nj))*normi(2)
4779 hsij(4,ni,nj,ci,cj) = hsij(4,ni,nj,ci,cj) + y* dwsj(1,nj) *normi(2)
4780 hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) + y*(-dwsj(2,nj) *normi(2)-dwsj(1,nj)*normi(1))
4781 hsij(6,ni,nj,ci,cj) = hsij(6,ni,nj,ci,cj) + y*(-dwsj(1,nj))*normi(1)
4782 hsij(7,ni,nj,ci,cj) = hsij(7,ni,nj,ci,cj) + y* dwsj(2,nj) *normi(1)
4787 y = x*wwsj(nj)/(2*sigmai)
4788 gsij(1,ni,nj,ci,cj) = gsij(1,ni,nj,ci,cj) + y*(-dwsi(2,ni))*normj(2)
4789 gsij(4,ni,nj,ci,cj) = gsij(4,ni,nj,ci,cj) + y* dwsi(2,ni) *normj(1)
4790 gsij(5,ni,nj,ci,cj) = gsij(5,ni,nj,ci,cj) + y*(-dwsi(2,ni) *normj(2)-dwsi(1,ni)*normj(1))
4791 gsij(6,ni,nj,ci,cj) = gsij(6,ni,nj,ci,cj) + y*(-dwsi(1,ni))*normj(1)
4792 gsij(7,ni,nj,ci,cj) = gsij(7,ni,nj,ci,cj) + y* dwsi(1,ni) *normj(2)
4810 i = interface_h_mu%jjs1(ni,ms)
4812 i = interface_h_mu%jjs2(ni,ms)
4814 ib = la_h%loc_to_glob(ki,i)
4815 ix = ni + n_wsi*((ki-1) + 3*(ci-1))
4821 j = h_mesh%jj(nj,m1)
4823 j = h_mesh%jj(nj,m2)
4825 jb = la_h%loc_to_glob(kj,j)
4826 jx = nj + n_wj*((kj-1) + 3*(cj-1))
4828 IF ((ki == 1) .AND. (kj == 1))
THEN 4829 mat_loc1(ix,jx) = hsij(1,ni,nj,ci,cj)
4830 mat_loc2(ix,jx) = hsij(1,ni,nj,ci,cj)
4831 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 4832 mat_loc1(ix,jx) = hsij(4,ni,nj,ci,cj)
4833 mat_loc2(ix,jx) = hsij(4,ni,nj,ci,cj)
4834 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 4835 mat_loc1(ix,jx) = hsij(7,ni,nj,ci,cj)
4836 mat_loc2(ix,jx) = hsij(7,ni,nj,ci,cj)
4837 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 4838 mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)
4839 mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)
4840 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 4841 mat_loc1(ix,jx) = hsij(6,ni,nj,ci,cj)
4842 mat_loc2(ix,jx) = hsij(6,ni,nj,ci,cj)
4852 CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_w2, jdxn(1:6*n_w2), &
4853 mat_loc1(1:6*n_ws1,1:6*n_w2), add_values, ierr)
4854 CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_w2, jdxn(1:6*n_w2), &
4855 mat_loc2(1:6*n_ws1,1:6*n_w2), add_values, ierr)
4863 i = h_mesh%jj(ni,m1)
4865 i = h_mesh%jj(ni,m2)
4867 ib = la_h%loc_to_glob(ki,i)
4868 ix = ni + n_wi*((ki-1) + 3*(ci-1))
4874 j = interface_h_mu%jjs1(nj,ms)
4876 j = interface_h_mu%jjs2(nj,ms)
4878 jb = la_h%loc_to_glob(kj,j)
4879 jx = nj + n_wsj*((kj-1) + 3*(cj-1))
4881 IF ((ki == 1) .AND. (kj == 1))
THEN 4882 mat_loc1(ix,jx) = gsij(1,ni,nj,ci,cj)
4883 mat_loc2(ix,jx) = gsij(1,ni,nj,ci,cj)
4884 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN 4885 mat_loc1(ix,jx) = gsij(4,ni,nj,ci,cj)
4886 mat_loc2(ix,jx) = gsij(4,ni,nj,ci,cj)
4887 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN 4888 mat_loc1(ix,jx) = gsij(7,ni,nj,ci,cj)
4889 mat_loc2(ix,jx) = gsij(7,ni,nj,ci,cj)
4890 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN 4891 mat_loc1(ix,jx) = gsij(5,ni,nj,ci,cj)
4892 mat_loc2(ix,jx) = gsij(5,ni,nj,ci,cj)
4893 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN 4894 mat_loc1(ix,jx) = gsij(6,ni,nj,ci,cj)
4895 mat_loc2(ix,jx) = gsij(6,ni,nj,ci,cj)
4904 CALL matsetvalues(h_p_phi_mat1, 6*n_w1, idxn(1:6*n_w1), 6*n_ws2, jdxn(1:6*n_ws2), &
4905 mat_loc1(1:6*n_w1,1:6*n_ws2), add_values, ierr)
4906 CALL matsetvalues(h_p_phi_mat2, 6*n_w1, idxn(1:6*n_w1), 6*n_ws2, jdxn(1:6*n_ws2), &
4907 mat_loc2(1:6*n_w1,1:6*n_ws2), add_values, ierr)
4911 CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
4912 CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
4913 CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
4914 CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
4921 SUBROUTINE courant_mu(H_mesh,interface_H_mu,sigma,mu_H_field,time,mode,nl, &
4922 la_h, vb_1, vb_2, b_ext, j_over_sigma_gauss, sigma_curl_gauss)
4929 #include "petsc/finclude/petsc.h" 4934 REAL(KIND=8),
INTENT(IN) :: time
4935 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN) :: sigma
4936 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
4937 INTEGER,
INTENT(IN) :: mode
4938 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: nl
4939 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: B_ext
4940 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: J_over_sigma_gauss
4941 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: sigma_curl_gauss
4942 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: src_H
4943 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_ws,H_mesh%gauss%l_Gs) :: w_cs
4944 REAL(KIND=8),
DIMENSION(2) :: normi, gaussp1, gaussp2
4945 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_ws) :: wwsi
4946 REAL(KIND=8) :: x, ray
4947 INTEGER :: i, ni, ms, k, ls, n_ws1, n_ws2, ms1, ms2, n_w1, n_w2, m1, m2, ci, n_wsi
4948 INTEGER :: mesh_id1, mesh_id2, index
4949 REAL(KIND=8),
DIMENSION(6) :: JsolH_anal, test, B_ext_l, J_over_sigma_l
4950 REAL(KIND=8) :: muhl1, muhl2, ref, diff
4957 INTEGER,
DIMENSION(H_mesh%np) :: idxn
4960 petscerrorcode :: ierr
4970 n_ws1 = h_mesh%gauss%n_ws
4971 n_ws2 = h_mesh%gauss%n_ws
4972 n_w1 = h_mesh%gauss%n_w
4973 n_w2 = h_mesh%gauss%n_w
4975 DO ms = 1, interface_h_mu%mes
4976 ms1 = interface_h_mu%mesh1(ms)
4977 ms2 = interface_h_mu%mesh2(ms)
4978 m1 = h_mesh%neighs(ms1)
4979 m2 = h_mesh%neighs(ms2)
4981 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
4982 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(1,ms2)))**2)
4983 IF (diff/ref .LT. 1.d-10)
THEN 4987 w_cs(1,ls)=
wws(2,ls)
4988 w_cs(2,ls)=
wws(1,ls)
4989 IF (n_ws1==3) w_cs(n_ws1,ls) =
wws(n_ws1,ls)
4990 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 4996 DO ms = 1, interface_h_mu%mes
4997 ms2 = interface_h_mu%mesh2(ms)
4998 ms1 = interface_h_mu%mesh1(ms)
4999 m2 = h_mesh%neighs(ms2)
5000 m1 = h_mesh%neighs(ms1)
5001 mesh_id1 = h_mesh%i_d(m1)
5002 mesh_id2 = h_mesh%i_d(m2)
5005 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
5010 b_ext_l(k) = sum(b_ext(h_mesh%jjs(:,ms1),k)*h_mesh%gauss%wws(:,ls))
5011 j_over_sigma_l(k) = j_over_sigma_gauss(index,k) + sigma_curl_gauss(index,k)
5013 muhl1=sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
5014 gaussp1(1) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms1))*w_cs(:,ls))
5015 gaussp1(2) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms1))*w_cs(:,ls))
5017 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 5019 jsolh_anal(k) = j_over_sigma_l(k) &
5020 + sum(nl(h_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
5024 jsolh_anal(k) = jexact_gauss(k, gaussp1, mode, one ,sigma(m1), &
5025 muhl1, time, mesh_id1, b_ext_l)/sigma(m1) &
5026 + sum(nl(h_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
5033 b_ext_l(k) = sum(b_ext(h_mesh%jjs(:,ms2),k)*h_mesh%gauss%wws(:,ls))
5034 j_over_sigma_l(k) = j_over_sigma_gauss(index,k) + sigma_curl_gauss(index,k)
5036 muhl2=sum(mu_h_field(h_mesh%jjs(:,ms2))*
wws(:,ls))
5037 gaussp2(1) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))*
wws(:,ls))
5038 gaussp2(2) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms2))*
wws(:,ls))
5039 IF (maxval(abs(gaussp1-gaussp2)) > 1.d-11)
THEN 5040 WRITE(*,*)
' BUG courant_mu ' 5044 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 5046 test(k) = j_over_sigma_l(k) &
5047 + sum(nl(h_mesh%jjs(1:n_ws2,ms2),k)*
wws(1:n_ws2,ls))
5048 jsolh_anal(k) = jsolh_anal(k) + test(k)
5052 test(k) = jexact_gauss(k, gaussp2, mode, one ,sigma(m2), &
5053 muhl2, time, mesh_id2, b_ext_l)/sigma(m2) &
5054 + sum(nl(h_mesh%jjs(1:n_ws2,ms2),k)*
wws(1:n_ws2,ls))
5055 jsolh_anal(k) = jsolh_anal(k) + test(k)
5073 i = interface_h_mu%jjs1(ni,ms)
5075 i = interface_h_mu%jjs2(ni,ms)
5077 x =
rjs(ls,ms2)*ray*wwsi(ni)/2
5078 src_h(i,1) = src_h(i,1)+x*(-jsolh_anal(3)*normi(2))
5079 src_h(i,2) = src_h(i,2)+x*(-jsolh_anal(4)*normi(2))
5080 src_h(i,3) = src_h(i,3)+x*(jsolh_anal(1)*normi(2)-jsolh_anal(5)*normi(1))
5081 src_h(i,4) = src_h(i,4)+x*(jsolh_anal(2)*normi(2)-jsolh_anal(6)*normi(1))
5082 src_h(i,5) = src_h(i,5)+x*(jsolh_anal(3)*normi(1))
5083 src_h(i,6) = src_h(i,6)+x*(jsolh_anal(4)*normi(1))
5089 IF (h_mesh%np /= 0)
THEN 5091 idxn = la_h%loc_to_glob(1,:)-1
5092 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,1), add_values, ierr)
5093 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,2), add_values, ierr)
5094 idxn = la_h%loc_to_glob(2,:)-1
5095 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,4), add_values, ierr)
5096 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,3), add_values, ierr)
5097 idxn = la_h%loc_to_glob(3,:)-1
5098 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,5), add_values, ierr)
5099 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,6), add_values, ierr)
5102 CALL vecassemblybegin(vb_1,ierr)
5103 CALL vecassemblyend(vb_1,ierr)
5104 CALL vecassemblybegin(vb_2,ierr)
5105 CALL vecassemblyend(vb_2,ierr)
5109 SUBROUTINE rhs_dirichlet(H_mesh,Dirichlet_bdy_H_sides,sigma,&
5110 mu_h_field,time,mode,nl,stab, la_h, vb_1, vb_2, b_ext, j_over_sigma_bdy, &
5117 #include "petsc/finclude/petsc.h" 5121 INTEGER,
DIMENSION(:),
INTENT(IN) :: Dirichlet_bdy_H_sides
5122 REAL(KIND=8),
INTENT(IN) :: time
5123 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN) :: sigma
5124 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
5125 INTEGER,
INTENT(IN) :: mode
5126 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: nl
5127 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
5128 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: B_ext
5129 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: J_over_sigma_bdy
5130 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: sigma_curl_bdy_in
5131 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: src_H
5132 REAL(KIND=8),
DIMENSION(2) :: gaussp1
5133 REAL(KIND=8) :: x, ray, stab_colle_H_mu
5134 INTEGER :: i, ni, ms, k, ls, m1, count
5136 REAL(KIND=8),
DIMENSION(6) :: JsolH_anal, B_ext_l
5137 REAL(KIND=8) :: muhl1, hm1
5138 REAL(KIND=8),
DIMENSION(6,H_mesh%gauss%l_Gs) :: Hloc, Hlocxn
5139 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%l_Gs) :: rloc
5140 REAL(KIND=8),
DIMENSION(1) :: muloc
5148 INTEGER,
DIMENSION(H_mesh%np) :: idxn
5151 petscerrorcode :: ierr
5165 stab_colle_h_mu = stab(3)
5168 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
5169 ms = dirichlet_bdy_h_sides(count)
5172 hm1 = stab_colle_h_mu/(sum(h_mesh%gauss%rjs(:,ms))*
inputs%sigma_min)
5174 m1 = h_mesh%neighs(ms)
5175 mesh_id1 = h_mesh%i_d(m1)
5176 muloc(1) = mu_h_field(h_mesh%jj(1,m1))
5177 DO ls = 1, h_mesh%gauss%l_Gs
5178 rloc(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
5179 rloc(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
5183 hloc(k,:) = hexact(h_mesh, k, rloc, mode, muloc, time)
5186 hlocxn(1,:) = hloc(3,:)*h_mesh%gauss%rnorms(2,:,ms)
5187 hlocxn(2,:) = hloc(4,:)*h_mesh%gauss%rnorms(2,:,ms)
5188 hlocxn(3,:) = hloc(5,:)*h_mesh%gauss%rnorms(1,:,ms)-hloc(1,:)*h_mesh%gauss%rnorms(2,:,ms)
5189 hlocxn(4,:) = hloc(6,:)*h_mesh%gauss%rnorms(1,:,ms)-hloc(2,:)*h_mesh%gauss%rnorms(2,:,ms)
5190 hlocxn(5,:) = -hloc(3,:)*h_mesh%gauss%rnorms(1,:,ms)
5191 hlocxn(6,:) = -hloc(4,:)*h_mesh%gauss%rnorms(1,:,ms)
5193 DO ls = 1, h_mesh%gauss%l_Gs
5198 b_ext_l(k) = sum(b_ext(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls))
5202 muhl1=sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
5203 gaussp1(1) = rloc(1,ls)
5204 gaussp1(2) = rloc(2,ls)
5212 IF (
inputs%if_level_set.AND.
inputs%variation_sigma_fluid)
THEN 5214 jsolh_anal(k) = j_over_sigma_bdy(index,k) &
5215 + sigma_curl_bdy_in(index,k) &
5216 + sum(nl(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls)) &
5221 jsolh_anal(k) = jexact_gauss(k, gaussp1, mode, one ,sigma(m1), muhl1, &
5222 time, mesh_id1, b_ext_l)/sigma(m1) &
5223 + sum(nl(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls)) &
5228 DO ni = 1, h_mesh%gauss%n_ws
5229 i = h_mesh%jjs(ni,ms)
5230 x = h_mesh%gauss%rjs(ls,ms)*ray*h_mesh%gauss%wws(ni,ls)
5232 src_h(i,1) = src_h(i,1)+x*(-jsolh_anal(3)*h_mesh%gauss%rnorms(2,ls,ms))
5233 src_h(i,2) = src_h(i,2)+x*(-jsolh_anal(4)*h_mesh%gauss%rnorms(2,ls,ms))
5234 src_h(i,3) = src_h(i,3)+x*(jsolh_anal(1)*h_mesh%gauss%rnorms(2,ls,ms)&
5235 -jsolh_anal(5)*h_mesh%gauss%rnorms(1,ls,ms))
5236 src_h(i,4) = src_h(i,4)+x*(jsolh_anal(2)*h_mesh%gauss%rnorms(2,ls,ms)&
5237 -jsolh_anal(6)*h_mesh%gauss%rnorms(1,ls,ms))
5238 src_h(i,5) = src_h(i,5)+x*(jsolh_anal(3)*h_mesh%gauss%rnorms(1,ls,ms))
5239 src_h(i,6) = src_h(i,6)+x*(jsolh_anal(4)*h_mesh%gauss%rnorms(1,ls,ms))
5246 IF (h_mesh%np /= 0)
THEN 5248 idxn = la_h%loc_to_glob(1,:)-1
5249 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,1), add_values, ierr)
5250 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,2), add_values, ierr)
5251 idxn = la_h%loc_to_glob(2,:)-1
5252 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,4), add_values, ierr)
5253 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,3), add_values, ierr)
5254 idxn = la_h%loc_to_glob(3,:)-1
5255 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,5), add_values, ierr)
5256 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,6), add_values, ierr)
5259 CALL vecassemblybegin(vb_1,ierr)
5260 CALL vecassemblyend(vb_1,ierr)
5261 CALL vecassemblybegin(vb_2,ierr)
5262 CALL vecassemblyend(vb_2,ierr)
5271 #include "petsc/finclude/petsc.h" 5276 INTEGER,
POINTER,
DIMENSION(:) :: js_D
5277 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: on_proc_loc, on_proc, not_cav_loc, not_cav
5278 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: is_ok, j_tmp
5279 INTEGER,
DIMENSION(1) :: loc
5280 INTEGER :: m, ms, i, nb_dom, idx, nb_cav, ni
5282 #include "petsc/finclude/petsc.h" 5283 mpi_comm,
INTENT(IN) :: communicator
5285 petscerrorcode :: ierr
5287 IF (
inputs%nb_dom_phi==0)
RETURN 5289 CALL mpi_comm_rank(communicator, rank, ierr)
5291 nb_dom =
inputs%nb_dom_phi
5292 ALLOCATE(on_proc_loc(nb_dom), on_proc(nb_dom))
5293 ALLOCATE(not_cav_loc(nb_dom), not_cav(nb_dom))
5301 IF (minval(abs(
inputs%list_dom_phi-i)) /= 0)
THEN 5302 WRITE(*,*)
'error in dirichlet cavities' 5304 loc = minloc(abs(
inputs%list_dom_phi-i))
5305 on_proc_loc(loc(1)) = rank
5307 IF (mesh%mes /= 0)
THEN 5308 ALLOCATE(is_ok(mesh%mes))
5309 is_ok = mesh%i_d(mesh%neighs)
5310 IF (interface_h_phi%mes /=0)
THEN 5311 is_ok(interface_h_phi%mesh2) = 0
5314 IF (sum(abs(mesh%rr(1,mesh%jjs(:,ms)))) .LT. 1.d-12*mesh%global_diameter)
THEN 5317 IF (
inputs%my_periodic%nb_periodic_pairs /=0)
THEN 5318 IF (minval(abs(
inputs%my_periodic%list_periodic-mesh%sides(ms))) == 0)
THEN 5325 IF (is_ok(ms) == 0) cycle
5327 IF (minval(abs(
inputs%list_dom_phi-i)) /= 0)
THEN 5328 WRITE(*,*)
'error in dirichlet cavities' 5330 loc = minloc(abs(
inputs%list_dom_phi-i))
5331 not_cav_loc(loc(1)) = rank
5334 CALL mpi_allreduce(on_proc_loc, on_proc, nb_dom, mpi_integer, mpi_max, communicator, ierr)
5335 CALL mpi_allreduce(not_cav_loc, not_cav, nb_dom, mpi_integer, mpi_max, communicator, ierr)
5337 ALLOCATE(j_tmp(
SIZE(js_d)+nb_dom))
5338 j_tmp(1:
SIZE(js_d)) = js_d
5341 IF ( (not_cav(i)==-1) .AND. (on_proc(i)==rank) )
THEN 5345 IF (mesh%i_d(m) ==
inputs%list_dom_phi(i))
THEN 5346 DO ni = 1, mesh%gauss%n_w
5347 IF (minval(abs(mesh%jjs-mesh%jj(ni,m))) /=0)
THEN 5348 j_tmp(idx) = mesh%jj(ni,m)
5354 WRITE(*,*)
'add ', j_tmp(idx),
'in dom ',
inputs%list_dom_phi(i),
' : proc ', rank
5355 WRITE(*,*)
'add ', mesh%rr(:,j_tmp(idx)), mesh%i_d(m)
5363 nb_cav = idx -
SIZE(js_d)
5364 IF (nb_cav /= 0)
THEN 5370 DEALLOCATE(on_proc_loc, on_proc, j_tmp)
5371 DEALLOCATE(not_cav_loc, not_cav)
5372 IF (
ALLOCATED(is_ok))
DEALLOCATE(is_ok)
5374 WRITE(*,
'(a,x,i2,x,a,x,i2)')
'I have detected', nb_cav,
' cavity(ies) on proc', rank
5384 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: angles
5385 INTEGER,
INTENT(IN) :: nb_angles
5386 INTEGER,
INTENT(IN) :: nb, ne
5387 REAL(KIND=8),
INTENT(IN) :: time
5388 REAL(KIND=8),
DIMENSION(nb_angles,ne-nb+1) :: vv
5390 vv = -1/mu_in_real_space(h_mesh,angles,nb_angles,nb,ne,time)
5395 FUNCTION one_over_mu(H_mesh,angles,nb_angles,nb,ne,time) RESULT(vv)
5400 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: angles
5401 INTEGER,
INTENT(IN) :: nb_angles
5402 INTEGER,
INTENT(IN) :: nb, ne
5403 REAL(KIND=8),
INTENT(IN) :: time
5404 REAL(KIND=8),
DIMENSION(nb_angles,ne-nb+1) :: vv
5406 vv = 1/mu_in_real_space(h_mesh,angles,nb_angles,nb,ne,time)
5410 SUBROUTINE smb_sigma_prod_curl(communicator, mesh, jj_v_to_H, list_mode, H_in, one_over_sigma_in, sigma_nj_m,&
5418 #include "petsc/finclude/petsc.h" 5422 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
5423 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5424 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: H_in
5425 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_in
5426 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%me),
INTENT(IN) :: sigma_nj_m
5427 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
5428 REAL(KIND=8),
DIMENSION(:,:,:) :: V_out
5429 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: H_gauss, RotH
5430 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: one_over_sigma_gauss
5431 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: RotH_bar
5432 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
5433 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
5434 INTEGER :: m, l , i, mode, index, k
5435 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,6) :: H_in_loc
5436 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: one_over_sigma_in_loc
5437 REAL(KIND=8) :: ray, sigma_np_gauss
5438 INTEGER :: nb_procs, bloc_size, m_max_pad, code
5439 mpi_comm :: communicator
5441 DO i = 1,
SIZE(list_mode)
5445 j_loc = mesh%jj(:,m)
5447 h_in_loc(:,k) = h_in(j_loc,k,i)
5450 one_over_sigma_in_loc(:,k) = one_over_sigma_in(j_loc,k,i)
5453 DO l = 1, mesh%gauss%l_G
5455 dw_loc = mesh%gauss%dw(:,:,l,m)
5458 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
5461 h_gauss(index,1,i) = sum(h_in_loc(:,1)*mesh%gauss%ww(:,l))
5462 h_gauss(index,3,i) = sum(h_in_loc(:,3)*mesh%gauss%ww(:,l))
5463 h_gauss(index,5,i) = sum(h_in_loc(:,5)*mesh%gauss%ww(:,l))
5465 h_gauss(index,2,i) = sum(h_in_loc(:,2)*mesh%gauss%ww(:,l))
5466 h_gauss(index,4,i) = sum(h_in_loc(:,4)*mesh%gauss%ww(:,l))
5467 h_gauss(index,6,i) = sum(h_in_loc(:,6)*mesh%gauss%ww(:,l))
5470 roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
5471 -sum(h_in_loc(:,3)*dw_loc(2,:))
5472 roth(index,4,i) = sum(h_in_loc(:,2)*dw_loc(2,:)) &
5473 -sum(h_in_loc(:,6)*dw_loc(1,:))
5474 roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
5475 +sum(h_in_loc(:,3)*dw_loc(1,:)) &
5476 -mode/ray*h_gauss(index,2,i)
5478 roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
5479 -sum(h_in_loc(:,4)*dw_loc(2,:))
5480 roth(index,3,i) = sum(h_in_loc(:,1)*dw_loc(2,:)) &
5481 -sum(h_in_loc(:,5)*dw_loc(1,:))
5482 roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
5483 +sum(h_in_loc(:,4)*dw_loc(1,:))&
5484 +mode/ray*h_gauss(index,1,i)
5486 IF (jj_v_to_h(mesh%jj(1,m)) == -1)
THEN 5487 one_over_sigma_gauss(index,1,i) = 0.d0
5488 one_over_sigma_gauss(index,2,i) = 0.d0
5490 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m)
5493 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc(:,1)*mesh%gauss%ww(:,l))
5494 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc(:,2)*mesh%gauss%ww(:,l))
5497 sigma_np_gauss = sum(sigma_nj_m(:,m)*mesh%gauss%ww(:,l))
5499 roth_bar(index,k,i) = roth(index,k,i)/sigma_np_gauss
5505 CALL mpi_comm_size(communicator, nb_procs, code)
5506 bloc_size =
SIZE(one_over_sigma_gauss,1)/nb_procs+1
5507 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5509 roth, one_over_sigma_gauss, v_out, 1, nb_procs, bloc_size, m_max_pad)
5510 v_out = roth_bar - v_out
5515 h_in, one_over_sigma_in, sigma_np, sigma, v_out)
5522 #include "petsc/finclude/petsc.h" 5526 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
5527 INTEGER,
DIMENSION(:),
INTENT(IN) :: Dirichlet_bdy_H_sides
5528 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5529 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: H_in
5530 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_in
5531 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_np
5532 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
5533 REAL(KIND=8),
DIMENSION(:,:,:) :: V_out
5534 REAL(KIND=8),
DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),6,SIZE(list_mode)) :: H_gauss, RotH
5535 REAL(KIND=8),
DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),6,SIZE(list_mode)) :: RotH_bar
5536 REAL(KIND=8),
DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),2,SIZE(list_mode)) :: one_over_sigma_gauss
5537 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
5538 INTEGER :: ms, ls , i, mode, index, k
5539 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc
5540 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,6) :: H_in_loc
5541 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: one_over_sigma_in_loc
5543 INTEGER :: nb_procs, bloc_size, m_max_pad, code, count, m1
5544 mpi_comm :: communicator
5546 DO i = 1,
SIZE(list_mode)
5549 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
5550 ms = dirichlet_bdy_h_sides(count)
5551 m1 = mesh%neighs(ms)
5553 j_loc = mesh%jjs(:,ms)
5555 h_in_loc(:,k) = h_in(j_loc,k,i)
5558 one_over_sigma_in_loc(:,k) = one_over_sigma_in(j_loc,k,i)
5561 DO ls = 1, mesh%gauss%l_Gs
5563 dw_loc = mesh%gauss%dw_s(:,:,ls,ms)
5566 ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
5569 h_gauss(index,1,i) = sum(h_in_loc(:,1)*mesh%gauss%wws(:,ls))
5570 h_gauss(index,3,i) = sum(h_in_loc(:,3)*mesh%gauss%wws(:,ls))
5571 h_gauss(index,5,i) = sum(h_in_loc(:,5)*mesh%gauss%wws(:,ls))
5573 h_gauss(index,2,i) = sum(h_in_loc(:,2)*mesh%gauss%wws(:,ls))
5574 h_gauss(index,4,i) = sum(h_in_loc(:,4)*mesh%gauss%wws(:,ls))
5575 h_gauss(index,6,i) = sum(h_in_loc(:,6)*mesh%gauss%wws(:,ls))
5578 roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
5579 -sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(2,:))
5581 roth(index,4,i) = sum(h_in(mesh%jj(:,m1),2,i)*dw_loc(2,:)) &
5582 -sum(h_in(mesh%jj(:,m1),6,i)*dw_loc(1,:))
5584 roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
5585 +sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(1,:)) &
5586 -mode/ray*h_gauss(index,2,i)
5589 roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
5590 -sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(2,:))
5592 roth(index,3,i) = sum(h_in(mesh%jj(:,m1),1,i)*dw_loc(2,:)) &
5593 -sum(h_in(mesh%jj(:,m1),5,i)*dw_loc(1,:))
5595 roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
5596 +sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(1,:))&
5597 +mode/ray*h_gauss(index,1,i)
5599 IF (jj_v_to_h(mesh%jj(1,m1)) == -1)
THEN 5600 one_over_sigma_gauss(index,1,i) = 0.d0
5601 one_over_sigma_gauss(index,2,i) = 0.d0
5603 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m1)
5606 roth_bar(index,k,i) = roth(index,k,i)/sigma(m1)
5609 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc(:,1)*mesh%gauss%wws(:,ls))
5610 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc(:,2)*mesh%gauss%wws(:,ls))
5612 roth_bar(index,k,i) = roth(index,k,i)/sum(sigma_np(mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
5619 IF (
SIZE(dirichlet_bdy_h_sides).GE.1)
THEN 5620 CALL mpi_comm_size(communicator, nb_procs, code)
5621 bloc_size =
SIZE(one_over_sigma_gauss,1)/nb_procs+1
5622 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5624 roth, one_over_sigma_gauss, v_out, 1, nb_procs, bloc_size, m_max_pad)
5626 v_out = roth_bar - v_out
5632 h_in, one_over_sigma_in, sigma_np, sigma, v_out)
5639 #include "petsc/finclude/petsc.h" 5643 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
5645 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5646 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: H_in
5647 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_in
5648 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_np
5649 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
5650 REAL(KIND=8),
DIMENSION(:,:,:) :: V_out
5651 REAL(KIND=8),
DIMENSION(2*mesh%gauss%l_Gs*interface_H_mu%mes,6,SIZE(list_mode)) :: H_gauss, RotH
5652 REAL(KIND=8),
DIMENSION(2*mesh%gauss%l_Gs*interface_H_mu%mes,6,SIZE(list_mode)) :: RotH_bar
5653 REAL(KIND=8),
DIMENSION(2*mesh%gauss%l_Gs*interface_H_mu%mes,2,SIZE(list_mode)) :: one_over_sigma_gauss
5654 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
5655 INTEGER :: ms, ls , i, mode, index, k
5656 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc1, j_loc2
5657 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,6) :: H_in_loc1, H_in_loc2
5658 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: one_over_sigma_in_loc1, one_over_sigma_in_loc2
5659 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,mesh%gauss%l_Gs) :: w_cs
5660 REAL(KIND=8) :: ray, diff, ref
5661 INTEGER :: nb_procs, bloc_size, m_max_pad, code
5662 INTEGER :: ms1, ms2, m1, m2, mesh_id1, mesh_id2
5663 mpi_comm :: communicator
5665 DO ms = 1, interface_h_mu%mes
5666 ms1 = interface_h_mu%mesh1(ms)
5667 ms2 = interface_h_mu%mesh2(ms)
5668 m1 = mesh%neighs(ms1)
5669 m2 = mesh%neighs(ms2)
5670 ref = 1.d-8+sum((mesh%rr(:,mesh%jjs(1,ms1)) - mesh%rr(:,mesh%jjs(2,ms1)))**2)
5671 diff = sum((mesh%rr(:,mesh%jjs(1,ms1)) - mesh%rr(:,mesh%jjs(1,ms2)))**2)
5672 IF (diff/ref .LT. 1.d-10)
THEN 5673 w_cs = mesh%gauss%wws
5675 DO ls = 1, mesh%gauss%l_Gs
5676 w_cs(1,ls)= mesh%gauss%wws(2,ls)
5677 w_cs(2,ls)= mesh%gauss%wws(1,ls)
5678 IF (mesh%gauss%n_ws==3) w_cs(mesh%gauss%n_ws,ls) = mesh%gauss%wws(mesh%gauss%n_ws,ls)
5679 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 5684 DO i = 1,
SIZE(list_mode)
5687 DO ms = 1, interface_h_mu%mes
5688 ms2 = interface_h_mu%mesh2(ms)
5689 ms1 = interface_h_mu%mesh1(ms)
5690 m2 = mesh%neighs(ms2)
5691 m1 = mesh%neighs(ms1)
5692 mesh_id1 = mesh%i_d(m1)
5693 mesh_id2 = mesh%i_d(m2)
5694 j_loc1 = mesh%jjs(:,ms1)
5695 j_loc2 = mesh%jjs(:,ms2)
5697 h_in_loc1(:,k) = h_in(j_loc1,k,i)
5698 h_in_loc2(:,k) = h_in(j_loc2,k,i)
5701 one_over_sigma_in_loc1(:,k) = one_over_sigma_in(j_loc1,k,i)
5702 one_over_sigma_in_loc2(:,k) = one_over_sigma_in(j_loc2,k,i)
5705 DO ls = 1, mesh%gauss%l_Gs
5708 dw_loc = mesh%gauss%dw_s(:,:,ls,ms1)
5710 ray = sum(mesh%rr(1,mesh%jjs(:,ms1))*w_cs(:,ls))
5712 h_gauss(index,1,i) = sum(h_in_loc1(:,1)*w_cs(:,ls))
5713 h_gauss(index,3,i) = sum(h_in_loc1(:,3)*w_cs(:,ls))
5714 h_gauss(index,5,i) = sum(h_in_loc1(:,5)*w_cs(:,ls))
5715 h_gauss(index,2,i) = sum(h_in_loc1(:,2)*w_cs(:,ls))
5716 h_gauss(index,4,i) = sum(h_in_loc1(:,4)*w_cs(:,ls))
5717 h_gauss(index,6,i) = sum(h_in_loc1(:,6)*w_cs(:,ls))
5720 roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
5721 -sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(2,:))
5722 roth(index,4,i) = sum(h_in(mesh%jj(:,m1),2,i)*dw_loc(2,:)) &
5723 -sum(h_in(mesh%jj(:,m1),6,i)*dw_loc(1,:))
5724 roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
5725 +sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(1,:)) &
5726 -mode/ray*h_gauss(index,2,i)
5728 roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
5729 -sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(2,:))
5730 roth(index,3,i) = sum(h_in(mesh%jj(:,m1),1,i)*dw_loc(2,:)) &
5731 -sum(h_in(mesh%jj(:,m1),5,i)*dw_loc(1,:))
5732 roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
5733 +sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(1,:))&
5734 +mode/ray*h_gauss(index,1,i)
5736 IF (jj_v_to_h(mesh%jj(1,m1)) == -1)
THEN 5737 one_over_sigma_gauss(index,1,i) = 0.d0
5738 one_over_sigma_gauss(index,2,i) = 0.d0
5740 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m1)
5743 roth_bar(index,k,i) = roth(index,k,i)/sigma(m1)
5746 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc1(:,1)*w_cs(:,ls))
5747 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc1(:,2)*w_cs(:,ls))
5749 roth_bar(index,k,i) = roth(index,k,i)/sum(sigma_np(mesh%jjs(:,ms1))*w_cs(:,ls))
5755 dw_loc = mesh%gauss%dw_s(:,:,ls,ms2)
5757 ray = sum(mesh%rr(1,mesh%jjs(:,ms2))*mesh%gauss%wws(:,ls))
5759 h_gauss(index,1,i) = sum(h_in_loc2(:,1)*mesh%gauss%wws(:,ls))
5760 h_gauss(index,3,i) = sum(h_in_loc2(:,3)*mesh%gauss%wws(:,ls))
5761 h_gauss(index,5,i) = sum(h_in_loc2(:,5)*mesh%gauss%wws(:,ls))
5762 h_gauss(index,2,i) = sum(h_in_loc2(:,2)*mesh%gauss%wws(:,ls))
5763 h_gauss(index,4,i) = sum(h_in_loc2(:,4)*mesh%gauss%wws(:,ls))
5764 h_gauss(index,6,i) = sum(h_in_loc2(:,6)*mesh%gauss%wws(:,ls))
5767 roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
5768 -sum(h_in(mesh%jj(:,m2),3,i)*dw_loc(2,:))
5769 roth(index,4,i) = sum(h_in(mesh%jj(:,m2),2,i)*dw_loc(2,:)) &
5770 -sum(h_in(mesh%jj(:,m2),6,i)*dw_loc(1,:))
5771 roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
5772 +sum(h_in(mesh%jj(:,m2),3,i)*dw_loc(1,:)) &
5773 -mode/ray*h_gauss(index,2,i)
5775 roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
5776 -sum(h_in(mesh%jj(:,m2),4,i)*dw_loc(2,:))
5777 roth(index,3,i) = sum(h_in(mesh%jj(:,m2),1,i)*dw_loc(2,:)) &
5778 -sum(h_in(mesh%jj(:,m2),5,i)*dw_loc(1,:))
5779 roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
5780 +sum(h_in(mesh%jj(:,m2),4,i)*dw_loc(1,:))&
5781 +mode/ray*h_gauss(index,1,i)
5783 IF (jj_v_to_h(mesh%jj(1,m2)) == -1)
THEN 5784 one_over_sigma_gauss(index,1,i) = 0.d0
5785 one_over_sigma_gauss(index,2,i) = 0.d0
5787 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m2)
5790 roth_bar(index,k,i) = roth(index,k,i)/sigma(m2)
5793 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc2(:,1)*mesh%gauss%wws(:,ls))
5794 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc2(:,2)*mesh%gauss%wws(:,ls))
5796 roth_bar(index,k,i) = roth(index,k,i)/sum(sigma_np(mesh%jjs(:,ms2))*mesh%gauss%wws(:,ls))
5803 IF (interface_h_mu%mes.GE.1)
THEN 5804 CALL mpi_comm_size(communicator, nb_procs, code)
5805 bloc_size =
SIZE(one_over_sigma_gauss,1)/nb_procs+1
5806 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5808 roth, one_over_sigma_gauss, v_out, 1, nb_procs, bloc_size, m_max_pad)
5810 v_out = roth_bar - v_out
5816 mu_h_field, mu_phi, one_over_sigma_tot, time, sigma, j_over_sigma_gauss)
5822 #include "petsc/finclude/petsc.h" 5826 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
5827 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5828 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: B_in
5829 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_tot
5830 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
5831 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
5832 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
5833 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: J_over_sigma_gauss
5834 REAL(KIND=8),
DIMENSION(mesh%me*mesh%gauss%l_G,6,SIZE(list_mode)) :: J_exact_gauss
5835 REAL(KIND=8),
DIMENSION(mesh%me*mesh%gauss%l_G,2,SIZE(list_mode)) :: one_over_sigma_gauss
5836 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
5837 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,6) :: B_in_loc
5838 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: one_over_sigma_in_loc
5839 REAL(KIND=8),
DIMENSION(6) :: B_ext_l
5840 REAL(KIND=8) :: muhl, ray
5841 REAL(KIND=8),
DIMENSION(2) :: gaussp
5842 INTEGER :: mode, mesh_id1, k, i, index, l, m, ni
5843 INTEGER :: nb_procs, bloc_size, m_max_pad, code
5844 mpi_comm :: communicator
5846 DO i = 1,
SIZE(list_mode)
5850 mesh_id1 = mesh%i_d(m)
5851 j_loc = mesh%jj(:,m)
5853 b_in_loc(:,k) = b_in(j_loc,k,i)
5856 one_over_sigma_in_loc(:,k) = one_over_sigma_tot(j_loc,k,i)
5858 DO l = 1, mesh%gauss%l_G
5862 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
5865 b_ext_l(1) = sum(b_in_loc(:,1)*mesh%gauss%ww(:,l))
5866 b_ext_l(3) = sum(b_in_loc(:,3)*mesh%gauss%ww(:,l))
5867 b_ext_l(5) = sum(b_in_loc(:,5)*mesh%gauss%ww(:,l))
5868 b_ext_l(2) = sum(b_in_loc(:,2)*mesh%gauss%ww(:,l))
5869 b_ext_l(4) = sum(b_in_loc(:,4)*mesh%gauss%ww(:,l))
5870 b_ext_l(6) = sum(b_in_loc(:,6)*mesh%gauss%ww(:,l))
5872 DO ni = 1, mesh%gauss%n_w
5873 gaussp = gaussp + mesh%rr(:,mesh%jj(ni,m))*mesh%gauss%ww(ni,l)
5875 muhl=sum(mu_h_field(mesh%jj(:,m))*mesh%gauss%ww(:,l))
5878 j_exact_gauss(index,k,i)=jexact_gauss(k, gaussp, mode, mu_phi, sigma(m), &
5879 muhl, time, mesh_id1, b_ext_l)
5882 IF (jj_v_to_h(mesh%jj(1,m)) == -1)
THEN 5883 one_over_sigma_gauss(index,1,i) = 0.d0
5884 one_over_sigma_gauss(index,2,i) = 0.d0
5886 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m)
5889 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc(:,1)*mesh%gauss%ww(:,l))
5890 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc(:,2)*mesh%gauss%ww(:,l))
5897 CALL mpi_comm_size(communicator, nb_procs, code)
5898 bloc_size =
SIZE(j_exact_gauss,1)/nb_procs+1
5899 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5901 j_exact_gauss, one_over_sigma_gauss, j_over_sigma_gauss,&
5902 1, nb_procs, bloc_size, m_max_pad)
5907 mu_h_field, mu_phi, one_over_sigma_tot, time, sigma, j_over_sigma_gauss)
5913 #include "petsc/finclude/petsc.h" 5917 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
5918 INTEGER,
DIMENSION(:),
INTENT(IN) :: Dirichlet_bdy_H_sides
5919 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5920 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: B_in
5921 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_tot
5922 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
5923 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
5924 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
5925 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: J_over_sigma_gauss
5926 REAL(KIND=8),
DIMENSION(SIZE(Dirichlet_bdy_H_sides)*mesh%gauss%l_Gs,6,SIZE(list_mode)) :: J_exact_gauss
5927 REAL(KIND=8),
DIMENSION(SIZE(Dirichlet_bdy_H_sides)*mesh%gauss%l_Gs,2,SIZE(list_mode)) :: one_over_sigma_gauss
5928 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc
5929 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,6) :: B_in_loc
5930 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: one_over_sigma_in_loc
5931 REAL(KIND=8),
DIMENSION(6) :: B_ext_l
5932 REAL(KIND=8) :: muhl, ray
5933 REAL(KIND=8),
DIMENSION(2) :: gaussp
5934 INTEGER :: mode, mesh_id1, k, i, count, index, ls, ms, ni
5935 INTEGER :: nb_procs, bloc_size, m_max_pad, code, m1
5936 mpi_comm :: communicator
5938 DO i = 1,
SIZE(list_mode)
5941 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
5942 ms = dirichlet_bdy_h_sides(count)
5943 m1 = mesh%neighs(ms)
5944 mesh_id1 = mesh%i_d(m1)
5945 j_loc = mesh%jjs(:,ms)
5947 b_in_loc(:,k) = b_in(j_loc,k,i)
5950 one_over_sigma_in_loc(:,k) = one_over_sigma_tot(j_loc,k,i)
5953 DO ls = 1, mesh%gauss%l_Gs
5957 ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
5960 b_ext_l(1) = sum(b_in_loc(:,1)*mesh%gauss%wws(:,ls))
5961 b_ext_l(3) = sum(b_in_loc(:,3)*mesh%gauss%wws(:,ls))
5962 b_ext_l(5) = sum(b_in_loc(:,5)*mesh%gauss%wws(:,ls))
5963 b_ext_l(2) = sum(b_in_loc(:,2)*mesh%gauss%wws(:,ls))
5964 b_ext_l(4) = sum(b_in_loc(:,4)*mesh%gauss%wws(:,ls))
5965 b_ext_l(6) = sum(b_in_loc(:,6)*mesh%gauss%wws(:,ls))
5967 DO ni = 1, mesh%gauss%n_ws
5968 gaussp = gaussp + mesh%rr(:,mesh%jjs(ni,ms))*mesh%gauss%wws(ni,ls)
5970 muhl=sum(mu_h_field(mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
5973 j_exact_gauss(index,k,i)=jexact_gauss(k, gaussp, mode, mu_phi, sigma(m1),&
5974 muhl, time, mesh_id1, b_ext_l)
5977 IF (jj_v_to_h(mesh%jj(1,m1)) == -1)
THEN 5978 one_over_sigma_gauss(index,1,i) = 0.d0
5979 one_over_sigma_gauss(index,2,i) = 0.d0
5981 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m1)
5984 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc(:,1)*mesh%gauss%wws(:,ls))
5985 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc(:,2)*mesh%gauss%wws(:,ls))
5991 IF (
SIZE(dirichlet_bdy_h_sides).GE.1)
THEN 5992 CALL mpi_comm_size(communicator, nb_procs, code)
5993 bloc_size =
SIZE(j_exact_gauss,1)/nb_procs+1
5994 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5996 j_exact_gauss, one_over_sigma_gauss, j_over_sigma_gauss,&
5997 1, nb_procs, bloc_size, m_max_pad)
6003 mu_h_field, mu_phi, one_over_sigma_tot, time, sigma, j_over_sigma_gauss)
6009 #include "petsc/finclude/petsc.h" 6013 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_H
6015 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
6016 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: B_in
6017 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_tot
6018 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_H_field
6019 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
6020 REAL(KIND=8),
DIMENSION(mesh%me),
INTENT(IN) :: sigma
6021 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: J_over_sigma_gauss
6022 REAL(KIND=8),
DIMENSION(2*interface_H_mu%mes*mesh%gauss%l_Gs,6,SIZE(list_mode)) :: J_exact_gauss
6023 REAL(KIND=8),
DIMENSION(2*interface_H_mu%mes*mesh%gauss%l_Gs,2,SIZE(list_mode)) :: one_over_sigma_gauss
6024 REAL(KIND=8),
DIMENSION(6) :: B_ext_l
6025 REAL(KIND=8) :: muhl, diff, ref, ray
6026 REAL(KIND=8),
DIMENSION(2) :: gaussp
6027 INTEGER :: mode, k, i, mesh_id1, mesh_id2, ni
6028 INTEGER :: nb_procs, bloc_size, m_max_pad, code
6029 INTEGER :: ms, ms1, ms2, m1, m2, ls, index
6030 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc1, j_loc2
6031 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,6) :: B_in_loc1, B_in_loc2
6032 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: one_over_sigma_in_loc1,one_over_sigma_in_loc2
6033 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,mesh%gauss%l_Gs) :: w_cs
6034 mpi_comm :: communicator
6037 DO ms = 1, interface_h_mu%mes
6038 ms1 = interface_h_mu%mesh1(ms)
6039 ms2 = interface_h_mu%mesh2(ms)
6040 m1 = mesh%neighs(ms1)
6041 m2 = mesh%neighs(ms2)
6042 ref = 1.d-8+sum((mesh%rr(:,mesh%jjs(1,ms1)) - mesh%rr(:,mesh%jjs(2,ms1)))**2)
6043 diff = sum((mesh%rr(:,mesh%jjs(1,ms1)) - mesh%rr(:,mesh%jjs(1,ms2)))**2)
6044 IF (diff/ref .LT. 1.d-10)
THEN 6045 w_cs = mesh%gauss%wws
6047 DO ls = 1, mesh%gauss%l_Gs
6048 w_cs(1,ls)= mesh%gauss%wws(2,ls)
6049 w_cs(2,ls)= mesh%gauss%wws(1,ls)
6050 IF (mesh%gauss%n_ws==3) w_cs(mesh%gauss%n_ws,ls) = mesh%gauss%wws(mesh%gauss%n_ws,ls)
6051 WRITE(*,*)
' Ouaps! oder of shape functions changed?' 6056 DO i = 1,
SIZE(list_mode)
6059 DO ms = 1, interface_h_mu%mes
6060 ms2 = interface_h_mu%mesh2(ms)
6061 ms1 = interface_h_mu%mesh1(ms)
6062 m2 = mesh%neighs(ms2)
6063 m1 = mesh%neighs(ms1)
6064 mesh_id1 = mesh%i_d(m1)
6065 mesh_id2 = mesh%i_d(m2)
6066 j_loc1 = mesh%jjs(:,ms1)
6067 j_loc2 = mesh%jjs(:,ms2)
6069 b_in_loc1(:,k) = b_in(j_loc1,k,i)
6070 b_in_loc2(:,k) = b_in(j_loc2,k,i)
6073 one_over_sigma_in_loc1(:,k) = one_over_sigma_tot(j_loc1,k,i)
6074 one_over_sigma_in_loc2(:,k) = one_over_sigma_tot(j_loc2,k,i)
6077 DO ls = 1, mesh%gauss%l_Gs
6081 ray = sum(mesh%rr(1,mesh%jjs(:,ms1))*w_cs(:,ls))
6083 b_ext_l(1) = sum(b_in_loc1(:,1)*w_cs(:,ls))
6084 b_ext_l(3) = sum(b_in_loc1(:,3)*w_cs(:,ls))
6085 b_ext_l(5) = sum(b_in_loc1(:,5)*w_cs(:,ls))
6086 b_ext_l(2) = sum(b_in_loc1(:,2)*w_cs(:,ls))
6087 b_ext_l(4) = sum(b_in_loc1(:,4)*w_cs(:,ls))
6088 b_ext_l(6) = sum(b_in_loc1(:,6)*w_cs(:,ls))
6090 DO ni = 1, mesh%gauss%n_ws
6091 gaussp = gaussp + mesh%rr(:,mesh%jjs(ni,ms1))*w_cs(ni,ls)
6093 muhl=sum(mu_h_field(mesh%jjs(:,ms1))*w_cs(:,ls))
6096 j_exact_gauss(index,k,i)=jexact_gauss(k, gaussp, mode, mu_phi, sigma(m1),&
6097 muhl, time, mesh_id1, b_ext_l)
6100 IF (jj_v_to_h(mesh%jj(1,m1)) == -1)
THEN 6101 one_over_sigma_gauss(index,1,i) = 0.d0
6102 one_over_sigma_gauss(index,2,i) = 0.d0
6104 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m1)
6107 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc1(:,1)*w_cs(:,ls))
6108 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc1(:,2)*w_cs(:,ls))
6114 ray = sum(mesh%rr(1,mesh%jjs(:,ms2))*mesh%gauss%wws(:,ls))
6116 b_ext_l(1) = sum(b_in_loc2(:,1)*mesh%gauss%wws(:,ls))
6117 b_ext_l(3) = sum(b_in_loc2(:,3)*mesh%gauss%wws(:,ls))
6118 b_ext_l(5) = sum(b_in_loc2(:,5)*mesh%gauss%wws(:,ls))
6119 b_ext_l(2) = sum(b_in_loc2(:,2)*mesh%gauss%wws(:,ls))
6120 b_ext_l(4) = sum(b_in_loc2(:,4)*mesh%gauss%wws(:,ls))
6121 b_ext_l(6) = sum(b_in_loc2(:,6)*mesh%gauss%wws(:,ls))
6123 DO ni = 1, mesh%gauss%n_ws
6124 gaussp = gaussp + mesh%rr(:,mesh%jjs(ni,ms2))*mesh%gauss%wws(ni,ls)
6126 muhl=sum(mu_h_field(mesh%jjs(:,ms2))*mesh%gauss%wws(:,ls))
6129 j_exact_gauss(index,k,i)=jexact_gauss(k, gaussp, mode, mu_phi, sigma(m2),&
6130 muhl, time, mesh_id2, b_ext_l)
6133 IF (jj_v_to_h(mesh%jj(1,m2)) == -1)
THEN 6134 one_over_sigma_gauss(index,1,i) = 0.d0
6135 one_over_sigma_gauss(index,2,i) = 0.d0
6137 one_over_sigma_gauss(index,1,i) = 1.d0/sigma(m2)
6140 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc2(:,1)*mesh%gauss%wws(:,ls))
6141 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc2(:,2)*mesh%gauss%wws(:,ls))
6147 IF (interface_h_mu%mes.GE.1)
THEN 6148 CALL mpi_comm_size(communicator, nb_procs, code)
6149 bloc_size =
SIZE(j_exact_gauss,1)/nb_procs+1
6150 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
6152 j_exact_gauss, one_over_sigma_gauss, j_over_sigma_gauss,&
6153 1, nb_procs, bloc_size, m_max_pad)
6158 SUBROUTINE smb_sigma_neumann(communicator, mesh, Neumann_bdy_H_sides, list_mode, &
6159 one_over_sigma_tot, sigma_tot_gauss_neumann)
6165 #include "petsc/finclude/petsc.h" 6169 INTEGER,
DIMENSION(:),
INTENT(IN) :: Neumann_bdy_H_sides
6170 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
6171 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: one_over_sigma_tot
6172 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: sigma_tot_gauss_Neumann
6173 REAL(KIND=8),
DIMENSION(SIZE(Neumann_bdy_H_sides)*mesh%gauss%l_Gs,2,SIZE(list_mode)):: one_over_sigma_gauss
6174 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc
6175 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: one_over_sigma_in_loc
6176 INTEGER :: mode, k, i, count, index, ls, ms
6177 INTEGER :: nb_procs, bloc_size, m_max_pad, code
6178 mpi_comm :: communicator
6180 DO i = 1,
SIZE(list_mode)
6183 DO count = 1,
SIZE(neumann_bdy_h_sides)
6184 ms = neumann_bdy_h_sides(count)
6185 j_loc = mesh%jjs(:,ms)
6187 one_over_sigma_in_loc(:,k) = one_over_sigma_tot(j_loc,k,i)
6190 DO ls = 1, mesh%gauss%l_Gs
6192 one_over_sigma_gauss(index,1,i) = sum(one_over_sigma_in_loc(:,1)*mesh%gauss%wws(:,ls))
6193 one_over_sigma_gauss(index,2,i) = sum(one_over_sigma_in_loc(:,2)*mesh%gauss%wws(:,ls))
6198 sigma_tot_gauss_neumann=one_over_sigma_gauss
6200 IF (
SIZE(neumann_bdy_h_sides).GE.1)
THEN 6201 CALL mpi_comm_size(communicator, nb_procs, code)
6202 bloc_size =
SIZE(one_over_sigma_gauss,1)/nb_procs+1
6203 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
6214 vv = x/ max(x*x, 1.d-20)
6218
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 solver(my_ksp, b, x, reinit, verbose)
subroutine dirichlet_cavities(communicator, interface_H_phi, mesh, js_D)
integer, dimension(:), allocatable neumann_bdy_pmag_sides
subroutine, public extract(xghost, ks, ke, LA, phi)
subroutine, public create_my_ghost(mesh, LA, ifrom)
subroutine, public fft_par_var_eta_prod_t_dcl(communicator, eta, H_mesh, c_in, c_out, nb_procs, bloc_size, m_max_pad, time, temps)
real(kind=8) function, dimension(nb_angles, ne-nb+1) minus_one_over_mu(H_mesh, angles, nb_angles, nb, ne, time)
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)
real(kind=8) function one_over_x(x)
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)
integer, dimension(:), allocatable neumann_bdy_h_sides
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)
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 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 create_local_petsc_matrix(communicator, LA, matrix, clean)
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 mat_dirichlet_maxwell(H_mesh, jj_v_to_H, Dirichlet_bdy_H_sides, mode, stab, mu_H_field, LA_H, H_p_phi_mat1, H_p_phi_mat2, sigma_np, sigma)
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)
subroutine, public maxwell_decouple_with_b(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)
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 init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
subroutine dirichlet_rhs(js_D, bs_D, b)
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_in)
subroutine dirichlet_nodes_parallel(mesh, list_dirichlet_sides, js_d)
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
real(kind=8) function, dimension(nb_angles, ne-nb+1) one_over_mu(H_mesh, angles, nb_angles, nb, ne, time)
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, H_pert, sigma_curl_gauss, J_over_sigma_gauss)
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 dirichlet_m_parallel(matrix, glob_js_D)
integer, dimension(:), allocatable neumann_bdy_phi_sides
type(my_verbose), public talk_to_me
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
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 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)
integer, dimension(:), allocatable dirichlet_bdy_h_sides
real(kind=8), dimension(:,:), pointer rjs
real(kind=8), dimension(:,:), pointer wws
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, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)
subroutine smb_sigma_neumann(communicator, mesh, Neumann_bdy_H_sides, list_mode, one_over_sigma_tot, sigma_tot_gauss_Neumann)
real(kind=8), dimension(:,:,:,:), pointer dw_s
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)