10 REAL(KIND=8),
PARAMETER,
PRIVATE :: alpha=0.6d0
18 interface_h_mu, hn, bn, phin, hn1, bn1, phin1, vel, stab_in, sigma_in, &
19 r_fourier, index_fourier, mu_h_field, mu_phi, time, dt_in, rem, list_mode, &
20 h_phi_per, la_h, la_pmag, la_phi, la_mhd, sigma_ns_in, jj_v_to_h)
37 TYPE(mesh_type),
INTENT(IN) :: h_mesh, phi_mesh, pmag_mesh
38 TYPE(interface_type),
INTENT(IN) :: interface_h_phi, interface_h_mu
39 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
40 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: vel
41 REAL(KIND=8),
DIMENSION(H_mesh%np,6,SIZE(list_mode)),
INTENT(INOUT) :: hn, hn1
42 REAL(KIND=8),
DIMENSION(H_mesh%np,6,SIZE(list_mode)),
INTENT(INOUT) :: bn, bn1
43 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: phin, phin1
44 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab_in
45 REAL(KIND=8),
INTENT(IN) :: r_fourier
46 INTEGER,
INTENT(IN) :: index_fourier
47 REAL(KIND=8),
INTENT(IN) :: mu_phi, time, dt_in, rem
48 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_in, mu_h_field
52 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: sigma_ns_in
53 INTEGER,
DIMENSION(:) ,
INTENT(IN) :: jj_v_to_h
56 REAL(KIND=8),
SAVE :: dt
57 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE,
SAVE :: sigma_ns_bar
59 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE,
SAVE :: sigma_np
60 REAL(KIND=8),
DIMENSION(:),
ALLOCATABLE,
SAVE :: sigma
61 REAL(KIND=8),
SAVE :: sigma_min
62 TYPE(dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: h_mode_global_js_d
63 TYPE(dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: h_global_d
65 TYPE(dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: pmag_mode_global_js_d
66 TYPE(dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: pmag_global_d
68 TYPE(dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: phi_mode_global_js_d
69 TYPE(dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: phi_global_d
70 INTEGER,
DIMENSION(:),
POINTER,
SAVE :: pmag_js_d, phi_js_d
71 INTEGER,
DIMENSION(:),
ALLOCATABLE,
SAVE :: dirichlet_bdy_h_sides
72 LOGICAL,
SAVE :: once=.true.
73 INTEGER,
SAVE :: m_max_c
75 REAL(KIND=8),
DIMENSION(3),
SAVE :: stab
76 INTEGER,
SAVE :: my_petscworld_rank
77 REAL(KIND=8),
ALLOCATABLE ,
DIMENSION(:,:,:),
SAVE :: sigma_curl_bdy
79 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_w,H_mesh%me) :: sigma_nj_m
80 REAL(KIND=8),
DIMENSION(H_mesh%gauss%l_G*H_mesh%me,6,SIZE(list_mode)) :: sigma_curl
81 REAL(KIND=8),
DIMENSION(H_mesh%np,6,SIZE(list_mode)) :: j_over_sigma
82 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),6,SIZE(Hn,3)) :: h_ns
83 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),2,SIZE(Hn,3)) :: sigma_tot
84 LOGICAL,
ALLOCATABLE,
DIMENSION(:) :: dir_pmag
85 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: rhs_h
86 REAL(KIND=8),
DIMENSION(phi_mesh%np,2) :: rhs_phi
89 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),6,SIZE(Hn,3)) :: h_pert
91 REAL(KIND=8),
DIMENSION(SIZE(Hn,1),6,SIZE(Hn,3)) :: nl, b_ext
92 REAL(KIND=8),
DIMENSION(3) :: temps_par
93 INTEGER,
POINTER,
DIMENSION(:) :: h_ifrom, pmag_ifrom, phi_ifrom, h_p_phi_ifrom
94 REAL(KIND=8),
DIMENSION(phi_mesh%np, 2) :: phin_p1
95 REAL(KIND=8),
DIMENSION(H_mesh%np, 6) :: hn_p1
96 INTEGER :: mode, k, i, n, m, ms, code, nj, j
97 INTEGER :: nb_procs, bloc_size, m_max_pad
98 REAL(KIND=8) :: tps, nr_vel, tps_tot, tps_cumul, norm
100 REAL(KIND=8) :: one_and_half
101 DATA one_and_half/1.5d0/
104 #include "petsc/finclude/petscsys.h"
105 #include "petsc/finclude/petscmat.h"
106 #include "petsc/finclude/petscksp.h"
107 #include "petsc/finclude/petscvec.h"
108 #include "petsc/finclude/petscvec.h90"
109 petscerrorcode :: ierr
110 mpi_comm,
DIMENSION(:),
POINTER :: comm_one_d
111 mat,
DIMENSION(:),
POINTER,
SAVE :: h_p_phi_mat1, h_p_phi_mat2
112 mat :: tampon1, tampon2, precond1, precond2
113 ksp,
DIMENSION(:),
POINTER,
SAVE :: h_p_phi_ksp1, h_p_phi_ksp2
114 vec,
SAVE :: vx_1, vb_1, vx_1_ghost, vx_2, vb_2, vx_2_ghost
127 CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
134 IF (inputs%if_quasi_static_approx)
THEN
140 n =
SIZE(h_ifrom)+
SIZE(pmag_ifrom)+
SIZE(phi_ifrom)
141 ALLOCATE(h_p_phi_ifrom(n))
142 IF (
SIZE(h_ifrom)/=0)
THEN
143 h_p_phi_ifrom(1:
SIZE(h_ifrom)) = h_ifrom
145 IF (
SIZE(pmag_ifrom)/=0)
THEN
146 h_p_phi_ifrom(
SIZE(h_ifrom)+1:
SIZE(h_ifrom)+
SIZE(pmag_ifrom)) = pmag_ifrom
148 IF (
SIZE(phi_ifrom)/=0)
THEN
149 h_p_phi_ifrom(
SIZE(h_ifrom)+
SIZE(pmag_ifrom)+1:)=phi_ifrom
152 n = 3*h_mesh%dom_np + pmag_mesh%dom_np + phi_mesh%dom_np
153 CALL veccreateghost(comm_one_d(1), n, &
154 petsc_determine,
SIZE(h_p_phi_ifrom), h_p_phi_ifrom, vx_1, ierr)
155 CALL vecghostgetlocalform(vx_1, vx_1_ghost, ierr)
156 CALL vecduplicate(vx_1, vb_1, ierr)
157 CALL veccreateghost(comm_one_d(1), n, &
158 petsc_determine,
SIZE(h_p_phi_ifrom), h_p_phi_ifrom, vx_2, ierr)
159 CALL vecghostgetlocalform(vx_2, vx_2_ghost, ierr)
160 CALL vecduplicate(vx_2, vb_2, ierr)
164 ALLOCATE(sigma(
SIZE(sigma_in)))
165 sigma = sigma_in * rem
168 CALL mpi_allreduce(minval(sigma),sigma_min,1,mpi_double_precision, mpi_min,comm_one_d(1), ierr)
174 IF (inputs%type_pb==
'mhd')
THEN
177 stab = stab_in*(1/sigma_min+1.d0)
181 IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid)
THEN
182 stab = stab_in*(1/(minval(inputs%sigma_fluid)*rem)+1.d0)
186 nr_vel =
norm_sf(comm_one_d,
'L2', h_mesh, list_mode, vel)
188 IF (nr_vel .LE. 1.d-10)
THEN
191 stab = stab_in*(1/sigma_min)
197 stab = stab_in*(1/sigma_min+1.d0)
207 m_max_c =
SIZE(list_mode)
211 IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid)
THEN
214 ALLOCATE(sigma_ns_bar(
SIZE(hn,1)))
215 DO n = 1,
SIZE(h_mesh%rr,2)
227 DO nj = 1, h_mesh%gauss%n_w
229 IF (jj_v_to_h(j) == -1)
THEN
230 sigma_nj_m(nj,m) = sigma(m)
234 sigma_nj_m(nj,m) = sigma_ns_bar(j)
241 sigma_nj_m(:,m) = sigma(m)
245 ALLOCATE(sigma_np(
SIZE(hn,1)))
248 DO nj = 1, h_mesh%gauss%n_w
249 sigma_np(h_mesh%jj(nj,m)) = sigma_nj_m(nj,m)
258 ALLOCATE (dir_pmag(maxval(pmag_mesh%sides)))
260 DO ms = 1,
SIZE(dir_pmag)
261 IF (minval(abs(inputs%list_dirichlet_sides_H-ms)) == 0)
THEN
262 dir_pmag(ms) = .true.
264 IF (minval(abs(inputs%list_inter_H_phi-ms)) == 0)
THEN
265 dir_pmag(ms) = .true.
269 CALL
dirichlet_nodes(pmag_mesh%jjs, pmag_mesh%sides, dir_pmag, pmag_js_d)
274 ALLOCATE(pmag_global_d(m_max_c))
276 ALLOCATE(pmag_global_d(i)%DRL(
SIZE(pmag_mode_global_js_d(i)%DIL)))
277 pmag_global_d(i)%DRL = 0.d0
284 DO ms = 1, h_mesh%mes
285 IF (minval(abs(h_mesh%sides(ms)-inputs%list_dirichlet_sides_H))/=0) cycle
286 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))) .LT.1d-12) cycle
289 ALLOCATE(dirichlet_bdy_h_sides(n))
291 DO ms = 1, h_mesh%mes
292 IF (minval(abs(h_mesh%sides(ms)-inputs%list_dirichlet_sides_H))/=0) cycle
293 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))) .LT.1d-12) cycle
295 dirichlet_bdy_h_sides(n) = ms
299 ALLOCATE(h_global_d(m_max_c))
301 ALLOCATE(h_global_d(i)%DRL(
SIZE(h_mode_global_js_d(i)%DIL)))
310 ALLOCATE(phi_global_d(m_max_c))
312 ALLOCATE(phi_global_d(i)%DRL(
SIZE(phi_mode_global_js_d(i)%DIL)))
313 phi_global_d(i)%DRL = 0.d0
318 ALLOCATE(h_p_phi_mat1(m_max_c),h_p_phi_ksp1(m_max_c))
319 ALLOCATE(h_p_phi_mat2(m_max_c),h_p_phi_ksp2(m_max_c))
320 IF (
SIZE(dirichlet_bdy_h_sides).GE.1)
THEN
321 ALLOCATE(sigma_curl_bdy(h_mesh%gauss%l_Gs*
SIZE(dirichlet_bdy_h_sides),6,
SIZE(list_mode)))
324 ALLOCATE(sigma_curl_bdy(0,0,0))
325 sigma_curl_bdy = 0.d0
348 mode,mu_h_field, mu_phi, one_and_half/dt, stab, r_fourier, index_fourier, &
349 la_h, la_pmag, la_phi, h_p_phi_mat1(i), h_p_phi_mat2(i), sigma_np)
356 mu_h_field, sigma, la_h, h_p_phi_mat1(i), h_p_phi_mat2(i))
362 mode, stab, mu_h_field, la_h, h_p_phi_mat1(i), h_p_phi_mat2(i), sigma_np)
364 IF (inputs%my_periodic%nb_periodic_pairs/=0)
THEN
366 h_phi_per%perlist, h_p_phi_mat1(i), la_mhd)
368 h_phi_per%perlist, h_p_phi_mat2(i), la_mhd)
385 CALL
init_solver(inputs%my_par_H_p_phi,h_p_phi_ksp1(i),h_p_phi_mat1(i),comm_one_d(1),&
386 solver=inputs%my_par_H_p_phi%solver,precond=inputs%my_par_H_p_phi%precond)
387 CALL
init_solver(inputs%my_par_H_p_phi,h_p_phi_ksp2(i),h_p_phi_mat2(i),comm_one_d(1),&
388 solver=inputs%my_par_H_p_phi%solver,precond=inputs%my_par_H_p_phi%precond)
392 CALL matdestroy(h_p_phi_mat1(i),ierr)
393 CALL matdestroy(h_p_phi_mat2(i),ierr)
401 CALL mpi_comm_rank(petsc_comm_world, my_petscworld_rank, code)
404 IF(inputs%if_arpack)
THEN
405 IF (h_mesh%me/=0)
THEN
406 IF (inputs%if_permeability_variable_in_theta)
THEN
407 CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
408 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
409 bloc_size =
SIZE(b_ext,1)/nb_procs+1
411 h_mesh, hn, bn, nb_procs, bloc_size, m_max_pad, time,temps_par)
413 h_mesh, hn1, bn1, nb_procs, bloc_size, m_max_pad, time,temps_par)
417 bn(:,k,i) = mu_h_field*hn(:,k,i)
418 bn1(:,k,i) = mu_h_field*hn1(:,k,i)
428 nr_vel =
norm_sf(comm_one_d,
'L2', h_mesh, list_mode, vel)
430 IF(inputs%if_quasi_static_approx)
THEN
439 IF (nr_vel .LE. 1.d-10)
THEN
442 CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
443 bloc_size =
SIZE(vel,1)/nb_procs+1
444 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
451 IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid)
THEN
455 DO nj = 1, h_mesh%gauss%n_w
458 IF (jj_v_to_h(j) /= -1)
THEN
459 h_ns(j,:,:) = 2*hn(j,:,:)- hn1(j,:,:)
460 sigma_tot(j,:,:) = sigma_ns_in(jj_v_to_h(j),:,:) * rem
462 DO i = 1,
SIZE(list_mode)
465 sigma_tot(j,1,i) = sigma(m)
473 sigma_ns_bar, sigma_tot, sigma_curl)
474 IF (
SIZE(dirichlet_bdy_h_sides).GE.1)
THEN
476 sigma_ns_bar, sigma_tot, sigma_curl_bdy)
478 sigma_curl_bdy = 0.d0
482 mu_h_field, mu_phi, sigma_tot, time, j_over_sigma)
485 sigma_curl_bdy = 0.d0
489 IF (inputs%if_permeability_variable_in_theta)
THEN
491 CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
492 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
493 bloc_size =
SIZE(b_ext,1)/nb_procs+1
496 h_mesh, b_ext, h_pert, nb_procs, bloc_size, m_max_pad, time,temps_par)
499 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
511 rhs_h(:,k) = (4*bn(:,k,i)-bn1(:,k,i))/(2*dt)
514 rhs_phi(:,k) = mu_phi*(4*phin(:,k,i)-phin1(:,k,i))/(2*dt)
519 rhs_h, nl(:,:,i), la_h, la_phi, vb_1, vb_2, b_ext(:,:,i), h_pert(:,:,i), &
520 sigma_curl(:,:,i), j_over_sigma(:,:,i))
527 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
531 CALL
courant_mu(h_mesh, interface_h_mu,sigma, mu_h_field, time,mode, nl(:,:,i), &
532 la_h, vb_1, vb_2, b_ext(:,:,i))
533 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
539 sigma, mu_h_field, time, mode, nl(:,:,i), stab, la_h, vb_1,vb_2,b_ext(:,:,i), &
540 j_over_sigma(:,:,i), sigma_curl_bdy(:,:,i))
545 CALL
surf_int(h_mesh,phi_mesh,interface_h_phi,interface_h_mu,inputs%list_dirichlet_sides_H, &
546 sigma,mu_phi,mu_h_field, time,mode,la_h, la_phi,vb_1,vb_2, r_fourier, index_fourier)
547 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
553 IF (inputs%my_periodic%nb_periodic_pairs/=0)
THEN
554 CALL
periodic_rhs_petsc(h_phi_per%n_bord, h_phi_per%list, h_phi_per%perlist, vb_1, la_mhd)
555 CALL
periodic_rhs_petsc(h_phi_per%n_bord, h_phi_per%list, h_phi_per%perlist, vb_2, la_mhd)
562 pmag_global_d(i)%DRL = 0.d0
563 CALL
dirichlet_rhs(pmag_mode_global_js_d(i)%DIL-1,pmag_global_d(i)%DRL,vb_1)
564 CALL
dirichlet_rhs(pmag_mode_global_js_d(i)%DIL-1,pmag_global_d(i)%DRL,vb_2)
566 IF (
SIZE(phi_js_d)>0)
THEN
571 phi_global_d(i)%DRL(1:n) =
phiexact(1,phi_mesh%rr(1:2,phi_js_d), mode, mu_phi, time)
572 IF (
SIZE(phi_global_d(i)%DRL)>n) phi_global_d(i)%DRL(n+1:)=0.d0
573 CALL
dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_1)
574 phi_global_d(i)%DRL(1:n) =
phiexact(2,phi_mesh%rr(1:2,phi_js_d), mode, mu_phi, time)
575 IF (
SIZE(phi_global_d(i)%DRL)>n) phi_global_d(i)%DRL(n+1:)=0.d0
576 CALL
dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_2)
580 phi_global_d(i)%DRL=0.d0
581 CALL
dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_1)
582 CALL
dirichlet_rhs(la_phi%loc_to_glob(1,phi_js_d)-1, phi_global_d(i)%DRL, vb_2)
586 h_global_d(i)%DRL = 0.d0
587 CALL
dirichlet_rhs(h_mode_global_js_d(i)%DIL-1,h_global_d(i)%DRL,vb_1)
588 CALL
dirichlet_rhs(h_mode_global_js_d(i)%DIL-1,h_global_d(i)%DRL,vb_2)
590 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
595 IF (inputs%my_par_H_p_phi%verbose .AND. (i==1))
WRITE(*,*)
'start solving'
598 CALL
solver(h_p_phi_ksp1(i),vb_1,vx_1,reinit=.false.,
verbose=inputs%my_par_H_p_phi%verbose)
600 CALL vecghostupdatebegin(vx_1,insert_values,scatter_forward,ierr)
601 CALL vecghostupdateend(vx_1,insert_values,scatter_forward,ierr)
602 IF (h_mesh%me/=0)
THEN
603 CALL
extract(vx_1_ghost,1,1,la_mhd,hn_p1(:,1))
604 CALL
extract(vx_1_ghost,2,2,la_mhd,hn_p1(:,4))
605 CALL
extract(vx_1_ghost,3,3,la_mhd,hn_p1(:,5))
607 IF (phi_mesh%me/=0)
THEN
608 CALL
extract(vx_1_ghost,5,5,la_mhd,phin_p1(:,1))
611 CALL
solver(h_p_phi_ksp2(i),vb_2,vx_2,reinit=.false.,
verbose=inputs%my_par_H_p_phi%verbose)
612 CALL vecghostupdatebegin(vx_2,insert_values,scatter_forward,ierr)
613 CALL vecghostupdateend(vx_2,insert_values,scatter_forward,ierr)
614 IF (h_mesh%me/=0)
THEN
615 CALL
extract(vx_2_ghost,1,1,la_mhd,hn_p1(:,2))
616 CALL
extract(vx_2_ghost,2,2,la_mhd,hn_p1(:,3))
617 CALL
extract(vx_2_ghost,3,3,la_mhd,hn_p1(:,6))
619 IF (phi_mesh%me/=0)
THEN
620 CALL
extract(vx_2_ghost,5,5,la_mhd,phin_p1(:,2))
623 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
631 IF (h_mesh%me /=0)
THEN
636 IF (phi_mesh%me /=0 )
THEN
649 IF (h_mesh%me /=0)
THEN
650 hn1(:,:,i) = hn(:,:,i)
660 bn1(:,:,i) = bn(:,:,i)
662 bn(:,1,i) = hn_p1(:,1)
663 bn(:,4,i) = hn_p1(:,4)
664 bn(:,5,i) = hn_p1(:,5)
666 bn(:,2,i) = hn_p1(:,2)
667 bn(:,3,i) = hn_p1(:,3)
668 bn(:,6,i) = hn_p1(:,6)
672 IF (phi_mesh%me /= 0)
THEN
673 phin1(:,:,i) = phin(:,:,i)
675 phin(:,1,i) = phin_p1(:,1)
677 phin(:,2,i) = phin_p1(:,2)
679 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
685 IF (h_mesh%me /=0)
THEN
686 IF (inputs%if_permeability_variable_in_theta)
THEN
687 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
688 bloc_size =
SIZE(bn,1)/nb_procs+1
690 h_mesh, bn, hn, nb_procs, bloc_size, m_max_pad, time,temps_par)
694 hn(:,k,i) = bn(:,k,i)/mu_h_field
701 IF (inputs%verbose_divergence)
THEN
702 norm =
norm_sf(comm_one_d,
'L2', h_mesh, list_mode, bn)
703 talk_to_me%div_B_L2 =
norm_sf(comm_one_d,
'div', h_mesh, list_mode, bn)/norm
714 mode, mu_h_field, mu_phi, c_mass, stab, r_fourier, index_fourier, &
715 la_h, la_pmag, la_phi, h_p_phi_mat1, h_p_phi_mat2, sigma_np)
726 INTEGER,
INTENT(IN) :: mode
727 REAL(KIND=8),
INTENT(IN) :: mu_phi, c_mass
728 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
729 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field
730 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_np
731 REAL(KIND=8),
OPTIONAL :: r_fourier
732 INTEGER,
OPTIONAL :: index_fourier
733 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_ws,phi_mesh%gauss%l_Gs) :: w_cs
734 REAL(KIND=8),
DIMENSION(2, H_mesh%gauss%n_w, phi_mesh%gauss%l_Gs, H_mesh%mes) :: dw_cs
735 INTEGER :: m, l, ms, ls, ni, nj, k, i, j, &
736 n_ws1, n_ws2, n_w2, n_w1, m1, m2, ki, kj,ib,jb, ms1, ms2
737 REAL(KIND=8) :: x, y, hm1, stab_div, stab_colle_h_phi
738 REAL(KIND=8) :: ray, error
739 LOGICAL :: mark=.false.
740 REAL(KIND=8),
DIMENSION(3,H_mesh%gauss%n_w,pmag_mesh%gauss%n_w) :: thpmag
741 REAL(KIND=8),
DIMENSION(pmag_mesh%gauss%n_w,pmag_mesh%gauss%n_w) :: tpmag
742 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: th
743 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: htob
744 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_w,phi_mesh%gauss%n_w):: tphi
772 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: hsij
773 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_w,phi_mesh%gauss%n_w) :: phisij
774 REAL(KIND=8),
DIMENSION(6,phi_mesh%gauss%n_w,phi_mesh%gauss%n_w) :: sij
775 REAL(KIND=8),
DIMENSION(6,phi_mesh%gauss%n_w,phi_mesh%gauss%n_w) :: smuij
796 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: ww_h
798 REAL(KIND=8),
DIMENSION(:,:,:,:),
POINTER :: dw_h
800 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: rj_h
802 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: ww_phi
804 REAL(KIND=8),
DIMENSION(:,:,:,:),
POINTER :: dw_phi
805 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%n_w,H_mesh%gauss%l_G) :: dwp
806 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_w,H_mesh%gauss%l_G) :: wwp
808 REAL(KIND=8),
DIMENSION(:,:),
POINTER :: rj_phi
809 REAL(KIND=8),
DIMENSION(2,phi_mesh%gauss%l_Gs) :: gauss1, gauss2
811 REAL(KIND=8) ::
ref, diff, mu_h, c_mu_phi, muhl, &
812 dzmuhl, drmuhl, c_div, hloc, viscolm, xij, eps
814 REAL(KIND=8) :: c_sym=.0d0
817 REAL(KIND=8) :: c_lap
822 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
823 INTEGER ,
DIMENSION(3*H_mesh%gauss%n_w+pmag_mesh%gauss%n_w+ &
phi_mesh%gauss%n_w) :: idxn, jdxn
826 INTEGER :: n_wpmag, n_wh, n_wphi, ix, jx
828 REAL(KIND=8),
DIMENSION(2) :: dmu_field
829 REAL(KIND=8),
DIMENSION(2) :: gauss_pt
830 INTEGER,
DIMENSION(1) :: gauss_pt_id
831 REAL(KIND=8),
DIMENSION(1) :: dummy_mu_bar
834 REAL(KIND=8) :: sigma_np_gauss
836 #include "petsc/finclude/petsc.h"
837 mat :: h_p_phi_mat1, h_p_phi_mat2
838 petscerrorcode :: ierr
839 CALL matzeroentries(h_p_phi_mat1, ierr)
840 CALL matzeroentries(h_p_phi_mat2, ierr)
841 CALL matsetoption(h_p_phi_mat1, mat_row_oriented, petsc_false, ierr)
842 CALL matsetoption(h_p_phi_mat2, mat_row_oriented, petsc_false, ierr)
846 stab_colle_h_phi = stab(2)
850 c_mu_phi = c_mass*mu_phi
852 ww_h => h_mesh%gauss%ww
853 dw_h => h_mesh%gauss%dw
854 rj_h => h_mesh%gauss%rj
855 ww_phi => phi_mesh%gauss%ww
856 dw_phi => phi_mesh%gauss%dw
857 rj_phi => phi_mesh%gauss%rj
859 n_wh = h_mesh%gauss%n_w
860 n_wpmag = pmag_mesh%gauss%n_w
861 n_wphi = phi_mesh%gauss%n_w
868 mesh_id = h_mesh%i_d(m)
873 DO l = 1, h_mesh%gauss%l_G
874 hloc = sqrt(sum(h_mesh%gauss%rj(:,m)))**(2*alpha)
878 IF(inputs%if_use_fem_integration_for_mu_bar)
THEN
879 muhl = sum(mu_h_field(h_mesh%jj(:,m))*ww_h(:,l))
880 drmuhl = sum(mu_h_field(h_mesh%jj(:,m))*dw_h(1,:,l,m))
881 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m))*dw_h(2,:,l,m))
883 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jj(:,m))*ww_h(:,l))
884 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jj(:,m))*ww_h(:,l))
893 sigma_np_gauss = sum(sigma_np(h_mesh%jj(:,m))*ww_h(:,l))
897 if ( muhl < 0.5d0)
then
898 write(*,*)
'!! Warning: mu_bar is almost negative in the matrix construction, value = ', muhl
902 c_div = stab_div*hloc
906 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
907 ray = ray + h_mesh%rr(1,i)*ww_h(ni,l)
910 DO ni = 1, h_mesh%gauss%n_w
911 DO nj = 1, h_mesh%gauss%n_w
914 th(1,ni,nj) = th(1,ni,nj) + rj_h(l,m) * ray* ( &
915 c_mass*ww_h(ni,l)*ww_h(nj,l) &
916 + (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) &
919 + c_div*(muhl*(ww_h(ni,l)/ray+dw_h(1,ni,l,m)) + ww_h(ni,l)*0) &
920 *(ww_h(nj,l)/ray+dw_h(1,nj,l,m))) &
921 + rj_h(l,m) * ray* ( &
922 (dw_h(2,ni,l,m)*dzmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
925 th(2,ni,nj) = th(2,ni,nj)+ rj_h(l,m) * ray* ( &
926 mode/ray**2 * ww_h(ni,l)*(ww_h(nj,l)+ray*dw_h(1,nj,l,m))/(sigma_np_gauss*muhl) &
929 + 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)) &
930 + rj_h(l,m) * ray* ( &
931 ((mode/ray)*ww_h(ni,l)*drmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
933 htob(2,ni,nj) = htob(2,ni,nj)+ rj_h(l,m) * ray* ( &
934 - mode/ray**2 * ww_h(ni,l)*(ww_h(nj,l)+ray*dw_h(1,nj,l,m))/(sigma_np_gauss*muhl) &
937 - 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)) &
938 - rj_h(l,m) * ray* ( &
939 ((mode/ray)*ww_h(ni,l)*drmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
941 th(3,ni,nj) = th(3,ni,nj)+ rj_h(l,m) * ray* ( &
942 - dw_h(2,ni,l,m)*dw_h(1,nj,l,m)/(sigma_np_gauss*muhl) &
945 + 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))) &
946 + rj_h(l,m) * ray* ( &
947 - dw_h(2,ni,l,m)*drmuhl*ww_h(nj,l)*(-1/(sigma_np_gauss*muhl**2)))
949 th(4,ni,nj) = th(4,ni,nj)+ rj_h(l,m) * ray* ( &
950 mode/ray**2 * ww_h(nj,l)*(ww_h(ni,l)+ray*dw_h(1,ni,l,m))/(sigma_np_gauss*muhl) &
953 + c_div*(mode/ray)*ww_h(ni,l)*muhl*(ww_h(nj,l)/ray+dw_h(1,nj,l,m)))
955 htob(4,ni,nj) = htob(4,ni,nj)+ rj_h(l,m) * ray* ( &
956 - mode/ray**2 * ww_h(nj,l)*(ww_h(ni,l)+ray*dw_h(1,ni,l,m))/(sigma_np_gauss*muhl) &
959 - c_div*(mode/ray)*ww_h(ni,l)*muhl*(ww_h(nj,l)/ray+dw_h(1,nj,l,m)))
961 th(5,ni,nj) = th(5,ni,nj) + rj_h(l,m) * ray* ( &
962 c_mass*ww_h(ni,l)*ww_h(nj,l) &
963 + (dw_h(2,ni,l,m)*dw_h(2,nj,l,m) &
964 + 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) &
967 +c_div*muhl*mode**2/ray**2*ww_h(ni,l)*ww_h(nj,l)) &
968 + rj_h(l,m)*ray*((dw_h(2,ni,l,m)*dzmuhl*ww_h(nj,l) &
969 + (ww_h(ni,l)/ray+dw_h(1,ni,l,m))*drmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
971 th(6,ni,nj) = th(6,ni,nj) + rj_h(l,m) * ray* (&
972 + mode/ray*dw_h(2,ni,l,m)*ww_h(nj,l)/(sigma_np_gauss*muhl) &
975 +c_div*mode/ray*muhl*ww_h(ni,l)*(dw_h(2,nj,l,m)))
977 htob(6,ni,nj) = htob(6,ni,nj) + rj_h(l,m) * ray* (&
978 - mode/ray*dw_h(2,ni,l,m)*ww_h(nj,l)/(sigma_np_gauss*muhl) &
981 - c_div*mode/ray*muhl*ww_h(ni,l)*(dw_h(2,nj,l,m)))
983 th(7,ni,nj) = th(7,ni,nj)+ rj_h(l,m) * ray* ( &
984 - dw_h(1,ni,l,m)*dw_h(2,nj,l,m)/(sigma_np_gauss*muhl) &
987 + 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))) &
988 + rj_h(l,m)*ray*((-dw_h(1,ni,l,m)*dzmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
990 th(8,ni,nj) = th(8,ni,nj) + rj_h(l,m) * ray* (&
991 + (mode/ray)*ww_h(ni,l)*dw_h(2,nj,l,m)/(sigma_np_gauss*muhl) &
994 + c_div*mode/ray*muhl*ww_h(nj,l)*(dw_h(2,ni,l,m))) &
995 + rj_h(l,m)*ray*(((mode/ray)*ww_h(ni,l)*dzmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
999 htob(8,ni,nj) = htob(8,ni,nj) + rj_h(l,m) * ray* (&
1000 - mode/ray*dw_h(2,nj,l,m)*ww_h(ni,l)/(sigma_np_gauss*muhl) &
1003 - c_div*mode/ray*muhl*ww_h(nj,l)*(dw_h(2,ni,l,m))) &
1004 - rj_h(l,m)*ray*(((mode/ray)*ww_h(ni,l)*dzmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
1008 th(9,ni,nj) = th(9,ni,nj) + rj_h(l,m) * ray* ( &
1009 c_mass*ww_h(ni,l)*ww_h(nj,l) &
1010 + (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) &
1013 + c_div*(muhl*dw_h(2,ni,l,m) + ww_h(ni,l)*0)*(dw_h(2,nj,l,m))) &
1014 + rj_h(l,m)*ray*((dw_h(1,ni,l,m)*drmuhl*ww_h(nj,l))*(-1/(sigma_np_gauss*muhl**2)))
1122 i = h_mesh%jj(ni, m)
1123 ib = la_h%loc_to_glob(ki,i)
1128 j = h_mesh%jj(nj, m)
1129 jb = la_h%loc_to_glob(kj,j)
1133 IF ((ki == 1) .AND. (kj == 1))
THEN
1134 mat_loc1(ix,jx) = th(1,ni,nj)
1135 mat_loc2(ix,jx) = th(1,ni,nj)
1136 ELSEIF ((ki == 1) .AND. (kj == 2))
THEN
1137 mat_loc1(ix,jx) = th(2,ni,nj)
1138 mat_loc2(ix,jx) = htob(2,ni,nj)
1139 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN
1140 mat_loc1(ix,jx) = th(3,ni,nj)
1141 mat_loc2(ix,jx) = th(3,ni,nj)
1142 ELSEIF ((ki == 2) .AND. (kj == 1))
THEN
1143 mat_loc1(ix,jx) = th(4,ni,nj)
1144 mat_loc2(ix,jx) = htob(4,ni,nj)
1145 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
1146 mat_loc1(ix,jx) = th(5,ni,nj)
1147 mat_loc2(ix,jx) = th(5,ni,nj)
1148 ELSEIF ((ki == 2) .AND. (kj == 3))
THEN
1149 mat_loc1(ix,jx) = th(6,ni,nj)
1150 mat_loc2(ix,jx) = htob(6,ni,nj)
1151 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN
1152 mat_loc1(ix,jx) = th(7,ni,nj)
1153 mat_loc2(ix,jx) = th(7,ni,nj)
1154 ELSEIF ((ki == 3) .AND. (kj == 2))
THEN
1155 mat_loc1(ix,jx) = th(8,ni,nj)
1156 mat_loc2(ix,jx) = htob(8,ni,nj)
1157 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN
1158 mat_loc1(ix,jx) = th(9,ni,nj)
1159 mat_loc2(ix,jx) = th(9,ni,nj)
1166 CALL matsetvalues(h_p_phi_mat1, 3*n_wh, idxn(1:3*n_wh), 3*n_wh, jdxn(1:3*n_wh), &
1167 mat_loc1(1:3*n_wh,1:3*n_wh), add_values, ierr)
1168 CALL matsetvalues(h_p_phi_mat2, 3*n_wh, idxn(1:3*n_wh), 3*n_wh, jdxn(1:3*n_wh), &
1169 mat_loc2(1:3*n_wh,1:3*n_wh), add_values, ierr)
1173 DO m = 1, pmag_mesh%me
1174 hloc = stab_div*sqrt(sum(pmag_mesh%gauss%rj(:,m)))**(2*(1-alpha))
1176 DO l = 1, pmag_mesh%gauss%l_G
1184 DO ni = 1, pmag_mesh%gauss%n_w
1185 i = pmag_mesh%jj(ni,m)
1186 ray = ray + pmag_mesh%rr(1,i)*pmag_mesh%gauss%ww(ni,l)
1188 viscolm = hloc*muhl*pmag_mesh%gauss%rj(l,m)
1189 DO nj = 1, pmag_mesh%gauss%n_w
1190 j = pmag_mesh%jj(nj, m)
1191 DO ni = 1, pmag_mesh%gauss%n_w
1192 i = pmag_mesh%jj(ni, m)
1196 xij = xij + pmag_mesh%gauss%dw(k,nj,l,m) * pmag_mesh%gauss%dw(k,ni,l,m)
1199 tpmag(ni,nj) = tpmag(ni,nj) + ray * viscolm* xij &
1200 + viscolm*mode**2*pmag_mesh%gauss%ww(ni,l)*pmag_mesh%gauss%ww(nj,l)/ray
1205 DO ni = 1, pmag_mesh%gauss%n_w
1206 i = pmag_mesh%jj(ni, m)
1207 ib = la_pmag%loc_to_glob(1,i)
1209 DO nj = 1, pmag_mesh%gauss%n_w
1210 j = pmag_mesh%jj(nj, m)
1211 jb = la_pmag%loc_to_glob(1,j)
1215 CALL matsetvalues(h_p_phi_mat1, n_wpmag, idxn(1:n_wpmag), n_wpmag, jdxn(1:n_wpmag), &
1216 tpmag(1:n_wpmag,1:n_wpmag), add_values, ierr)
1217 CALL matsetvalues(h_p_phi_mat2, n_wpmag, idxn(1:n_wpmag), n_wpmag, jdxn(1:n_wpmag), &
1218 tpmag(1:n_wpmag,1:n_wpmag), add_values, ierr)
1223 DO m = 1, pmag_mesh%me
1224 mesh_id = h_mesh%i_d(m)
1225 IF (h_mesh%gauss%n_w==3)
THEN
1226 dwp=h_mesh%gauss%dw(:,:,:,m)
1229 dwp(:,1,:) = h_mesh%gauss%dw(:,1,:,m) + 0.5d0*(h_mesh%gauss%dw(:,5,:,m)+h_mesh%gauss%dw(:,6,:,m))
1230 dwp(:,2,:) = h_mesh%gauss%dw(:,2,:,m) + 0.5d0*(h_mesh%gauss%dw(:,6,:,m)+h_mesh%gauss%dw(:,4,:,m))
1231 dwp(:,3,:) = h_mesh%gauss%dw(:,3,:,m) + 0.5d0*(h_mesh%gauss%dw(:,4,:,m)+h_mesh%gauss%dw(:,5,:,m))
1232 wwp(1,:) = h_mesh%gauss%ww(1,:) + 0.5d0*(h_mesh%gauss%ww(5,:)+h_mesh%gauss%ww(6,:))
1233 wwp(2,:) = h_mesh%gauss%ww(2,:) + 0.5d0*(h_mesh%gauss%ww(6,:)+h_mesh%gauss%ww(4,:))
1234 wwp(3,:) = h_mesh%gauss%ww(3,:) + 0.5d0*(h_mesh%gauss%ww(4,:)+h_mesh%gauss%ww(5,:))
1238 DO l = 1, h_mesh%gauss%l_G
1240 DO ni = 1, h_mesh%gauss%n_w
1242 ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
1246 IF(inputs%if_use_fem_integration_for_mu_bar)
THEN
1247 muhl = stab_div*ray*h_mesh%gauss%rj(l,m)*sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
1250 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jj(:,m))*ww_h(:,l))
1251 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jj(:,m))*ww_h(:,l))
1252 gauss_pt_id(1)=mesh_id
1254 muhl=dummy_mu_bar(1)
1255 muhl = stab_div*ray*h_mesh%gauss%rj(l,m)*muhl
1257 DO nj = 1, pmag_mesh%gauss%n_w
1258 j = pmag_mesh%jj(nj, m)
1259 DO ni = 1, h_mesh%gauss%n_w
1260 i = h_mesh%jj(ni, m)
1261 thpmag(1,ni,nj) = thpmag(1,ni,nj) + muhl*dwp(1,nj,l)*h_mesh%gauss%ww(ni,l)
1262 thpmag(2,ni,nj) = thpmag(2,ni,nj) - muhl*mode*wwp(nj,l)*h_mesh%gauss%ww(ni,l)/ray
1263 thpmag(3,ni,nj) = thpmag(3,ni,nj) + muhl*dwp(2,nj,l)*h_mesh%gauss%ww(ni,l)
1273 i = h_mesh%jj(ni, m)
1280 ib = la_h%loc_to_glob(k,i)
1281 ix = (k-1)*n_wh + ni
1284 j = pmag_mesh%jj(nj, m)
1285 jb = la_pmag%loc_to_glob(1,j)
1288 mat_loc1(ix,jx) = thpmag(k,ni,nj)
1289 mat_loc2(ix,jx) = eps*thpmag(k,ni,nj)
1294 CALL matsetvalues(h_p_phi_mat1, 3*n_wh, idxn(1:3*n_wh), n_wpmag, jdxn(1:n_wpmag), &
1295 mat_loc1(1:3*n_wh,1:n_wpmag), add_values, ierr)
1296 CALL matsetvalues(h_p_phi_mat2, 3*n_wh, idxn(1:3*n_wh), n_wpmag, jdxn(1:n_wpmag), &
1297 mat_loc2(1:3*n_wh,1:n_wpmag), add_values, ierr)
1301 DO l = 1, h_mesh%gauss%l_G
1303 DO ni = 1, h_mesh%gauss%n_w
1305 ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
1307 x = stab_div*ray*h_mesh%gauss%rj(l,m)
1308 DO nj = 1, h_mesh%gauss%n_w
1309 j = h_mesh%jj(nj, m)
1310 DO ni = 1, pmag_mesh%gauss%n_w
1311 i = pmag_mesh%jj(ni, m)
1312 thpmag(1,nj,ni) = thpmag(1,nj,ni) - x*dwp(1,ni,l)*h_mesh%gauss%ww(nj,l)
1313 thpmag(2,nj,ni) = thpmag(2,nj,ni) + x*mode*wwp(ni,l)*h_mesh%gauss%ww(nj,l)/ray
1314 thpmag(3,nj,ni) = thpmag(3,nj,ni) - x*dwp(2,ni,l)*h_mesh%gauss%ww(nj,l)
1322 i = pmag_mesh%jj(ni, m)
1323 ib = la_pmag%loc_to_glob(1,i)
1333 j = h_mesh%jj(nj, m)
1334 jb = la_h%loc_to_glob(k,j)
1335 jx = (k-1)*n_wh + nj
1337 mat_loc1(ix,jx) = thpmag(k,nj,ni)
1338 mat_loc2(ix,jx) = eps*thpmag(k,nj,ni)
1342 CALL matsetvalues(h_p_phi_mat1, n_wpmag, idxn(1:n_wpmag), 3*n_wh, jdxn(1:3*n_wh), &
1343 mat_loc1(1:n_wpmag,1:3*n_wh), add_values, ierr)
1344 CALL matsetvalues(h_p_phi_mat2, n_wpmag, idxn(1:n_wpmag), 3*n_wh, jdxn(1:3*n_wh), &
1345 mat_loc2(1:n_wpmag,1:3*n_wh), add_values, ierr)
1352 DO m = 1,phi_mesh%me
1356 DO l = 1, phi_mesh%gauss%l_G
1360 DO ni = 1, phi_mesh%gauss%n_w; i = phi_mesh%jj(ni,m)
1361 ray = ray + phi_mesh%rr(1,i)*ww_phi(ni,l)
1364 DO ni = 1, phi_mesh%gauss%n_w
1365 DO nj = 1, phi_mesh%gauss%n_w
1373 tphi(ni,nj) = tphi(ni,nj) + rj_phi(l,m) * ray* (c_mass+c_lap)*mu_phi &
1374 *(dw_phi(1,ni,l,m)*dw_phi(1,nj,l,m)+dw_phi(2,ni,l,m)*dw_phi(2,nj,l,m) &
1375 +mode**2/ray**2*ww_phi(ni,l)*ww_phi(nj,l))
1382 DO ni = 1, phi_mesh%gauss%n_w
1383 i = phi_mesh%jj(ni, m)
1384 ib = la_phi%loc_to_glob(1,i)
1386 DO nj = 1, phi_mesh%gauss%n_w
1387 j = phi_mesh%jj(nj, m)
1388 jb = la_phi%loc_to_glob(1,j)
1392 CALL matsetvalues(h_p_phi_mat1, n_wphi, idxn(1:n_wphi), n_wphi, jdxn(1:n_wphi), &
1393 tphi(1:n_wphi,1:n_wphi), add_values, ierr)
1394 CALL matsetvalues(h_p_phi_mat2, n_wphi, idxn(1:n_wphi), n_wphi, jdxn(1:n_wphi), &
1395 tphi(1:n_wphi,1:n_wphi), add_values, ierr)
1403 CALL
gauss(phi_mesh)
1404 n_ws1 = h_mesh%gauss%n_ws
1405 n_ws2 = phi_mesh%gauss%n_ws
1406 n_w1 = h_mesh%gauss%n_w
1407 n_w2 = phi_mesh%gauss%n_w
1409 IF (h_mesh%gauss%n_ws == n_ws)
THEN
1411 DO ms = 1, interface_h_phi%mes
1413 ms2 = interface_h_phi%mesh2(ms)
1414 m2 = phi_mesh%neighs(ms2)
1415 ms1 = interface_h_phi%mesh1(ms)
1416 m1 = h_mesh%neighs(ms1)
1418 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
1419 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - phi_mesh%rr(:,phi_mesh%jjs(1,ms2)))**2)
1420 IF (diff/
ref .LT. 1.d-10)
THEN
1424 w_cs(1,ls)= wws(2,ls)
1425 w_cs(2,ls)= wws(1,ls)
1426 w_cs(3,ls)= wws(3,ls)
1427 WRITE(*,*)
' Ouaps! oder of shape functions changed?'
1432 gauss2(1,ls) = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))*phi_mesh%gauss%wws(:,ls))
1433 gauss2(2,ls) = sum(phi_mesh%rr(2,phi_mesh%jjs(:,ms2))*phi_mesh%gauss%wws(:,ls))
1434 gauss1(1,ls) = sum( h_mesh%rr(1, h_mesh%jjs(:,ms1))* h_mesh%gauss%wws(:,ls))
1435 gauss1(2,ls) = sum( h_mesh%rr(2, h_mesh%jjs(:,ms1))* h_mesh%gauss%wws(:,ls))
1439 ref = sqrt(1.d-8+sum(gauss2(:,ls2)**2))
1442 diff = sqrt(sum((gauss2(:,ls2)-gauss1(:,ls1))**2))
1443 IF (diff .LT. 1.d-10)
THEN
1444 dw_cs(:,:,ls2,ms1) = h_mesh%gauss%dw_s(:,:,ls1,ms1)
1449 IF (.NOT.mark)
WRITE(*,*)
' BUG '
1455 DO ms = 1, interface_h_phi%mes
1457 ms2 = interface_h_phi%mesh2(ms)
1458 m2 = phi_mesh%neighs(ms2)
1459 ms1 = interface_h_phi%mesh1(ms)
1460 m1 = h_mesh%neighs(ms1)
1462 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
1463 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - phi_mesh%rr(:,phi_mesh%jjs(1,ms2)))**2)
1464 IF (diff/
ref .LT. 1.d-10)
THEN
1466 w_cs(1,ls)= wws(1,ls)+0.5*wws(3,ls)
1467 w_cs(2,ls)= wws(2,ls)+0.5*wws(3,ls)
1472 w_cs(1,ls)= wws(2,ls)+0.5*wws(3,ls)
1473 w_cs(2,ls)= wws(1,ls)+0.5*wws(3,ls)
1475 WRITE(*,*)
' Ouaps! oder of shape functions changed?'
1480 dw_cs(1,:,ls,ms1) = h_mesh%gauss%dw(1,:,1,m1)
1481 dw_cs(2,:,ls,ms1) = h_mesh%gauss%dw(2,:,1,m1)
1488 DO ms = 1, interface_h_phi%mes
1490 ms2 = interface_h_phi%mesh2(ms)
1491 ms1 = interface_h_phi%mesh1(ms)
1492 m2 = phi_mesh%neighs(ms2)
1493 m1 = h_mesh%neighs(ms1)
1494 mu_h = sum(mu_h_field(h_mesh%jj(:,m1)))/h_mesh%gauss%n_w
1496 hm1 = stab_colle_h_phi/sum(rjs(:,ms2))
1510 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1511 x = hm1*rjs(ls,ms2)*ray
1513 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1519 y = x * w_cs(ni,ls)*w_cs(nj,ls)/muhl
1521 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(rnorms(2,ls,ms2)**2)
1522 hsij(4,ni,nj) = hsij(4,ni,nj) - y*rnorms(1,ls,ms2)*rnorms(2,ls,ms2)
1523 hsij(5,ni,nj) = hsij(5,ni,nj) + y
1524 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(rnorms(1,ls,ms2)**2)
1535 i = interface_h_phi%jjs1(ni,ms)
1536 ib = la_h%loc_to_glob(ki,i)
1537 ix = (ki-1)*n_ws1+ni
1541 j = interface_h_phi%jjs1(nj,ms)
1542 jb = la_h%loc_to_glob(kj,j)
1543 jx = (kj-1)*n_ws1+nj
1545 IF ((ki == 1) .AND. (kj == 1))
THEN
1546 mat_loc1(ix,jx) = hsij(1,ni,nj)
1547 mat_loc2(ix,jx) = hsij(1,ni,nj)
1548 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN
1549 mat_loc1(ix,jx) = hsij(4,ni,nj)
1550 mat_loc2(ix,jx) = hsij(4,ni,nj)
1551 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN
1552 mat_loc1(ix,jx) = hsij(4,nj,ni)
1553 mat_loc2(ix,jx) = hsij(4,nj,ni)
1554 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
1555 mat_loc1(ix,jx) = hsij(5,ni,nj)
1556 mat_loc2(ix,jx) = hsij(5,ni,nj)
1557 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN
1558 mat_loc1(ix,jx) = hsij(6,ni,nj)
1559 mat_loc2(ix,jx) = hsij(6,ni,nj)
1565 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1566 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1567 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1568 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1574 DO ls = 1, phi_mesh%gauss%l_Gs
1577 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1580 x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1583 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1584 drmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(1,:,ls,ms1))
1585 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(2,:,ls,ms1))
1591 y = x*w_cs(ni,ls)*w_cs(nj,ls)
1598 hsij(2,ni,nj) = hsij(2,ni,nj) + y * (-mode/ray)*(-rnorms(1,ls,ms2))/muhl
1599 hsij(3,ni,nj) = hsij(3,ni,nj) + y * mode/ray *(-rnorms(1,ls,ms2))/muhl
1600 hsij(5,ni,nj) = hsij(5,ni,nj) + y * (-1/ray) *(-rnorms(1,ls,ms2))/muhl
1601 hsij(8,ni,nj) = hsij(8,ni,nj) + y * (-mode/ray)*(-rnorms(2,ls,ms2))/muhl
1602 hsij(9,ni,nj) = hsij(9,ni,nj) + y * mode/ray *(-rnorms(2,ls,ms2))/muhl
1604 hsij(1,ni,nj) = hsij(1,ni,nj) - y*(-rnorms(2,ls,ms2))*(-(dzmuhl/muhl**2))
1605 hsij(4,ni,nj) = hsij(4,ni,nj) - y*(-rnorms(2,ls,ms2))*( (drmuhl/muhl**2))
1606 hsij(5,ni,nj) = hsij(5,ni,nj) &
1607 + y*(-rnorms(1,ls,ms2)*drmuhl-rnorms(2,ls,ms2)*dzmuhl)/muhl**2
1608 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(-rnorms(1,ls,ms2))*( (drmuhl/muhl**2))
1609 hsij(7,ni,nj) = hsij(7,ni,nj) + y*(-rnorms(1,ls,ms2))*(-(dzmuhl/muhl**2))
1624 i = interface_h_phi%jjs1(ni,ms)
1625 ib = la_h%loc_to_glob(ki,i)
1626 ix = (ki-1)*n_ws1 + ni
1630 j = interface_h_phi%jjs1(nj,ms)
1631 jb = la_h%loc_to_glob(kj,j)
1632 jx = (kj-1)*n_ws1 + nj
1635 IF ( (ki == 1) .AND. (kj == 1))
THEN
1636 mat_loc1(ix,jx) = hsij(1,ni,nj)
1637 mat_loc2(ix,jx) = hsij(1,ni,nj)
1638 ELSE IF ( (ki == 1) .AND. (kj == 3))
THEN
1639 mat_loc1(ix,jx) = hsij(4,ni,nj)
1640 mat_loc2(ix,jx) = hsij(4,ni,nj)
1642 ELSE IF ( (ki == 2) .AND. (kj == 1))
THEN
1643 mat_loc1(ix,jx) = hsij(2,ni,nj)
1644 mat_loc2(ix,jx) = hsij(3,ni,nj)
1645 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
1646 mat_loc1(ix,jx) = hsij(5,ni,nj)
1647 mat_loc2(ix,jx) = hsij(5,ni,nj)
1648 ELSEIF ( (ki == 2) .AND. (kj == 3))
THEN
1649 mat_loc1(ix,jx) = hsij(8,ni,nj)
1650 mat_loc2(ix,jx) = hsij(9,ni,nj)
1652 ELSE IF ( (ki == 3) .AND. (kj == 1))
THEN
1653 mat_loc1(ix,jx) = hsij(7,ni,nj)
1654 mat_loc2(ix,jx) = hsij(7,ni,nj)
1655 ELSE IF ( (ki == 3) .AND. (kj == 3))
THEN
1656 mat_loc1(ix,jx) = hsij(6,ni,nj)
1657 mat_loc2(ix,jx) = hsij(6,ni,nj)
1664 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1665 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1666 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1667 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1676 i = interface_h_phi%jjs1(ni,ms)
1677 ib = la_h%loc_to_glob(ki,i)
1678 ix = (ki-1)*n_ws1 + ni
1682 j = interface_h_phi%jjs1(nj,ms)
1683 jb = la_h%loc_to_glob(kj,j)
1684 jx = (kj-1)*n_ws1 + nj
1686 IF ( (kj == 2) .AND. (ki == 1))
THEN
1687 mat_loc1(ix,jx) = hsij(2,nj,ni)
1688 mat_loc2(ix,jx) = hsij(3,nj,ni)
1689 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN
1690 mat_loc1(ix,jx) = hsij(5,nj,ni)
1691 mat_loc2(ix,jx) = hsij(5,nj,ni)
1692 ELSEIF ( (kj == 2) .AND. (ki == 3))
THEN
1693 mat_loc1(ix,jx) = hsij(8,nj,ni)
1694 mat_loc2(ix,jx) = hsij(9,nj,ni)
1701 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1702 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1703 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
1704 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
1707 DO ls = 1, phi_mesh%gauss%l_Gs
1710 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1713 x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1716 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1729 hsij(1,ni,nj) = hsij(1,ni,nj) - y*(-rnorms(2,ls,ms2))*(dw_cs(2,nj,ls,ms1)/muhl)
1730 hsij(4,ni,nj) = hsij(4,ni,nj) - y*(-rnorms(2,ls,ms2))*(-dw_cs(1,nj,ls,ms1)/muhl)
1731 hsij(5,ni,nj) = hsij(5,ni,nj) &
1732 + y*(-rnorms(2,ls,ms2))*(-dw_cs(2,nj,ls,ms1)/muhl) &
1733 - y*(-rnorms(1,ls,ms2))*( dw_cs(1,nj,ls,ms1)/muhl)
1734 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(-rnorms(1,ls,ms2))*(-dw_cs(1,nj,ls,ms1)/muhl)
1735 hsij(7,ni,nj) = hsij(7,ni,nj) + y*(-rnorms(1,ls,ms2))*(dw_cs(2,nj,ls,ms1)/muhl)
1749 i = interface_h_phi%jjs1(ni,ms)
1750 ib = la_h%loc_to_glob(ki,i)
1751 ix = (ki-1)*n_ws1 + ni
1755 j = h_mesh%jj(nj,m1)
1756 jb = la_h%loc_to_glob(kj,j)
1757 jx = (kj-1)*n_w1 + nj
1759 IF ((ki == 1) .AND. (kj == 1))
THEN
1760 mat_loc1(ix,jx) = hsij(1,ni,nj)
1761 mat_loc2(ix,jx) = hsij(1,ni,nj)
1762 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN
1763 mat_loc1(ix,jx) = hsij(4,ni,nj)
1764 mat_loc2(ix,jx) = hsij(4,ni,nj)
1765 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
1766 mat_loc1(ix,jx) = hsij(5,ni,nj)
1767 mat_loc2(ix,jx) = hsij(5,ni,nj)
1768 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN
1769 mat_loc1(ix,jx) = hsij(6,ni,nj)
1770 mat_loc2(ix,jx) = hsij(6,ni,nj)
1771 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN
1772 mat_loc1(ix,jx) = hsij(7,ni,nj)
1773 mat_loc2(ix,jx) = hsij(7,ni,nj)
1780 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
1781 mat_loc1(1:3*n_ws1,1:3*n_w1), add_values, ierr)
1782 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
1783 mat_loc2(1:3*n_ws1,1:3*n_w1), add_values, ierr)
1791 i = h_mesh%jj(ni,m1)
1792 ib = la_h%loc_to_glob(ki,i)
1793 ix = (ki-1)*n_w1 + ni
1797 j = interface_h_phi%jjs1(nj,ms)
1798 jb = la_h%loc_to_glob(kj,j)
1799 jx = (kj-1)*n_ws1 + nj
1801 IF ((kj == 1) .AND. (ki == 1))
THEN
1802 mat_loc1(ix,jx) = hsij(1,nj,ni)
1803 mat_loc2(ix,jx) = hsij(1,nj,ni)
1804 ELSEIF ((kj == 1) .AND. (ki == 3))
THEN
1805 mat_loc1(ix,jx) = hsij(4,nj,ni)
1806 mat_loc2(ix,jx) = hsij(4,nj,ni)
1807 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN
1808 mat_loc1(ix,jx) = hsij(5,nj,ni)
1809 mat_loc2(ix,jx) = hsij(5,nj,ni)
1810 ELSEIF ((kj == 3) .AND. (ki == 3))
THEN
1811 mat_loc1(ix,jx) = hsij(6,nj,ni)
1812 mat_loc2(ix,jx) = hsij(6,nj,ni)
1813 ELSEIF ((kj == 3) .AND. (ki == 1))
THEN
1814 mat_loc1(ix,jx) = hsij(7,nj,ni)
1815 mat_loc2(ix,jx) = hsij(7,nj,ni)
1821 CALL matsetvalues(h_p_phi_mat1, 3*n_w1, idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
1822 mat_loc1(1:3*n_w1,1:3*n_ws1), add_values, ierr)
1823 CALL matsetvalues(h_p_phi_mat2, 3*n_w1, idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
1824 mat_loc2(1:3*n_w1,1:3*n_ws1), add_values, ierr)
1837 DO ls = 1, phi_mesh%gauss%l_Gs
1840 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1841 x = hm1*rjs(ls,ms2)*ray
1846 phisij(ni,nj) = phisij(ni,nj) + x*mode**2/ray**2*wws(ni,ls)*wws(nj,ls)
1857 i = interface_h_phi%jjs2(ni,ms)
1858 ib = la_phi%loc_to_glob(1,i)
1861 j = interface_h_phi%jjs2(nj,ms)
1862 jb = la_phi%loc_to_glob(1,j)
1866 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
1867 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
1868 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
1869 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
1875 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1876 x = hm1*rjs(ls,ms2)*ray
1881 phisij(ni,nj) = phisij(ni,nj) + x*( &
1882 (dw_s(2,ni,ls,ms2)*rnorms(1,ls,ms2) - dw_s(1,ni,ls,ms2)*rnorms(2,ls,ms2))* &
1883 (dw_s(2,nj,ls,ms2)*rnorms(1,ls,ms2) - dw_s(1,nj,ls,ms2)*rnorms(2,ls,ms2)))
1893 i = phi_mesh%jj(ni, m2)
1894 ib = la_phi%loc_to_glob(1,i)
1897 j = phi_mesh%jj(nj, m2)
1898 jb = la_phi%loc_to_glob(1,j)
1902 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), n_w2, jdxn(1:n_w2), &
1903 phisij(1:n_w2,1:n_w2), add_values, ierr)
1904 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), n_w2, jdxn(1:n_w2), &
1905 phisij(1:n_w2,1:n_w2), add_values, ierr)
1919 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1920 x = hm1*rjs(ls,ms2)*ray
1922 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
1927 sij(3,ni,nj) = sij(3,ni,nj) + x*(mode/ray)*w_cs(ni,ls)*wws(nj,ls)
1928 smuij(3,ni,nj) = smuij(3,ni,nj) + x*(mode/ray)*w_cs(ni,ls)*wws(nj,ls)/muhl
1932 sij(4,:,:) = -sij(3,:,:)
1934 smuij(4,:,:) = -smuij(3,:,:)
1944 i = interface_h_phi%jjs1(ni,ms)
1945 ib = la_h%loc_to_glob(ki,i)
1948 j = interface_h_phi%jjs2(nj,ms)
1949 jb = la_phi%loc_to_glob(1,j)
1953 CALL matsetvalues(h_p_phi_mat1, n_ws1, idxn(1:n_ws1), n_ws2, jdxn(1:n_ws2), &
1954 sij(3,1:n_ws1,1:n_ws2), add_values, ierr)
1955 CALL matsetvalues(h_p_phi_mat2, n_ws1, idxn(1:n_ws1), n_ws2, jdxn(1:n_ws2), &
1956 sij(4,1:n_ws1,1:n_ws2), add_values, ierr)
1965 i = interface_h_phi%jjs2(ni,ms)
1966 ib = la_phi%loc_to_glob(1,i)
1969 j = interface_h_phi%jjs1(nj,ms)
1970 jb = la_h%loc_to_glob(kj,j)
1975 mat_loc1(ni,nj) = smuij(3,nj,ni)
1976 mat_loc2(ni,nj) = smuij(4,nj,ni)
1981 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws1, jdxn(1:n_ws1), &
1982 mat_loc1(1:n_ws2,1:n_ws1), add_values, ierr)
1983 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws1, jdxn(1:n_ws1), &
1984 mat_loc2(1:n_ws2,1:n_ws1), add_values, ierr)
1994 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
1995 x = hm1*rjs(ls,ms2)*ray
1998 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2004 sij(1,ni,nj) = sij(1,ni,nj) + &
2005 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))
2007 smuij(1,ni,nj) = smuij(1,ni,nj) + &
2008 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
2010 sij(5,ni,nj) = sij(5,ni,nj) + &
2011 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))
2013 smuij(5,ni,nj) = smuij(5,ni,nj) + &
2014 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
2028 i = interface_h_phi%jjs1(ni,ms)
2029 ib = la_h%loc_to_glob(ki,i)
2030 ix = (ki-1)*n_ws1 + ni
2033 j = phi_mesh%jj(nj,m2)
2034 jb = la_phi%loc_to_glob(1,j)
2038 mat_loc1(ix,jx) = sij(1,ni,nj)
2039 mat_loc2(ix,jx) = sij(1,ni,nj)
2040 ELSEIF (ki == 3)
THEN
2041 mat_loc1(ix,jx) = sij(5,ni,nj)
2042 mat_loc2(ix,jx) = sij(5,ni,nj)
2047 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), n_w2, jdxn(1:n_w2), &
2048 mat_loc1(1:3*n_ws1,1:n_w2), add_values, ierr)
2049 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), n_w2, jdxn(1:n_w2), &
2050 mat_loc2(1:3*n_ws1,1:n_w2), add_values, ierr)
2058 i = phi_mesh%jj(ni,m2)
2059 ib = la_phi%loc_to_glob(1,i)
2064 j = interface_h_phi%jjs1(nj,ms)
2065 jb = la_h%loc_to_glob(kj,j)
2066 jx = (kj-1)*n_ws1 + nj
2072 mat_loc1(ix,jx) = smuij(1,nj,ni)
2073 mat_loc2(ix,jx) = smuij(1,nj,ni)
2075 ELSEIF (kj == 3)
THEN
2079 mat_loc1(ix,jx) = smuij(5,nj,ni)
2080 mat_loc2(ix,jx) = smuij(5,nj,ni)
2087 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
2088 mat_loc1(1:n_w2,1:3*n_ws1), add_values, ierr)
2089 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
2090 mat_loc2(1:n_w2,1:3*n_ws1), add_values, ierr)
2105 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2108 x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2111 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2112 drmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(1,:,ls,ms1))
2113 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(2,:,ls,ms1))
2120 y = x * wws(ni,ls)*w_cs(nj,ls)
2126 sij(1,ni,nj) = sij(1,ni,nj) + y*( mode/ray)**2*rnorms(1,ls,ms2)/muhl
2127 sij(3,ni,nj) = sij(3,ni,nj) + y*( mode/ray**2)*rnorms(1,ls,ms2)/muhl
2128 sij(4,ni,nj) = sij(4,ni,nj) + y*(-mode/ray**2)*rnorms(1,ls,ms2)/muhl
2129 sij(5,ni,nj) = sij(5,ni,nj) + y*( mode/ray)**2*rnorms(2,ls,ms2)/muhl
2131 sij(3,ni,nj) = sij(3,ni,nj) &
2132 + y*(mode/ray)*(drmuhl*rnorms(1,ls,ms2)+dzmuhl*rnorms(2,ls,ms2))*(-1/muhl**2)
2133 sij(4,ni,nj) = sij(4,ni,nj) &
2134 - y*(mode/ray)*(drmuhl*rnorms(1,ls,ms2)+dzmuhl*rnorms(2,ls,ms2))*(-1/muhl**2)
2148 i = interface_h_phi%jjs2(ni,ms)
2149 ib = la_phi%loc_to_glob(1,i)
2154 j = interface_h_phi%jjs1(nj,ms)
2155 jb = la_h%loc_to_glob(kj,j)
2156 jx = (kj-1)*n_ws1 + nj
2159 mat_loc1(ix,jx) = sij(1,ni,nj)
2160 mat_loc2(ix,jx) = sij(1,ni,nj)
2161 ELSEIF (kj == 2)
THEN
2162 mat_loc1(ix,jx) = sij(3,ni,nj)
2163 mat_loc2(ix,jx) = sij(4,ni,nj)
2164 ELSEIF (kj == 3)
THEN
2165 mat_loc1(ix,jx) = sij(5,ni,nj)
2166 mat_loc2(ix,jx) = sij(5,ni,nj)
2171 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), 3*n_ws1, jdxn(1:3*n_ws1), &
2172 mat_loc1(1:n_ws2,1:3*n_ws1), add_values, ierr)
2173 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), 3*n_ws1, jdxn(1:3*n_ws1), &
2174 mat_loc2(1:n_ws2,1:3*n_ws1), add_values, ierr)
2182 i = interface_h_phi%jjs1(ni,ms)
2183 ib = la_h%loc_to_glob(ki,i)
2184 ix = (ki-1)*n_ws1 + ni
2187 j = interface_h_phi%jjs2(nj,ms)
2188 jb = la_phi%loc_to_glob(1,j)
2192 mat_loc1(ix,jx) = sij(1,nj,ni)
2193 mat_loc2(ix,jx) = sij(1,nj,ni)
2194 ELSEIF (ki == 2)
THEN
2195 mat_loc1(ix,jx) = sij(3,nj,ni)
2196 mat_loc2(ix,jx) = sij(4,nj,ni)
2197 ELSEIF (ki == 3)
THEN
2198 mat_loc1(ix,jx) = sij(5,nj,ni)
2199 mat_loc2(ix,jx) = sij(5,nj,ni)
2204 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1, idxn(1:3*n_ws1), n_ws2, jdxn(1:n_ws2), &
2205 mat_loc1(1:3*n_ws1,1:n_ws2), add_values, ierr)
2206 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1, idxn(1:3*n_ws1), n_ws2, jdxn(1:n_ws2), &
2207 mat_loc2(1:3*n_ws1,1:n_ws2), add_values, ierr)
2214 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2217 x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2220 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2224 y = x*wws(ni,ls)*mode/ray
2226 sij(3,ni,nj) = sij(3,ni,nj) + &
2229 y*(dw_cs(2,nj,ls,ms1)*rnorms(2,ls,ms2)/muhl + dw_cs(1,nj,ls,ms1)*rnorms(1,ls,ms2)/muhl)
2234 sij(4,:,:) = -sij(3,:,:)
2240 i = interface_h_phi%jjs2(ni,ms)
2241 ib = la_phi%loc_to_glob(1,i)
2244 j = h_mesh%jj(nj,m1)
2245 jb = la_h%loc_to_glob(kj,j)
2249 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_w1, jdxn(1:n_w1), &
2250 sij(3,1:n_ws2,1:n_w1), add_values, ierr)
2251 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_w1, jdxn(1:n_w1), &
2252 sij(4,1:n_ws2,1:n_w1), add_values, ierr)
2261 i = h_mesh%jj(ni,m1)
2262 ib = la_h%loc_to_glob(ki,i)
2265 j = interface_h_phi%jjs2(nj,ms)
2266 jb = la_phi%loc_to_glob(1,j)
2268 mat_loc1(ix,jx) = sij(3,nj,ni)
2269 mat_loc2(ix,jx) = sij(4,nj,ni)
2272 CALL matsetvalues(h_p_phi_mat1, n_w1, idxn(1:n_w1), n_ws2, jdxn(1:n_ws2), &
2273 mat_loc1(1:n_w1,1:n_ws2), add_values, ierr)
2274 CALL matsetvalues(h_p_phi_mat2, n_w1, idxn(1:n_w1), n_ws2, jdxn(1:n_ws2), &
2275 mat_loc2(1:n_w1,1:n_ws2), add_values, ierr)
2282 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2285 x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2288 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2292 y = x*(dw_s(2,ni,ls,ms2)*rnorms(1,ls,ms2) - dw_s(1,ni,ls,ms2)*rnorms(2,ls,ms2))
2297 sij(1,ni,nj) = sij(1,ni,nj) + y * (dw_cs(2,nj,ls,ms1)/muhl)
2298 sij(5,ni,nj) = sij(5,ni,nj) + (-y)*(dw_cs(1,nj,ls,ms1)/muhl)
2312 i = phi_mesh%jj(ni,m2)
2313 ib = la_phi%loc_to_glob(1,i)
2317 j = h_mesh%jj(nj,m1)
2319 jb = la_h%loc_to_glob(kj,j)
2320 jx = (kj-1)*n_w1 + nj
2323 mat_loc1(ix,jx) = sij(1,ni,nj)
2324 mat_loc2(ix,jx) = sij(1,ni,nj)
2325 ELSEIF (kj == 3)
THEN
2326 mat_loc1(ix,jx) = sij(5,ni,nj)
2327 mat_loc2(ix,jx) = sij(5,ni,nj)
2332 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 3*n_w1, jdxn(1:3*n_w1), &
2333 mat_loc1(1:n_w2,1:3*n_w1), add_values, ierr)
2334 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 3*n_w1, jdxn(1:3*n_w1), &
2335 mat_loc2(1:n_w2,1:3*n_w1), add_values, ierr)
2341 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2344 x = rjs(ls,ms2)*ray/sum(sigma_np(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2346 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2347 drmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(1,:,ls,ms1))
2348 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(2,:,ls,ms1))
2351 y = x*(dw_s(2,ni,ls,ms2)*rnorms(1,ls,ms2) - dw_s(1,ni,ls,ms2)*rnorms(2,ls,ms2))
2353 sij(1,ni,nj) = sij(1,ni,nj) + y *( (-1/muhl**2)*dzmuhl)*w_cs(nj,ls)
2354 sij(5,ni,nj) = sij(5,ni,nj) + y *(-(-1/muhl**2)*drmuhl)*w_cs(nj,ls)
2366 i = phi_mesh%jj(ni,m2)
2367 ib = la_phi%loc_to_glob(1,i)
2372 j = interface_h_phi%jjs1(nj,ms)
2373 jb = la_h%loc_to_glob(kj,j)
2374 jx = (kj-1)*n_ws1 + nj
2377 mat_loc1(ix,jx) = sij(1,ni,nj)
2378 mat_loc2(ix,jx) = sij(1,ni,nj)
2379 ELSEIF (kj == 3)
THEN
2380 mat_loc1(ix,jx) = sij(5,ni,nj)
2381 mat_loc2(ix,jx) = sij(5,ni,nj)
2386 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
2387 mat_loc1(1:n_w2,1:3*n_ws1), add_values, ierr)
2388 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 3*n_ws1, jdxn(1:3*n_ws1), &
2389 mat_loc2(1:n_w2,1:3*n_ws1), add_values, ierr)
2398 i = h_mesh%jj(ni,m1)
2399 ib = la_h%loc_to_glob(ki,i)
2400 ix = (ki-1)*n_w1 + ni
2403 j = phi_mesh%jj(nj,m2)
2404 jb = la_phi%loc_to_glob(1,j)
2408 mat_loc1(ix,jx) = sij(1,nj,ni)
2409 mat_loc2(ix,jx) = sij(1,nj,ni)
2410 ELSEIF (ki == 3)
THEN
2411 mat_loc1(ix,jx) = sij(5,nj,ni)
2412 mat_loc2(ix,jx) = sij(5,nj,ni)
2417 CALL matsetvalues(h_p_phi_mat1, 3*n_w1, idxn(1:3*n_w1), n_w2, jdxn(1:n_w2), &
2418 mat_loc1(1:3*n_w1,1:n_w2), add_values, ierr)
2419 CALL matsetvalues(h_p_phi_mat2, 3*n_w1, idxn(1:3*n_w1), n_w2, jdxn(1:n_w2), &
2420 mat_loc2(1:3*n_w1,1:n_w2), add_values, ierr)
2427 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* phi_mesh%gauss%wws(:,ls))
2428 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2431 x = c_lap*rjs(ls,ms2)*ray
2435 sij(1,ni,nj) = sij(1,ni,nj) - x*w_cs(nj,ls)*wws(ni,ls)*rnorms(1,ls,ms2)
2436 sij(5,ni,nj) = sij(5,ni,nj) - x*w_cs(nj,ls)*wws(ni,ls)*rnorms(2,ls,ms2)
2447 i = interface_h_phi%jjs2(ni,ms)
2448 ib = la_phi%loc_to_glob(1,i)
2452 j = interface_h_phi%jjs1(nj,ms)
2453 jb = la_h%loc_to_glob(1,j)
2456 mat_loc1(ix,jx) = sij(1,ni,nj)
2457 mat_loc2(ix,jx) = sij(1,ni,nj)
2459 jb = la_h%loc_to_glob(3,j)
2462 mat_loc1(ix,jx) = sij(5,ni,nj)
2463 mat_loc2(ix,jx) = sij(5,ni,nj)
2466 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), 2*n_ws1, jdxn(1:2*n_ws1), &
2467 mat_loc1(1:n_ws2,1:2*n_ws1), add_values, ierr)
2468 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), 2*n_ws1, jdxn(1:2*n_ws1), &
2469 mat_loc2(1:n_ws2,1:2*n_ws1), add_values, ierr)
2480 IF (stab(2) > 1.d-12)
THEN
2489 ms2 = interface_h_phi%mesh2(ms)
2490 hm1 = sum(rjs(:,ms2))**(2*alpha-1)
2495 muhl = sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
2498 DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,ms2)
2499 ray = ray + phi_mesh%rr(1,i)* phi_mesh%gauss%wws(ni,ls)
2504 ray = stab_div*ray*hm1*rjs(ls,ms2)
2512 x = muhl*w_cs(ni,ls)*w_cs(nj,ls)*ray
2514 hsij(1,ni,nj) = hsij(1,ni,nj) + x*rnorms(1,ls,ms2)**2
2515 hsij(4,ni,nj) = hsij(4,ni,nj) + x*rnorms(1,ls,ms2)*rnorms(2,ls,ms2)
2516 hsij(6,ni,nj) = hsij(6,ni,nj) + x*rnorms(2,ls,ms2)**2
2520 x = muhl*mu_phi*w_cs(ni,ls)*(dw_s(1,nj,ls,ms2)*rnorms(1,ls,ms2) +&
2521 dw_s(2,nj,ls,ms2)*rnorms(2,ls,ms2))*ray
2522 sij(1,ni,nj) = sij(1,ni,nj) - x*rnorms(1,ls,ms2)
2523 sij(5,ni,nj) = sij(5,ni,nj) - x*rnorms(2,ls,ms2)
2530 x = mu_phi*w_cs(nj,ls)*(dw_s(1,ni,ls,ms2)*rnorms(1,ls,ms2) +&
2531 dw_s(2,ni,ls,ms2)*rnorms(2,ls,ms2))*ray
2532 smuij(1,ni,nj) = smuij(1,ni,nj) - x*rnorms(1,ls,ms2)
2533 smuij(5,ni,nj) = smuij(5,ni,nj) - x*rnorms(2,ls,ms2)
2541 x = mu_phi**2*(dw_s(1,ni,ls,ms2)*rnorms(1,ls,ms2) + dw_s(2,ni,ls,ms2)*rnorms(2,ls,ms2))* &
2542 (dw_s(1,nj,ls,ms2)*rnorms(1,ls,ms2) + dw_s(2,nj,ls,ms2)*rnorms(2,ls,ms2))*ray
2543 phisij(ni,nj) = phisij(ni,nj) + x
2548 sij(2,:,:) = sij(1,:,:)
2549 sij(6,:,:) = sij(5,:,:)
2555 i = h_mesh%jjs(ni,ms1)
2557 ib = la_h%loc_to_glob(ki,i)
2558 ix = (ki/2)*n_ws1 + ni
2561 j = h_mesh%jjs(nj,ms1)
2563 jb = la_h%loc_to_glob(kj,j)
2564 jx = (kj/2)*n_ws1 + nj
2567 mat_loc1(ix,jx) = hsij(1,ni,nj)
2568 mat_loc2(ix,jx) = hsij(1,ni,nj)
2569 ELSE IF (ki*kj==9)
THEN
2570 mat_loc1(ix,jx) = hsij(6,ni,nj)
2571 mat_loc2(ix,jx) = hsij(6,ni,nj)
2572 ELSE IF (ki*kj==3)
THEN
2573 mat_loc1(ix,jx) = hsij(4,ni,nj)
2574 mat_loc2(ix,jx) = hsij(4,ni,nj)
2580 j = phi_mesh%jj(nj,m2)
2581 jb = la_phi%loc_to_glob(1,j)
2584 mat_loc1(ix,jx) = sij(2*ki-1,ni,nj)
2585 mat_loc2(ix,jx) = sij(2*ki-1,ni,nj)
2589 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), &
2590 mat_loc1(1:2*n_ws1,1:2*n_ws1+n_w2), add_values, ierr)
2591 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), &
2592 mat_loc2(1:2*n_ws1,1:2*n_ws1+n_w2), add_values, ierr)
2597 i = phi_mesh%jj(ni,m2)
2598 ib = la_phi%loc_to_glob(1,i)
2602 j = h_mesh%jjs(nj,ms1)
2604 jb = la_h%loc_to_glob(kj,j)
2605 jx = (kj/2)*n_ws1 + nj
2610 mat_loc1(ix,jx) = smuij(2*kj-1,ni,nj)
2611 mat_loc2(ix,jx) = smuij(2*kj-1,ni,nj)
2617 j = phi_mesh%jj(nj,m2)
2618 jb = la_phi%loc_to_glob(1,j)
2621 mat_loc1(ix,jx) = phisij(ni,nj)
2622 mat_loc2(ix,jx) = phisij(ni,nj)
2625 CALL matsetvalues(h_p_phi_mat1, n_w2, idxn(1:n_w2), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2626 mat_loc1(1:n_w2,1:2*n_ws1+n_w2), add_values, ierr)
2627 CALL matsetvalues(h_p_phi_mat2, n_w2, idxn(1:n_w2), 2*n_ws1+n_w2, jdxn(1:2*n_ws1+n_w2), &
2628 mat_loc2(1:n_w2,1:2*n_ws1+n_w2), add_values, ierr)
2639 IF (.NOT.present(index_fourier) .OR. .NOT.present(r_fourier))
RETURN
2640 IF (r_fourier.GT.0.d0)
THEN
2642 DO ms = 1, phi_mesh%mes
2643 IF (phi_mesh%sides(ms) /= index_fourier) cycle
2647 DO ls = 1, phi_mesh%gauss%l_Gs
2650 ray = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms))* phi_mesh%gauss%wws(:,ls))
2652 x = c_mu_phi*rjs(ls,ms)*ray/r_fourier
2654 DO ni=1, phi_mesh%gauss%n_ws
2655 DO nj=1, phi_mesh%gauss%n_ws
2656 phisij(ni,nj) = phisij(ni,nj) + x*wws(ni,ls)*wws(nj,ls)
2663 DO ni = 1, phi_mesh%gauss%n_ws
2664 i = phi_mesh%jjs(ni,ms)
2665 ib = la_phi%loc_to_glob(1,i)
2667 DO nj = 1, phi_mesh%gauss%n_ws
2668 j = phi_mesh%jjs(nj,ms)
2669 jb = la_phi%loc_to_glob(1,j)
2673 CALL matsetvalues(h_p_phi_mat1, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
2674 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
2675 CALL matsetvalues(h_p_phi_mat2, n_ws2, idxn(1:n_ws2), n_ws2, jdxn(1:n_ws2), &
2676 phisij(1:n_ws2,1:n_ws2), add_values, ierr)
2680 CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
2681 CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
2682 CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
2683 CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
2690 mode, stab, mu_h_field, la_h, h_p_phi_mat1, h_p_phi_mat2, sigma_np)
2699 INTEGER,
DIMENSION(:),
INTENT(IN) :: dirichlet_bdy_h_sides
2700 INTEGER,
INTENT(IN) :: mode
2701 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
2702 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field, sigma_np
2703 INTEGER :: ms, ls, ni, nj, i, j, &
2704 n_ws1, n_w1, m1, ki, kj, ib, jb
2705 REAL(KIND=8) :: x, y, hm1
2706 REAL(KIND=8) :: ray, error, stab_colle_h_mu
2707 REAL(KIND=8),
DIMENSION(9,H_mesh%gauss%n_w,H_mesh%gauss%n_w) :: hsij
2727 REAL(KIND=8) :: muhl, dzmuhl, drmuhl
2729 REAL(KIND=8) :: c_sym=.0d0
2734 REAL(KIND=8),
DIMENSION(3*H_mesh%gauss%n_w,3*H_mesh%gauss%n_w) :: mat_loc1, mat_loc2
2735 INTEGER ,
DIMENSION(3*H_mesh%gauss%n_w) :: idxn, jdxn
2741 REAL(KIND=8),
DIMENSION(2) :: dmu_field
2742 REAL(KIND=8),
DIMENSION(2) :: gauss_pt
2743 INTEGER,
DIMENSION(1) :: gauss_pt_id
2744 REAL(KIND=8),
DIMENSION(1) :: dummy_mu_bar
2747 #include "petsc/finclude/petsc.h"
2748 petscerrorcode :: ierr
2749 mat :: h_p_phi_mat1, h_p_phi_mat2
2752 stab_colle_h_mu = stab(3)
2759 n_ws1 = h_mesh%gauss%n_ws
2760 n_w1 = h_mesh%gauss%n_w
2768 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
2769 ms = dirichlet_bdy_h_sides(count)
2770 hm1 = stab_colle_h_mu/sum(h_mesh%gauss%rjs(:,ms))
2771 m1 = h_mesh%neighs(ms)
2772 mesh_id = h_mesh%i_d(m1)
2781 DO ls = 1, h_mesh%gauss%l_Gs
2783 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2786 IF(inputs%if_use_fem_integration_for_mu_bar)
THEN
2787 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2789 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2790 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2791 gauss_pt_id(1)=mesh_id
2793 muhl=dummy_mu_bar(1)
2797 x = hm1*h_mesh%gauss%rjs(ls,ms)*ray /muhl
2799 DO ni = 1, h_mesh%gauss%n_ws
2800 DO nj = 1, h_mesh%gauss%n_ws
2801 y = x * h_mesh%gauss%wws(ni,ls)*h_mesh%gauss%wws(nj,ls)
2803 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(h_mesh%gauss%rnorms(2,ls,ms)**2)
2804 hsij(4,ni,nj) = hsij(4,ni,nj) - y*h_mesh%gauss%rnorms(1,ls,ms)*h_mesh%gauss%rnorms(2,ls,ms)
2805 hsij(5,ni,nj) = hsij(5,ni,nj) + y
2806 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(h_mesh%gauss%rnorms(1,ls,ms)**2)
2821 i = h_mesh%jjs(ni,ms)
2822 ib = la_h%loc_to_glob(ki,i)
2823 ix = ni + (ki-1)*n_ws1
2827 j = h_mesh%jjs(nj,ms)
2828 jb = la_h%loc_to_glob(kj,j)
2829 jx = nj + (kj-1)*n_ws1
2831 IF ((ki == 1) .AND. (kj == 1))
THEN
2832 mat_loc1(ix,jx) = hsij(1,ni,nj)
2833 mat_loc2(ix,jx) = hsij(1,ni,nj)
2834 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN
2835 mat_loc1(ix,jx) = hsij(4,ni,nj)
2836 mat_loc2(ix,jx) = hsij(4,ni,nj)
2837 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN
2838 mat_loc1(ix,jx) = hsij(4,nj,ni)
2839 mat_loc2(ix,jx) = hsij(4,nj,ni)
2840 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
2841 mat_loc1(ix,jx) = hsij(5,ni,nj)
2842 mat_loc2(ix,jx) = hsij(5,ni,nj)
2843 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN
2844 mat_loc1(ix,jx) = hsij(6,ni,nj)
2845 mat_loc2(ix,jx) = hsij(6,ni,nj)
2852 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2853 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2854 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2855 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2866 DO ls = 1, h_mesh%gauss%l_Gs
2868 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
2871 IF(inputs%if_use_fem_integration_for_mu_bar)
THEN
2872 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2873 drmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_s(1,:,ls,ms))
2874 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m1))*dw_s(2,:,ls,ms))
2876 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
2877 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
2878 gauss_pt_id(1)=mesh_id
2880 muhl=dummy_mu_bar(1)
2882 drmuhl =dmu_field(1)
2883 dzmuhl =dmu_field(2)
2887 x = h_mesh%gauss%rjs(ls,ms)*ray/sum(sigma_np(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
2894 y = x*h_mesh%gauss%wws(ni,ls)*h_mesh%gauss%wws(nj,ls)
2896 hsij(2,ni,nj) = hsij(2,ni,nj) + y * (-mode/ray)*(rnorms(1,ls,ms))/muhl
2897 hsij(3,ni,nj) = hsij(3,ni,nj) + y * mode/ray *(rnorms(1,ls,ms))/muhl
2898 hsij(5,ni,nj) = hsij(5,ni,nj) + y * (-1/ray) *(rnorms(1,ls,ms))/muhl
2899 hsij(8,ni,nj) = hsij(8,ni,nj) + y * (-mode/ray)*(rnorms(2,ls,ms))/muhl
2900 hsij(9,ni,nj) = hsij(9,ni,nj) + y * mode/ray *(rnorms(2,ls,ms))/muhl
2902 hsij(1,ni,nj) = hsij(1,ni,nj) - y*(rnorms(2,ls,ms))*(-(dzmuhl/muhl**2))
2903 hsij(4,ni,nj) = hsij(4,ni,nj) - y*(rnorms(2,ls,ms))*( (drmuhl/muhl**2))
2904 hsij(5,ni,nj) = hsij(5,ni,nj) &
2905 + y*(rnorms(1,ls,ms)*drmuhl+rnorms(2,ls,ms)*dzmuhl)/muhl**2
2906 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(rnorms(1,ls,ms))*( (drmuhl/muhl**2))
2907 hsij(7,ni,nj) = hsij(7,ni,nj) + y*(rnorms(1,ls,ms))*(-(dzmuhl/muhl**2))
2924 i = h_mesh%jjs(ni,ms)
2925 ib = la_h%loc_to_glob(ki,i)
2926 ix = ni + (ki-1)*n_ws1
2930 j = h_mesh%jjs(nj,ms)
2931 jb = la_h%loc_to_glob(kj,j)
2932 jx = nj + (kj-1)*n_ws1
2934 IF ((ki == 1) .AND. (kj == 1))
THEN
2935 mat_loc1(ix,jx) = hsij(1,ni,nj)
2936 mat_loc2(ix,jx) = hsij(1,ni,nj)
2937 ELSE IF ( (ki == 1) .AND. (kj == 3))
THEN
2938 mat_loc1(ix,jx) = hsij(4,ni,nj)
2939 mat_loc2(ix,jx) = hsij(4,ni,nj)
2940 ELSE IF ( (ki == 2) .AND. (kj == 1))
THEN
2941 mat_loc1(ix,jx) = hsij(2,ni,nj)
2942 mat_loc2(ix,jx) = hsij(3,ni,nj)
2943 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
2944 mat_loc1(ix,jx) = hsij(5,ni,nj)
2945 mat_loc2(ix,jx) = hsij(5,ni,nj)
2946 ELSEIF ( (ki == 2) .AND. (kj == 3))
THEN
2947 mat_loc1(ix,jx) = hsij(8,ni,nj)
2948 mat_loc2(ix,jx) = hsij(9,ni,nj)
2949 ELSEIF ( (ki == 3) .AND. (kj == 1))
THEN
2950 mat_loc1(ix,jx) = hsij(7,ni,nj)
2951 mat_loc2(ix,jx) = hsij(7,ni,nj)
2952 ELSEIF ( (ki == 3) .AND. (kj == 3))
THEN
2953 mat_loc1(ix,jx) = hsij(6,ni,nj)
2954 mat_loc2(ix,jx) = hsij(6,ni,nj)
2961 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2962 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2963 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2964 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2972 i = h_mesh%jjs(ni,ms)
2973 ib = la_h%loc_to_glob(ki,i)
2974 ix = ni + (ki-1)*n_ws1
2978 j = h_mesh%jjs(nj,ms)
2979 jb = la_h%loc_to_glob(kj,j)
2980 jx = nj + (kj-1)*n_ws1
2982 IF ( (kj == 2) .AND. (ki == 1))
THEN
2983 mat_loc1(ix,jx) = hsij(2,nj,ni)
2984 mat_loc2(ix,jx) = hsij(3,nj,ni)
2985 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN
2986 mat_loc1(ix,jx) = hsij(5,nj,ni)
2987 mat_loc2(ix,jx) = hsij(5,nj,ni)
2988 ELSEIF ( (kj == 2) .AND. (ki == 3))
THEN
2989 mat_loc1(ix,jx) = hsij(8,nj,ni)
2990 mat_loc2(ix,jx) = hsij(9,nj,ni)
2996 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2997 mat_loc1(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
2998 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_ws1, jdxn(1:3*n_ws1), &
2999 mat_loc2(1:3*n_ws1,1:3*n_ws1), add_values, ierr)
3005 DO ls = 1, h_mesh%gauss%l_Gs
3007 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3010 IF(inputs%if_use_fem_integration_for_mu_bar)
THEN
3011 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3013 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3014 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jjs(:,ms))* h_mesh%gauss%wws(:,ls))
3015 gauss_pt_id(1)=mesh_id
3017 muhl=dummy_mu_bar(1)
3021 x = h_mesh%gauss%rjs(ls,ms)*ray/(sum(sigma_np(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))*muhl)
3029 y = x*h_mesh%gauss%wws(ni,ls)
3031 hsij(1,ni,nj) = hsij(1,ni,nj) + y*(-h_mesh%gauss%dw_s(2,nj,ls,ms))*(rnorms(2,ls,ms))
3032 hsij(4,ni,nj) = hsij(4,ni,nj) + y* h_mesh%gauss%dw_s(1,nj,ls,ms) *(rnorms(2,ls,ms))
3033 hsij(5,ni,nj) = hsij(5,ni,nj) + &
3034 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)))
3035 hsij(6,ni,nj) = hsij(6,ni,nj) + y*(-h_mesh%gauss%dw_s(1,nj,ls,ms))*(rnorms(1,ls,ms))
3036 hsij(7,ni,nj) = hsij(7,ni,nj) + y* h_mesh%gauss%dw_s(2,nj,ls,ms) *(rnorms(1,ls,ms))
3050 i = h_mesh%jjs(ni,ms)
3051 ib = la_h%loc_to_glob(ki,i)
3052 ix = ni + (ki-1)*n_ws1
3056 j = h_mesh%jj(nj,m1)
3057 jb = la_h%loc_to_glob(kj,j)
3058 jx = nj + (kj-1)*n_w1
3060 IF ((ki == 1) .AND. (kj == 1))
THEN
3061 mat_loc1(ix,jx) = hsij(1,ni,nj)
3062 mat_loc2(ix,jx) = hsij(1,ni,nj)
3063 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN
3064 mat_loc1(ix,jx) = hsij(4,ni,nj)
3065 mat_loc2(ix,jx) = hsij(4,ni,nj)
3066 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
3067 mat_loc1(ix,jx) = hsij(5,ni,nj)
3068 mat_loc2(ix,jx) = hsij(5,ni,nj)
3069 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN
3070 mat_loc1(ix,jx) = hsij(6,ni,nj)
3071 mat_loc2(ix,jx) = hsij(6,ni,nj)
3072 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN
3073 mat_loc1(ix,jx) = hsij(7,ni,nj)
3074 mat_loc2(ix,jx) = hsij(7,ni,nj)
3081 CALL matsetvalues(h_p_phi_mat1, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
3082 mat_loc1(1:3*n_ws1,1:3*n_w1), add_values, ierr)
3083 CALL matsetvalues(h_p_phi_mat2, 3*n_ws1 , idxn(1:3*n_ws1), 3*n_w1, jdxn(1:3*n_w1), &
3084 mat_loc2(1:3*n_ws1,1:3*n_w1), add_values, ierr)
3092 i = h_mesh%jj(ni,m1)
3093 ib = la_h%loc_to_glob(ki,i)
3094 ix = ni + (ki-1)*n_w1
3098 j = h_mesh%jjs(nj,ms)
3099 jb = la_h%loc_to_glob(kj,j)
3100 jx = nj + (kj-1)*n_ws1
3102 IF ((kj == 1) .AND. (ki == 1))
THEN
3103 mat_loc1(ix,jx) = hsij(1,nj,ni)
3104 mat_loc2(ix,jx) = hsij(1,nj,ni)
3105 ELSEIF ((kj == 1) .AND. (ki == 3))
THEN
3106 mat_loc1(ix,jx) = hsij(4,nj,ni)
3107 mat_loc2(ix,jx) = hsij(4,nj,ni)
3108 ELSEIF ((kj == 2) .AND. (ki == 2))
THEN
3109 mat_loc1(ix,jx) = hsij(5,nj,ni)
3110 mat_loc2(ix,jx) = hsij(5,nj,ni)
3111 ELSEIF ((kj == 3) .AND. (ki == 3))
THEN
3112 mat_loc1(ix,jx) = hsij(6,nj,ni)
3113 mat_loc2(ix,jx) = hsij(6,nj,ni)
3114 ELSEIF ((kj == 3) .AND. (ki == 1))
THEN
3115 mat_loc1(ix,jx) = hsij(7,nj,ni)
3116 mat_loc2(ix,jx) = hsij(7,nj,ni)
3123 CALL matsetvalues(h_p_phi_mat1, 3*n_w1 , idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
3124 mat_loc1(1:3*n_w1,1:3*n_ws1), add_values, ierr)
3125 CALL matsetvalues(h_p_phi_mat2, 3*n_w1 , idxn(1:3*n_w1), 3*n_ws1, jdxn(1:3*n_ws1), &
3126 mat_loc2(1:3*n_w1,1:3*n_ws1), add_values, ierr)
3130 CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
3131 CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
3132 CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
3133 CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
3143 SUBROUTINE courant_int_by_parts(H_mesh,phi_mesh,interface_H_phi,sigma,mu_phi,mu_H_field,time,mode,&
3144 rhs_h,nl, la_h, la_phi, vb_1, vb_2, b_ext, h_pert, sigma_curl, j_over_sigma)
3155 TYPE(mesh_type),
INTENT(IN) :: h_mesh, phi_mesh
3157 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
3158 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma
3159 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field
3160 INTEGER,
INTENT(IN) :: mode
3161 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: nl
3162 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: b_ext
3163 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: rhs_h
3164 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: h_pert
3165 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: j_over_sigma
3166 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: sigma_curl
3168 REAL(KIND=8),
DIMENSION(6) :: j_over_sigma_gauss
3170 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: src_h
3171 REAL(KIND=8),
DIMENSION(phi_mesh%np,2) :: src_phi
3172 REAL(KIND=8),
DIMENSION(6) :: c
3173 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_ws,phi_mesh%gauss%l_Gs) :: w_cs
3174 REAL(KIND=8),
DIMENSION(2) :: gaussp
3176 INTEGER :: m, l, i, ni, k, ms, ls, n_ws1, n_ws2, ms1, ms2, h_bloc_size, n_w2, m1
3178 REAL(KIND=8),
DIMENSION(6) :: jsolh_anal, rhs_hl
3179 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%n_w) :: dwh
3182 REAL(KIND=8) :: ray_rjl, muhl, drmuhl, dzmuhl
3187 INTEGER,
DIMENSION(H_mesh%np) :: idxn_h
3188 INTEGER,
DIMENSION(phi_mesh%np) :: idxn_phi
3190 REAL(KIND=8),
DIMENSION(2) :: dmu_field
3191 REAL(KIND=8),
DIMENSION(2) :: gauss_pt
3192 INTEGER,
DIMENSION(1) :: gauss_pt_id
3193 REAL(KIND=8),
DIMENSION(1) :: dummy_mu_bar
3196 REAL(KIND=8),
DIMENSION(6) :: b_ext_l
3197 #include "petsc/finclude/petsc.h"
3198 petscerrorcode :: ierr
3203 CALL veczeroentries(vb_1, ierr)
3204 CALL veczeroentries(vb_2, ierr)
3216 mesh_id1 = h_mesh%i_d(m)
3217 DO l = 1, h_mesh%gauss%l_G
3222 muhl=sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
3224 dwh = h_mesh%gauss%dw(:,:,l,m)
3227 b_ext_l(k) = sum(b_ext(h_mesh%jj(:,m),k)*h_mesh%gauss%ww(:,l))
3233 j_over_sigma_gauss = 0.d0
3234 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
3235 gaussp = gaussp + h_mesh%rr(:,i)*h_mesh%gauss%ww(ni,l)
3236 jsolh_anal(:) = jsolh_anal(:) + nl(i,:)*h_mesh%gauss%ww(ni,l)
3237 rhs_hl(:) = rhs_hl(:) + rhs_h(i,:)*h_mesh%gauss%ww(ni,l)
3238 j_over_sigma_gauss(:) = j_over_sigma_gauss(:) + j_over_sigma(i,:)*h_mesh%gauss%ww(ni,l)
3241 ray_rjl = h_mesh%gauss%rj(l,m)*ray
3243 IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid)
THEN
3245 jsolh_anal(k) = jsolh_anal(k) + j_over_sigma_gauss(k) + sigma_curl(index,k)
3250 jsolh_anal(k) = jsolh_anal(k) + &
3251 jexact_gauss(k, gaussp, mode, mu_phi, sigma(m), muhl, time, mesh_id1, b_ext_l)/sigma(m)
3256 IF (inputs%if_permeability_variable_in_theta)
THEN
3258 IF(inputs%if_use_fem_integration_for_mu_bar)
THEN
3259 muhl = sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
3260 drmuhl = sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%dw(1,:,l,m))
3261 dzmuhl = sum(mu_h_field(h_mesh%jj(:,m))*h_mesh%gauss%dw(2,:,l,m))
3263 gauss_pt(1)=sum(h_mesh%rr(1,h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
3264 gauss_pt(2)=sum(h_mesh%rr(2,h_mesh%jj(:,m))*h_mesh%gauss%ww(:,l))
3265 gauss_pt_id(1)=mesh_id1
3267 muhl=dummy_mu_bar(1)
3269 drmuhl =dmu_field(1)
3270 dzmuhl =dmu_field(2)
3274 DO ni = 1,h_mesh%gauss%n_w
3280 c(1) = c(1) + ( mode/ray*h_pert(i,6)*h_mesh%gauss%ww(ni,l) &
3281 - h_pert(i,3)*h_mesh%gauss%dw(2,ni,l,m)) &
3282 + (1/muhl)*( mode/ray*b_ext(i,6)*h_mesh%gauss%ww(ni,l) &
3283 - b_ext(i,3)*h_mesh%gauss%dw(2,ni,l,m)) &
3284 -(1/muhl**2)*(-dzmuhl*b_ext(i,3)*h_mesh%gauss%ww(ni,l))
3286 c(2) = c(2) + (-mode/ray*h_pert(i,5)*h_mesh%gauss%ww(ni,l) &
3287 - h_pert(i,4)*h_mesh%gauss%dw(2,ni,l,m)) &
3288 + (1/muhl)*(-mode/ray*b_ext(i,5)*h_mesh%gauss%ww(ni,l) &
3289 - b_ext(i,4)*h_mesh%gauss%dw(2,ni,l,m)) &
3290 -(1/muhl**2)*(-dzmuhl*b_ext(i,4)*h_mesh%gauss%ww(ni,l))
3293 c(3) = c(3) + (h_pert(i,1)*h_mesh%gauss%dw(2,ni,l,m) &
3294 - h_pert(i,5)*h_mesh%gauss%dw(1,ni,l,m)) &
3295 + (1/muhl)*(b_ext(i,1)*h_mesh%gauss%dw(2,ni,l,m) &
3296 - b_ext(i,5)*h_mesh%gauss%dw(1,ni,l,m)) &
3297 -(1/muhl**2)*(dzmuhl*b_ext(i,1)-drmuhl*b_ext(i,5))*h_mesh%gauss%ww(ni,l)
3299 c(4) = c(4) + (h_pert(i,2)*h_mesh%gauss%dw(2,ni,l,m) &
3300 - h_pert(i,6)*h_mesh%gauss%dw(1,ni,l,m)) &
3301 + (1/muhl)*(b_ext(i,2)*h_mesh%gauss%dw(2,ni,l,m) &
3302 - b_ext(i,6)*h_mesh%gauss%dw(1,ni,l,m)) &
3303 -(1/muhl**2)*(dzmuhl*b_ext(i,2)-drmuhl*b_ext(i,6))*h_mesh%gauss%ww(ni,l)
3306 c(5) = c(5) + (h_pert(i,3)*h_mesh%gauss%dw(1,ni,l,m) &
3307 + h_pert(i,3)*h_mesh%gauss%ww(ni,l)/ray &
3308 - mode/ray*h_pert(i,2)*h_mesh%gauss%ww(ni,l)) &
3309 + (1/muhl)*(b_ext(i,3)*h_mesh%gauss%dw(1,ni,l,m) &
3310 + b_ext(i,3)*h_mesh%gauss%ww(ni,l)/ray &
3311 - mode/ray*b_ext(i,2)*h_mesh%gauss%ww(ni,l)) &
3312 -(1/muhl**2)*(drmuhl*b_ext(i,3)*h_mesh%gauss%ww(ni,l))
3314 c(6) = c(6) + (h_pert(i,4)*h_mesh%gauss%dw(1,ni,l,m) &
3315 + h_pert(i,4)*h_mesh%gauss%ww(ni,l)/ray &
3316 + mode/ray*h_pert(i,1)*h_mesh%gauss%ww(ni,l)) &
3317 + (1/muhl)*(b_ext(i,4)*h_mesh%gauss%dw(1,ni,l,m) &
3318 + b_ext(i,4)*h_mesh%gauss%ww(ni,l)/ray &
3319 + mode/ray*b_ext(i,1)*h_mesh%gauss%ww(ni,l)) &
3320 -(1/muhl**2)*(drmuhl*b_ext(i,4)*h_mesh%gauss%ww(ni,l))
3324 jsolh_anal= jsolh_anal + c
3328 DO ni = 1,h_mesh%gauss%n_w
3333 src_h(i,1) = src_h(i,1)+ ray_rjl &
3334 *(jsolh_anal(3)*dwh(2,ni) &
3335 + mode/ray*jsolh_anal(6)*h_mesh%gauss%ww(ni,l) &
3336 + rhs_hl(1)*h_mesh%gauss%ww(ni,l))
3338 src_h(i,2) = src_h(i,2)+ ray_rjl &
3339 *(jsolh_anal(4)*dwh(2,ni) &
3340 - mode/ray*jsolh_anal(5)*h_mesh%gauss%ww(ni,l) &
3341 + rhs_hl(2)*h_mesh%gauss%ww(ni,l))
3344 src_h(i,3) = src_h(i,3)+ ray_rjl &
3345 * (-jsolh_anal(1)*dwh(2,ni) &
3346 + 1/ray*jsolh_anal(5)*(ray*dwh(1,ni) + h_mesh%gauss%ww(ni,l)) &
3347 + rhs_hl(3)*h_mesh%gauss%ww(ni,l))
3349 src_h(i,4) = src_h(i,4)+ ray_rjl &
3350 * (-jsolh_anal(2)*dwh(2,ni) &
3351 + 1/ray*jsolh_anal(6)*(ray*dwh(1,ni) + h_mesh%gauss%ww(ni,l)) &
3352 + rhs_hl(4)*h_mesh%gauss%ww(ni,l))
3355 src_h(i,5) = src_h(i,5)+ ray_rjl* &
3356 (-mode/ray*jsolh_anal(2)*h_mesh%gauss%ww(ni,l) &
3357 - jsolh_anal(3)*dwh(1,ni) &
3358 + rhs_hl(5)*h_mesh%gauss%ww(ni,l))
3360 src_h(i,6) = src_h(i,6)+ ray_rjl* &
3361 (mode/ray*jsolh_anal(1)*h_mesh%gauss%ww(ni,l) &
3362 - jsolh_anal(4)*dwh(1,ni) &
3363 + rhs_hl(6)*h_mesh%gauss%ww(ni,l))
3420 CALL
gauss(phi_mesh)
3422 n_ws1 = h_mesh%gauss%n_ws
3423 n_ws2 = phi_mesh%gauss%n_ws
3424 n_w2 = phi_mesh%gauss%n_w
3426 h_bloc_size = h_mesh%np
3428 IF (interface_h_phi%mes /=0)
THEN
3429 IF (h_mesh%gauss%n_ws == n_ws)
THEN
3433 w_cs(1,ls)= wws(1,ls)+0.5*wws(3,ls)
3434 w_cs(2,ls)= wws(2,ls)+0.5*wws(3,ls)
3442 DO ms = 1, interface_h_phi%mes
3444 ms2 = interface_h_phi%mesh2(ms)
3445 ms1 = interface_h_phi%mesh1(ms)
3446 m = phi_mesh%neighs(ms2)
3447 m1 = h_mesh%neighs(ms1)
3448 mesh_id1 = h_mesh%i_d(m1)
3452 muhl=sum(mu_h_field(interface_h_phi%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
3455 b_ext_l(k) = sum(b_ext(interface_h_phi%jjs1(1:n_ws1,ms),k)*w_cs(1:n_ws1,ls))
3456 j_over_sigma_gauss(k) = sum(j_over_sigma(interface_h_phi%jjs1(1:n_ws1,ms),k)*w_cs(1:n_ws1,ls))
3461 DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,interface_h_phi%mesh2(ms))
3462 ray = ray + phi_mesh%rr(1,i)* wws(ni,ls)
3467 i=phi_mesh%jjs(ni,ms2)
3468 gaussp = gaussp + phi_mesh%rr(:,i)*phi_mesh%gauss%wws(ni,ls)
3472 jsolh_anal(k) =
jexact_gauss(k, gaussp, mode, mu_phi ,sigma(m1), &
3473 muhl, time, mesh_id1, b_ext_l)/sigma(m1) &
3474 + sum(nl(h_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
3494 i = interface_h_phi%jjs1(ni,ms)
3495 src_h(i,1) = src_h(i,1)+rjs(ls,ms2)*ray*( &
3496 -jsolh_anal(3)*w_cs(ni,ls)*(-rnorms(2,ls,ms2)))
3498 src_h(i,2) = src_h(i,2)+rjs(ls,ms2)*ray*( &
3499 -jsolh_anal(4)*w_cs(ni,ls)*(-rnorms(2,ls,ms2)))
3501 src_h(i,3) = src_h(i,3)+rjs(ls,ms2)*ray*( &
3502 jsolh_anal(1)*w_cs(ni,ls)*(-rnorms(2,ls,ms2)) &
3503 -jsolh_anal(5)*w_cs(ni,ls)*(-rnorms(1,ls,ms2)))
3505 src_h(i,4) = src_h(i,4)+rjs(ls,ms2)*ray*( &
3506 jsolh_anal(2)*w_cs(ni,ls)*(-rnorms(2,ls,ms2)) &
3507 -jsolh_anal(6)*w_cs(ni,ls)*(-rnorms(1,ls,ms2)))
3509 src_h(i,5) = src_h(i,5)+rjs(ls,ms2)*ray*( &
3510 jsolh_anal(3)*w_cs(ni,ls)*(-rnorms(1,ls,ms2)))
3512 src_h(i,6) = src_h(i,6)+rjs(ls,ms2)*ray*( &
3513 jsolh_anal(4)*w_cs(ni,ls)*(-rnorms(1,ls,ms2)))
3519 i = interface_h_phi%jjs2(ni,ms)
3522 src_phi(i,1) = src_phi(i,1)+rjs(ls,ms2)*( &
3523 - mode*jsolh_anal(2)*wws(ni,ls) * rnorms(2,ls,ms2) &
3524 + mode*jsolh_anal(6)*wws(ni,ls) * rnorms(1,ls,ms2))
3526 src_phi(i,2) = src_phi(i,2)+rjs(ls,ms2)*( &
3527 + mode*jsolh_anal(1)*wws(ni,ls) * rnorms(2,ls,ms2) &
3528 - mode*jsolh_anal(5)*wws(ni,ls) * rnorms(1,ls,ms2))
3534 i = phi_mesh%jj(ni,m)
3535 src_phi(i,1) = src_phi(i,1)+rjs(ls,ms2)*ray*( &
3536 + jsolh_anal(3) *(dw_s(2,ni,ls,ms2) * rnorms(1,ls,ms2)&
3537 -dw_s(1,ni,ls,ms2) * rnorms(2,ls,ms2)))
3539 src_phi(i,2) = src_phi(i,2)+rjs(ls,ms2)*ray*( &
3540 + jsolh_anal(4)*(dw_s(2,ni,ls,ms2) * rnorms(1,ls,ms2)&
3541 -dw_s(1,ni,ls,ms2) * rnorms(2,ls,ms2)))
3548 i = interface_h_phi%jjs1(ni,ms)
3549 rhs_hl(:) = rhs_hl(:) + rhs_h(i,:)*w_cs(ni,ls)
3553 i = interface_h_phi%jjs2(ni,ms)
3554 src_phi(i,1) = src_phi(i,1)+rjs(ls,ms2)*ray*wws(ni,ls)*( &
3555 rhs_hl(1)*rnorms(1,ls,ms2) + rhs_hl(5)*rnorms(2,ls,ms2))
3556 src_phi(i,2) = src_phi(i,2)+rjs(ls,ms2)*ray*wws(ni,ls)*( &
3557 rhs_hl(2)*rnorms(1,ls,ms2) + rhs_hl(6)*rnorms(2,ls,ms2))
3565 IF (h_mesh%np /= 0)
THEN
3567 idxn_h = la_h%loc_to_glob(1,:)-1
3568 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,1), add_values, ierr)
3569 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,2), add_values, ierr)
3570 idxn_h = la_h%loc_to_glob(2,:)-1
3571 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,4), add_values, ierr)
3572 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,3), add_values, ierr)
3573 idxn_h = la_h%loc_to_glob(3,:)-1
3574 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,5), add_values, ierr)
3575 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,6), add_values, ierr)
3578 IF (phi_mesh%np /=0)
THEN
3580 idxn_phi = la_phi%loc_to_glob(1,:)-1
3581 CALL vecsetvalues(vb_1, phi_mesh%np, idxn_phi, src_phi(:,1), add_values, ierr)
3582 CALL vecsetvalues(vb_2, phi_mesh%np, idxn_phi, src_phi(:,2), add_values, ierr)
3586 CALL vecassemblybegin(vb_1,ierr)
3587 CALL vecassemblyend(vb_1,ierr)
3588 CALL vecassemblybegin(vb_2,ierr)
3589 CALL vecassemblyend(vb_2,ierr)
3598 SUBROUTINE surf_int(H_mesh,phi_mesh, interface_H_phi, interface_H_mu, &
3599 list_dirichlet_sides_h, sigma,mu_phi, mu_h_field,time,mode, &
3600 la_h, la_phi, vb_1, vb_2, r_fourier, index_fourier)
3607 TYPE(mesh_type),
INTENT(IN) :: h_mesh, phi_mesh
3608 TYPE(interface_type),
INTENT(IN) :: interface_h_phi, interface_h_mu
3609 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_dirichlet_sides_h
3610 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
3611 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN):: sigma
3612 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field
3613 INTEGER,
INTENT(IN) :: mode
3614 REAL(KIND=8),
OPTIONAL :: r_fourier
3615 INTEGER,
OPTIONAL :: index_fourier
3616 REAL(KIND=8),
DIMENSION(H_mesh%np, 6) :: src_h
3617 REAL(KIND=8),
DIMENSION(phi_mesh%np, 2) :: src_phi
3621 REAL(KIND=8) :: ray, muhl
3622 INTEGER :: ms, ls, ns, i, k, m, n, ni
3623 REAL(KIND=8),
DIMENSION(2) :: gaussp
3624 REAL(KIND=8),
DIMENSION(6) :: esolh_anal, esolphi_anal
3627 INTEGER,
DIMENSION(H_mesh%np) :: idxn_h
3628 INTEGER,
DIMENSION(phi_mesh%np) :: idxn_phi
3632 INTEGER,
DIMENSION(:),
ALLOCATABLE,
SAVE :: neumann_bdy_h_sides
3633 INTEGER,
DIMENSION(:),
ALLOCATABLE,
SAVE :: neumann_bdy_phi_sides
3634 LOGICAL,
SAVE :: once_neumann_bdy=.true.
3635 LOGICAL,
DIMENSION(:),
ALLOCATABLE :: virgin1
3636 LOGICAL,
DIMENSION(:),
ALLOCATABLE :: virgin2
3637 #include "petsc/finclude/petsc.h"
3638 petscerrorcode :: ierr
3645 IF (once_neumann_bdy)
THEN
3646 once_neumann_bdy = .false.
3647 ALLOCATE(virgin1(h_mesh%mes))
3648 ALLOCATE(virgin2(phi_mesh%mes))
3651 IF (interface_h_phi%mes/=0)
THEN
3652 virgin1(interface_h_phi%mesh1) = .false.
3653 virgin2(interface_h_phi%mesh2) = .false.
3655 IF (interface_h_mu%mes/=0)
THEN
3656 virgin1(interface_h_mu%mesh1) = .false.
3657 virgin1(interface_h_mu%mesh2) = .false.
3661 DO ms = 1, h_mesh%mes
3662 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))).LT.1d-12) cycle
3663 IF (.NOT.virgin1(ms)) cycle
3664 IF(minval(abs(h_mesh%sides(ms)-list_dirichlet_sides_h))==0) cycle
3667 ALLOCATE(neumann_bdy_h_sides(count))
3669 DO ms = 1, h_mesh%mes
3670 IF (maxval(abs(h_mesh%rr(1,h_mesh%jjs(:,ms)))).LT.1d-12) cycle
3671 IF (.NOT.virgin1(ms)) cycle
3672 IF(minval(abs(h_mesh%sides(ms)-list_dirichlet_sides_h))==0) cycle
3674 neumann_bdy_h_sides(count) = ms
3678 DO ms = 1, phi_mesh%mes
3679 IF (present(index_fourier))
THEN
3680 IF (phi_mesh%sides(ms)==index_fourier) cycle
3682 IF (.NOT.virgin2(ms)) cycle
3683 IF (maxval(abs(phi_mesh%rr(1,phi_mesh%jjs(:,ms)))).LT.1d-12) cycle
3684 IF (minval(abs(phi_mesh%sides(ms)-inputs%phi_list_dirichlet_sides))==0) cycle
3687 ALLOCATE(neumann_bdy_phi_sides(count))
3689 DO ms = 1, phi_mesh%mes
3690 IF (present(index_fourier))
THEN
3691 IF (phi_mesh%sides(ms)==index_fourier) cycle
3693 IF (.NOT.virgin2(ms)) cycle
3694 IF (maxval(abs(phi_mesh%rr(1,phi_mesh%jjs(:,ms)))).LT.1d-12) cycle
3695 IF (minval(abs(phi_mesh%sides(ms)-inputs%phi_list_dirichlet_sides))==0) cycle
3697 neumann_bdy_phi_sides(count) = ms
3706 DO count = 1,
SIZE(neumann_bdy_h_sides)
3707 ms = neumann_bdy_h_sides(count)
3708 m = h_mesh%neighs(ms)
3709 DO ls = 1, h_mesh%gauss%l_Gs
3711 muhl = sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
3716 DO ni = 1, h_mesh%gauss%n_ws; i = h_mesh%jjs(ni,ms)
3717 ray = ray + h_mesh%rr(1,i)* h_mesh%gauss%wws(ni,ls)
3720 IF (ray.LT.1.d-15) cycle
3723 DO ns=1, h_mesh%gauss%n_ws
3725 gaussp = gaussp + h_mesh%rr(:,i)*h_mesh%gauss%wws(ns,ls)
3729 esolh_anal(k) =
eexact_gauss(k,gaussp,mode,mu_phi,sigma(m),muhl, time)
3735 DO ns=1, h_mesh%gauss%n_ws
3736 i = h_mesh%jjs(ns,ms)
3737 src_h(i,1) = src_h(i,1)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3738 -esolh_anal(3)*h_mesh%gauss%wws(ns,ls)* &
3739 (h_mesh%gauss%rnorms(2,ls,ms)))
3741 src_h(i,2) = src_h(i,2)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3742 -esolh_anal(4)*h_mesh%gauss%wws(ns,ls)* &
3743 (h_mesh%gauss%rnorms(2,ls,ms)))
3745 src_h(i,3) = src_h(i,3)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3746 esolh_anal(1)*h_mesh%gauss%wws(ns,ls)* &
3747 (h_mesh%gauss%rnorms(2,ls,ms)) - &
3748 esolh_anal(5)*h_mesh%gauss%wws(ns,ls) * &
3749 (h_mesh%gauss%rnorms(1,ls,ms)))
3751 src_h(i,4) = src_h(i,4)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3752 esolh_anal(2)*h_mesh%gauss%wws(ns,ls) * &
3753 (h_mesh%gauss%rnorms(2,ls,ms)) - &
3754 esolh_anal(6)*h_mesh%gauss%wws(ns,ls) * &
3755 (h_mesh%gauss%rnorms(1,ls,ms)))
3757 src_h(i,5) = src_h(i,5)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3758 esolh_anal(3)*h_mesh%gauss%wws(ns,ls)* &
3759 (h_mesh%gauss%rnorms(1,ls,ms)))
3761 src_h(i,6) = src_h(i,6)-h_mesh%gauss%rjs(ls,ms)*ray*( &
3762 esolh_anal(4)*h_mesh%gauss%wws(ns,ls) * &
3763 (h_mesh%gauss%rnorms(1,ls,ms)))
3771 DO count = 1,
SIZE(neumann_bdy_phi_sides)
3772 ms = neumann_bdy_phi_sides(count)
3773 m = phi_mesh%neighs(ms)
3774 DO ls = 1, phi_mesh%gauss%l_Gs
3777 DO ni = 1, phi_mesh%gauss%n_ws; i = phi_mesh%jjs(ni,ms)
3778 ray = ray + phi_mesh%rr(1,i)* phi_mesh%gauss%wws(ni,ls)
3782 DO ns=1, phi_mesh%gauss%n_ws
3783 i=phi_mesh%jjs(ns,ms)
3784 gaussp = gaussp + phi_mesh%rr(:,i)*phi_mesh%gauss%wws(ns,ls)
3790 esolphi_anal(k) =
eexact_gauss(k,gaussp,mode,mu_phi,sigma(1),mu_h_field(1),time)
3794 DO ns=1, phi_mesh%gauss%n_ws
3795 i = phi_mesh%jjs(ns,ms)
3796 DO n = 1, phi_mesh%gauss%n_w
3797 IF (phi_mesh%jj(n,m) == i)
EXIT
3800 src_phi(i,1) = src_phi(i,1)-phi_mesh%gauss%rjs(ls,ms)*ray*( &
3801 +esolphi_anal(3)*(phi_mesh%gauss%dw_s(2,n,ls,ms)*phi_mesh%gauss%rnorms(1,ls,ms) &
3802 -phi_mesh%gauss%dw_s(1,n,ls,ms)*phi_mesh%gauss%rnorms(2,ls,ms))) &
3803 -phi_mesh%gauss%rjs(ls,ms)*(&
3804 -mode*esolphi_anal(2)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(2,ls,ms) &
3805 +mode*esolphi_anal(6)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(1,ls,ms))
3807 src_phi(i,2) = src_phi(i,2)-phi_mesh%gauss%rjs(ls,ms)*ray*( &
3808 +esolphi_anal(4)*(phi_mesh%gauss%dw_s(2,n,ls,ms)*phi_mesh%gauss%rnorms(1,ls,ms) &
3809 -phi_mesh%gauss%dw_s(1,n,ls,ms)*phi_mesh%gauss%rnorms(2,ls,ms))) &
3810 -phi_mesh%gauss%rjs(ls,ms)*(&
3811 mode*esolphi_anal(1)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(2,ls,ms) &
3812 -mode*esolphi_anal(5)*phi_mesh%gauss%wws(ns,ls)*phi_mesh%gauss%rnorms(1,ls,ms))
3829 IF (present(r_fourier))
THEN
3830 IF (r_fourier.GE.0.d0) CALL
error_petsc(
'maxwell_update_time_with_B: R_fourier should be -1')
3853 IF (h_mesh%mes/=0)
THEN
3855 idxn_h = la_h%loc_to_glob(1,:)-1
3856 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,1), add_values, ierr)
3857 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,2), add_values, ierr)
3858 idxn_h = la_h%loc_to_glob(2,:)-1
3859 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,4), add_values, ierr)
3860 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,3), add_values, ierr)
3861 idxn_h = la_h%loc_to_glob(3,:)-1
3862 CALL vecsetvalues(vb_1, h_mesh%np, idxn_h, src_h(:,5), add_values, ierr)
3863 CALL vecsetvalues(vb_2, h_mesh%np, idxn_h, src_h(:,6), add_values, ierr)
3866 IF (phi_mesh%mes/=0)
THEN
3868 idxn_phi = la_phi%loc_to_glob(1,:)-1
3869 CALL vecsetvalues(vb_1, phi_mesh%np, idxn_phi, src_phi(:,1), add_values, ierr)
3870 CALL vecsetvalues(vb_2, phi_mesh%np, idxn_phi, src_phi(:,2), add_values, ierr)
3874 CALL vecassemblybegin(vb_1,ierr)
3875 CALL vecassemblyend(vb_1,ierr)
3876 CALL vecassemblybegin(vb_2,ierr)
3877 CALL vecassemblyend(vb_2,ierr)
3884 mu_h_field, sigma, la_h, h_p_phi_mat1, h_p_phi_mat2)
3891 INTEGER,
INTENT(IN) :: mode
3892 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
3893 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma, mu_h_field
3894 INTEGER :: ms, ls, ni, nj, i, j, &
3895 n_ws1, n_ws2, n_w2, n_w1, m1, m2, ki, kj,ib,jb, ms1, ms2
3896 REAL(KIND=8) :: x, y,
z, norm, hm1
3897 REAL(KIND=8) :: ray, stab_colle_h_mu
3898 LOGICAL :: mark=.false.
3899 REAL(KIND=8),
DIMENSION(9,SIZE(H_mesh%jj,1),SIZE(H_mesh%jj,1),2,2) :: hsij, gsij
3914 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_ws,H_mesh%gauss%l_Gs) :: w_cs
3915 REAL(KIND=8),
DIMENSION(2, H_mesh%gauss%n_w, H_mesh%gauss%l_Gs, H_mesh%mes) :: dw_cs
3916 REAL(KIND=8),
DIMENSION(2, H_mesh%gauss%n_w) :: dwsi,dwsj
3917 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%l_Gs) :: gauss1, gauss2
3923 REAL(KIND=8),
DIMENSION(2) :: normi, normj
3924 REAL(KIND=8),
DIMENSION(SIZE(H_mesh%jjs,1)) :: wwsi, wwsj
3925 INTEGER :: n_wsi, n_wsj, ci, cj, n_wi, n_wj
3928 REAL(KIND=8) ::
ref, diff, mu_h, muhl1, muhl2, muhi, muhj, sigmai, sigmaj
3930 REAL(KIND=8) :: c_sym =.0d0
3931 REAL(KIND=8) :: wwiwwj, normt, stab_div
3933 REAL(KIND=8) :: drmuhl1, drmuhl2, drmuhj, dzmuhl1, dzmuhl2, dzmuhj
3939 REAL(KIND=8),
DIMENSION(6*H_mesh%gauss%n_w,6*H_mesh%gauss%n_w) :: mat_loc1, mat_loc2
3940 INTEGER ,
DIMENSION(6*H_mesh%gauss%n_w) :: idxn, jdxn
3944 #include "petsc/finclude/petsc.h"
3945 mat :: h_p_phi_mat1, h_p_phi_mat2
3946 petscerrorcode :: ierr
3949 stab_colle_h_mu = stab(3)
3959 n_ws1 = h_mesh%gauss%n_ws
3960 n_ws2 = h_mesh%gauss%n_ws
3961 n_w1 = h_mesh%gauss%n_w
3962 n_w2 = h_mesh%gauss%n_w
3974 DO ms = 1, interface_h_mu%mes
3976 ms2 = interface_h_mu%mesh2(ms)
3977 m2 = h_mesh%neighs(ms2)
3978 ms1 = interface_h_mu%mesh1(ms)
3979 m1 = h_mesh%neighs(ms1)
3981 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
3982 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(1,ms2)))**2)
3983 IF (diff/
ref .LT. 1.d-10)
THEN
3987 w_cs(1,ls)= wws(2,ls)
3988 w_cs(2,ls)= wws(1,ls)
3989 IF (n_ws1==3) w_cs(n_ws1,ls) = wws(n_ws1,ls)
3990 WRITE(*,*)
' Ouaps! oder of shape functions changed?'
3995 gauss2(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))*h_mesh%gauss%wws(:,ls))
3996 gauss2(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms2))*h_mesh%gauss%wws(:,ls))
3997 gauss1(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms1))*h_mesh%gauss%wws(:,ls))
3998 gauss1(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms1))*h_mesh%gauss%wws(:,ls))
4002 ref = sqrt(1.d-8+sum(gauss2(:,ls2)**2))
4005 diff = sqrt(sum((gauss2(:,ls2)-gauss1(:,ls1))**2))
4006 IF (diff .LT. 1.d-10)
THEN
4007 dw_cs(:,:,ls2,ms1) = h_mesh%gauss%dw_s(:,:,ls1,ms1)
4012 IF (.NOT.mark)
WRITE(*,*)
' BUG '
4017 DO ms = 1, interface_h_mu%mes
4019 ms2 = interface_h_mu%mesh2(ms)
4020 ms1 = interface_h_mu%mesh1(ms)
4021 m2 = h_mesh%neighs(ms2)
4022 m1 = h_mesh%neighs(ms1)
4023 mu_h = sum(mu_h_field(h_mesh%jj(:,m1)))/h_mesh%gauss%n_w
4026 hm1 = 1/sum(rjs(:,ms2))
4038 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
4039 x = hm1*rjs(ls,ms2)*ray
4042 muhl1 = sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4043 muhl2 = sum(mu_h_field(h_mesh%jjs(:,ms2))* wws(:,ls))
4045 normt =stab_colle_h_mu
4046 norm = stab_div*sum(rjs(:,ms2))**(2*alpha)
4056 normi = rnorms(:,ls,ms1)
4061 normi = rnorms(:,ls,ms2)
4068 normj = rnorms(:,ls,ms1)
4073 normj = rnorms(:,ls,ms2)
4081 wwiwwj = x * wwsi(ni)*wwsj(nj)
4087 y = normt * wwiwwj/muhj
4088 z = norm * muhi * wwiwwj
4090 hsij(1,ni,nj,ci,cj) = hsij(1,ni,nj,ci,cj) + y*normi(2)*normj(2) &
4091 +
z*normi(1)*normj(1)
4092 hsij(4,ni,nj,ci,cj) = hsij(4,ni,nj,ci,cj) - y*normj(1)*normi(2) &
4093 +
z*normi(1)*normj(2)
4095 hsij(7,ni,nj,ci,cj) = hsij(7,ni,nj,ci,cj) - y*normj(2)*normi(1) &
4096 +
z*normi(2)*normj(1)
4098 hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) + y*(normi(1)*normj(1) + normi(2)*normj(2))
4099 hsij(6,ni,nj,ci,cj) = hsij(6,ni,nj,ci,cj) + y*normi(1)*normj(1) &
4100 +
z*normi(2)*normj(2)
4113 i = interface_h_mu%jjs1(ni,ms)
4115 i = interface_h_mu%jjs2(ni,ms)
4117 ib = la_h%loc_to_glob(ki,i)
4118 ix = ni + n_ws1*((ki-1) + 3*(ci-1))
4125 j = interface_h_mu%jjs1(nj,ms)
4127 j = interface_h_mu%jjs2(nj,ms)
4129 jb = la_h%loc_to_glob(kj,j)
4130 jx = nj + n_ws2*((kj-1) + 3*(cj-1))
4132 IF ((ki == 1) .AND. (kj == 1))
THEN
4133 mat_loc1(ix,jx) = hsij(1,ni,nj,ci,cj)
4134 mat_loc2(ix,jx) = hsij(1,ni,nj,ci,cj)
4135 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN
4136 mat_loc1(ix,jx) = hsij(4,ni,nj,ci,cj)
4137 mat_loc2(ix,jx) = hsij(4,ni,nj,ci,cj)
4138 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN
4143 mat_loc1(ix,jx) = hsij(7,ni,nj,ci,cj)
4144 mat_loc2(ix,jx) = hsij(7,ni,nj,ci,cj)
4145 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
4146 mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)
4147 mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)
4148 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN
4149 mat_loc1(ix,jx) = hsij(6,ni,nj,ci,cj)
4150 mat_loc2(ix,jx) = hsij(6,ni,nj,ci,cj)
4159 CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
4160 mat_loc1(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
4161 CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
4162 mat_loc2(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
4171 DO ls = 1, h_mesh%gauss%l_Gs
4174 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
4177 muhl1 = sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4178 muhl2 = sum(mu_h_field(h_mesh%jjs(:,ms2))* wws(:,ls))
4179 drmuhl1 = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(1,:,ls,ms1))
4180 drmuhl2 = sum(mu_h_field(h_mesh%jj(:,m2))* dw_s(1,:,ls,ms2))
4181 dzmuhl1 = sum(mu_h_field(h_mesh%jj(:,m1))*dw_cs(2,:,ls,ms1))
4182 dzmuhl2 = sum(mu_h_field(h_mesh%jj(:,m2))* dw_s(2,:,ls,ms2))
4189 normi = rnorms(:,ls,ms1)
4197 normi = rnorms(:,ls,ms2)
4207 normj = rnorms(:,ls,ms1)
4220 normj = rnorms(:,ls,ms2)
4238 y = x*wwsi(ni)*wwsj(nj)/(2*sigmaj*muhj)
4240 hsij(2,ni,nj,ci,cj) = hsij(2,ni,nj,ci,cj) + y * (-mode/ray)*normi(1)
4241 hsij(3,ni,nj,ci,cj) = hsij(3,ni,nj,ci,cj) + y * mode/ray *normi(1)
4244 hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) - y * normi(1) *((+1/ray)+(-1/muhj)*drmuhj) &
4245 + y * normi(2)*(1/muhj)*dzmuhj
4247 hsij(8,ni,nj,ci,cj) = hsij(8,ni,nj,ci,cj) + y * (-mode/ray)*normi(2)
4248 hsij(9,ni,nj,ci,cj) = hsij(9,ni,nj,ci,cj) + y * mode/ray *normi(2)
4250 hsij(1,ni,nj,ci,cj) = hsij(1,ni,nj,ci,cj) - y * normi(2)*(-1/muhj)*dzmuhj
4251 hsij(4,ni,nj,ci,cj) = hsij(4,ni,nj,ci,cj) - y * normi(2)* (1/muhj)*drmuhj
4252 hsij(7,ni,nj,ci,cj) = hsij(7,ni,nj,ci,cj) + y * normi(1)*(-1/muhj)*dzmuhj
4253 hsij(6,ni,nj,ci,cj) = hsij(6,ni,nj,ci,cj) + y * normi(1)* (1/muhj)*drmuhj
4256 y = x*wwsi(ni)*wwsj(nj)/(2*sigmai)
4257 gsij(2,ni,nj,ci,cj) = gsij(2,ni,nj,ci,cj) + y * (-mode/ray)*normj(1)
4258 gsij(3,ni,nj,ci,cj) = gsij(3,ni,nj,ci,cj) + y * ( mode/ray)*normj(1)
4259 gsij(5,ni,nj,ci,cj) = gsij(5,ni,nj,ci,cj) + y * (-1/ray) *normj(1)
4260 gsij(8,ni,nj,ci,cj) = gsij(8,ni,nj,ci,cj) + y * (-mode/ray)*normj(2)
4261 gsij(9,ni,nj,ci,cj) = gsij(9,ni,nj,ci,cj) + y * mode/ray *normj(2)
4279 i = interface_h_mu%jjs1(ni,ms)
4281 i = interface_h_mu%jjs2(ni,ms)
4283 ib = la_h%loc_to_glob(ki,i)
4284 ix = ni + n_wsi*((ki-1) + 3*(ci-1))
4290 j = interface_h_mu%jjs1(nj,ms)
4292 j = interface_h_mu%jjs2(nj,ms)
4294 jb = la_h%loc_to_glob(kj,j)
4295 jx = nj + n_wsj*((kj-1) + 3*(cj-1))
4297 IF ((ki == 1) .AND. (kj == 1))
THEN
4298 mat_loc1(ix,jx) = hsij(1,ni,nj,ci,cj)
4299 mat_loc2(ix,jx) = hsij(1,ni,nj,ci,cj)
4300 ELSE IF((ki == 1) .AND. (kj == 2))
THEN
4301 mat_loc1(ix,jx) = gsij(2,ni,nj,ci,cj)
4302 mat_loc2(ix,jx) = gsij(3,ni,nj,ci,cj)
4303 ELSE IF ((ki == 1) .AND. (kj == 3))
THEN
4304 mat_loc1(ix,jx) = hsij(4,ni,nj,ci,cj)
4305 mat_loc2(ix,jx) = hsij(4,ni,nj,ci,cj)
4306 ELSE IF ((ki == 2) .AND. (kj == 1))
THEN
4307 mat_loc1(ix,jx) = hsij(2,ni,nj,ci,cj)
4308 mat_loc2(ix,jx) = hsij(3,ni,nj,ci,cj)
4309 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
4310 mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)+gsij(5,ni,nj,ci,cj)
4311 mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)+gsij(5,ni,nj,ci,cj)
4312 ELSEIF ((ki == 2) .AND. (kj == 3))
THEN
4313 mat_loc1(ix,jx) = hsij(8,ni,nj,ci,cj)
4314 mat_loc2(ix,jx) = hsij(9,ni,nj,ci,cj)
4315 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN
4316 mat_loc1(ix,jx) = hsij(7,ni,nj,ci,cj)
4317 mat_loc2(ix,jx) = hsij(7,ni,nj,ci,cj)
4318 ELSEIF ((ki == 3) .AND. (kj == 2))
THEN
4319 mat_loc1(ix,jx) = gsij(8,ni,nj,ci,cj)
4320 mat_loc2(ix,jx) = gsij(9,ni,nj,ci,cj)
4321 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN
4322 mat_loc1(ix,jx) = hsij(6,ni,nj,ci,cj)
4323 mat_loc2(ix,jx) = hsij(6,ni,nj,ci,cj)
4332 CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
4333 mat_loc1(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
4334 CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_ws2, jdxn(1:6*n_ws2), &
4335 mat_loc2(1:6*n_ws1,1:6*n_ws2), add_values, ierr)
4340 DO ls = 1, h_mesh%gauss%l_Gs
4343 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
4346 muhl1 = sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4347 muhl2 = sum(mu_h_field(h_mesh%jjs(:,ms2))* wws(:,ls))
4355 normi = rnorms(:,ls,ms1)
4357 dwsi = dw_cs(:,:,ls,ms1)
4365 normi = rnorms(:,ls,ms2)
4367 dwsi = dw_s(:,:,ls,ms2)
4377 normj = rnorms(:,ls,ms1)
4379 dwsj = dw_cs(:,:,ls,ms1)
4388 normj = rnorms(:,ls,ms2)
4390 dwsj = dw_s(:,:,ls,ms2)
4405 y = x*wwsi(ni)/(2*sigmaj*muhj)
4407 hsij(1,ni,nj,ci,cj) = hsij(1,ni,nj,ci,cj) + y*(-dwsj(2,nj))*normi(2)
4408 hsij(4,ni,nj,ci,cj) = hsij(4,ni,nj,ci,cj) + y* dwsj(1,nj) *normi(2)
4409 hsij(5,ni,nj,ci,cj) = hsij(5,ni,nj,ci,cj) + y*(-dwsj(2,nj) *normi(2)-dwsj(1,nj)*normi(1))
4410 hsij(6,ni,nj,ci,cj) = hsij(6,ni,nj,ci,cj) + y*(-dwsj(1,nj))*normi(1)
4411 hsij(7,ni,nj,ci,cj) = hsij(7,ni,nj,ci,cj) + y* dwsj(2,nj) *normi(1)
4416 y = x*wwsj(nj)/(2*sigmai)
4417 gsij(1,ni,nj,ci,cj) = gsij(1,ni,nj,ci,cj) + y*(-dwsi(2,ni))*normj(2)
4418 gsij(4,ni,nj,ci,cj) = gsij(4,ni,nj,ci,cj) + y* dwsi(2,ni) *normj(1)
4419 gsij(5,ni,nj,ci,cj) = gsij(5,ni,nj,ci,cj) + y*(-dwsi(2,ni) *normj(2)-dwsi(1,ni)*normj(1))
4420 gsij(6,ni,nj,ci,cj) = gsij(6,ni,nj,ci,cj) + y*(-dwsi(1,ni))*normj(1)
4421 gsij(7,ni,nj,ci,cj) = gsij(7,ni,nj,ci,cj) + y* dwsi(1,ni) *normj(2)
4439 i = interface_h_mu%jjs1(ni,ms)
4441 i = interface_h_mu%jjs2(ni,ms)
4443 ib = la_h%loc_to_glob(ki,i)
4444 ix = ni + n_wsi*((ki-1) + 3*(ci-1))
4450 j = h_mesh%jj(nj,m1)
4452 j = h_mesh%jj(nj,m2)
4454 jb = la_h%loc_to_glob(kj,j)
4455 jx = nj + n_wj*((kj-1) + 3*(cj-1))
4457 IF ((ki == 1) .AND. (kj == 1))
THEN
4458 mat_loc1(ix,jx) = hsij(1,ni,nj,ci,cj)
4459 mat_loc2(ix,jx) = hsij(1,ni,nj,ci,cj)
4460 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN
4461 mat_loc1(ix,jx) = hsij(4,ni,nj,ci,cj)
4462 mat_loc2(ix,jx) = hsij(4,ni,nj,ci,cj)
4463 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN
4464 mat_loc1(ix,jx) = hsij(7,ni,nj,ci,cj)
4465 mat_loc2(ix,jx) = hsij(7,ni,nj,ci,cj)
4466 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
4467 mat_loc1(ix,jx) = hsij(5,ni,nj,ci,cj)
4468 mat_loc2(ix,jx) = hsij(5,ni,nj,ci,cj)
4469 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN
4470 mat_loc1(ix,jx) = hsij(6,ni,nj,ci,cj)
4471 mat_loc2(ix,jx) = hsij(6,ni,nj,ci,cj)
4481 CALL matsetvalues(h_p_phi_mat1, 6*n_ws1, idxn(1:6*n_ws1), 6*n_w2, jdxn(1:6*n_w2), &
4482 mat_loc1(1:6*n_ws1,1:6*n_w2), add_values, ierr)
4483 CALL matsetvalues(h_p_phi_mat2, 6*n_ws1, idxn(1:6*n_ws1), 6*n_w2, jdxn(1:6*n_w2), &
4484 mat_loc2(1:6*n_ws1,1:6*n_w2), add_values, ierr)
4492 i = h_mesh%jj(ni,m1)
4494 i = h_mesh%jj(ni,m2)
4496 ib = la_h%loc_to_glob(ki,i)
4497 ix = ni + n_wi*((ki-1) + 3*(ci-1))
4503 j = interface_h_mu%jjs1(nj,ms)
4505 j = interface_h_mu%jjs2(nj,ms)
4507 jb = la_h%loc_to_glob(kj,j)
4508 jx = nj + n_wsj*((kj-1) + 3*(cj-1))
4510 IF ((ki == 1) .AND. (kj == 1))
THEN
4511 mat_loc1(ix,jx) = gsij(1,ni,nj,ci,cj)
4512 mat_loc2(ix,jx) = gsij(1,ni,nj,ci,cj)
4513 ELSEIF ((ki == 1) .AND. (kj == 3))
THEN
4514 mat_loc1(ix,jx) = gsij(4,ni,nj,ci,cj)
4515 mat_loc2(ix,jx) = gsij(4,ni,nj,ci,cj)
4516 ELSEIF ((ki == 3) .AND. (kj == 1))
THEN
4517 mat_loc1(ix,jx) = gsij(7,ni,nj,ci,cj)
4518 mat_loc2(ix,jx) = gsij(7,ni,nj,ci,cj)
4519 ELSEIF ((ki == 2) .AND. (kj == 2))
THEN
4520 mat_loc1(ix,jx) = gsij(5,ni,nj,ci,cj)
4521 mat_loc2(ix,jx) = gsij(5,ni,nj,ci,cj)
4522 ELSEIF ((ki == 3) .AND. (kj == 3))
THEN
4523 mat_loc1(ix,jx) = gsij(6,ni,nj,ci,cj)
4524 mat_loc2(ix,jx) = gsij(6,ni,nj,ci,cj)
4533 CALL matsetvalues(h_p_phi_mat1, 6*n_w1, idxn(1:6*n_w1), 6*n_ws2, jdxn(1:6*n_ws2), &
4534 mat_loc1(1:6*n_w1,1:6*n_ws2), add_values, ierr)
4535 CALL matsetvalues(h_p_phi_mat2, 6*n_w1, idxn(1:6*n_w1), 6*n_ws2, jdxn(1:6*n_ws2), &
4536 mat_loc2(1:6*n_w1,1:6*n_ws2), add_values, ierr)
4540 CALL matassemblybegin(h_p_phi_mat1,mat_final_assembly,ierr)
4541 CALL matassemblyend(h_p_phi_mat1,mat_final_assembly,ierr)
4542 CALL matassemblybegin(h_p_phi_mat2,mat_final_assembly,ierr)
4543 CALL matassemblyend(h_p_phi_mat2,mat_final_assembly,ierr)
4550 SUBROUTINE courant_mu(H_mesh,interface_H_mu,sigma,mu_H_field,time,mode,nl, &
4551 la_h, vb_1, vb_2, b_ext)
4560 REAL(KIND=8),
INTENT(IN) :: time
4561 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN) :: sigma
4562 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field
4563 INTEGER,
INTENT(IN) :: mode
4564 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: nl
4565 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: b_ext
4566 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: src_h
4567 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_ws,H_mesh%gauss%l_Gs) :: w_cs
4568 REAL(KIND=8),
DIMENSION(2) :: normi, gaussp1, gaussp2
4569 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_ws) :: wwsi
4570 REAL(KIND=8) :: x, ray
4571 INTEGER :: i, ni, ms, k, ls, n_ws1, n_ws2, ms1, ms2, n_w1, n_w2, m1, m2, ci, n_wsi
4572 INTEGER :: mesh_id1, mesh_id2
4573 REAL(KIND=8),
DIMENSION(6) :: jsolh_anal, test, b_ext_l
4574 REAL(KIND=8) :: muhl1, muhl2,
ref, diff
4581 INTEGER,
DIMENSION(H_mesh%np) :: idxn
4584 #include "petsc/finclude/petsc.h"
4585 petscerrorcode :: ierr
4595 n_ws1 = h_mesh%gauss%n_ws
4596 n_ws2 = h_mesh%gauss%n_ws
4597 n_w1 = h_mesh%gauss%n_w
4598 n_w2 = h_mesh%gauss%n_w
4600 DO ms = 1, interface_h_mu%mes
4601 ms1 = interface_h_mu%mesh1(ms)
4602 ms2 = interface_h_mu%mesh2(ms)
4603 m1 = h_mesh%neighs(ms1)
4604 m2 = h_mesh%neighs(ms2)
4606 ref = 1.d-8+sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(2,ms1)))**2)
4607 diff = sum((h_mesh%rr(:,h_mesh%jjs(1,ms1)) - h_mesh%rr(:,h_mesh%jjs(1,ms2)))**2)
4608 IF (diff/
ref .LT. 1.d-10)
THEN
4612 w_cs(1,ls)= wws(2,ls)
4613 w_cs(2,ls)= wws(1,ls)
4614 IF (n_ws1==3) w_cs(n_ws1,ls) = wws(n_ws1,ls)
4615 WRITE(*,*)
' Ouaps! oder of shape functions changed?'
4620 DO ms = 1, interface_h_mu%mes
4621 ms2 = interface_h_mu%mesh2(ms)
4622 ms1 = interface_h_mu%mesh1(ms)
4623 m2 = h_mesh%neighs(ms2)
4624 m1 = h_mesh%neighs(ms1)
4625 mesh_id1 = h_mesh%i_d(m1)
4626 mesh_id2 = h_mesh%i_d(m2)
4629 ray = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* h_mesh%gauss%wws(:,ls))
4634 b_ext_l(k) = sum(b_ext(h_mesh%jjs(:,ms1),k)*h_mesh%gauss%wws(:,ls))
4636 muhl1=sum(mu_h_field(h_mesh%jjs(:,ms1))*w_cs(:,ls))
4637 gaussp1(1) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms1))*w_cs(:,ls))
4638 gaussp1(2) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms1))*w_cs(:,ls))
4640 jsolh_anal(k) =
jexact_gauss(k, gaussp1, mode, one ,sigma(m1), &
4641 muhl1, time, mesh_id1, b_ext_l)/sigma(m1) &
4642 + sum(nl(h_mesh%jjs(1:n_ws1,ms1),k)*w_cs(1:n_ws1,ls))
4650 b_ext_l(k) = sum(b_ext(h_mesh%jjs(:,ms2),k)*h_mesh%gauss%wws(:,ls))
4652 muhl2=sum(mu_h_field(h_mesh%jjs(:,ms2))*wws(:,ls))
4653 gaussp2(1) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))*wws(:,ls))
4654 gaussp2(2) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms2))*wws(:,ls))
4655 IF (maxval(abs(gaussp1-gaussp2)) > 1.d-11)
THEN
4656 WRITE(*,*)
' BUG courant_mu '
4660 test(k) =
jexact_gauss(k, gaussp2, mode, one ,sigma(m2), &
4661 muhl2, time, mesh_id2, b_ext_l)/sigma(m2) &
4662 + sum(nl(h_mesh%jjs(1:n_ws2,ms2),k)*wws(1:n_ws2,ls))
4663 jsolh_anal(k) = jsolh_anal(k) + test(k)
4673 normi = rnorms(:,ls,ms1)
4677 normi = rnorms(:,ls,ms2)
4683 i = interface_h_mu%jjs1(ni,ms)
4685 i = interface_h_mu%jjs2(ni,ms)
4687 x = rjs(ls,ms2)*ray*wwsi(ni)/2
4688 src_h(i,1) = src_h(i,1)+x*(-jsolh_anal(3)*normi(2))
4689 src_h(i,2) = src_h(i,2)+x*(-jsolh_anal(4)*normi(2))
4690 src_h(i,3) = src_h(i,3)+x*(jsolh_anal(1)*normi(2)-jsolh_anal(5)*normi(1))
4691 src_h(i,4) = src_h(i,4)+x*(jsolh_anal(2)*normi(2)-jsolh_anal(6)*normi(1))
4692 src_h(i,5) = src_h(i,5)+x*(jsolh_anal(3)*normi(1))
4693 src_h(i,6) = src_h(i,6)+x*(jsolh_anal(4)*normi(1))
4699 IF (h_mesh%np /= 0)
THEN
4701 idxn = la_h%loc_to_glob(1,:)-1
4702 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,1), add_values, ierr)
4703 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,2), add_values, ierr)
4704 idxn = la_h%loc_to_glob(2,:)-1
4705 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,4), add_values, ierr)
4706 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,3), add_values, ierr)
4707 idxn = la_h%loc_to_glob(3,:)-1
4708 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,5), add_values, ierr)
4709 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,6), add_values, ierr)
4712 CALL vecassemblybegin(vb_1,ierr)
4713 CALL vecassemblyend(vb_1,ierr)
4714 CALL vecassemblybegin(vb_2,ierr)
4715 CALL vecassemblyend(vb_2,ierr)
4719 SUBROUTINE rhs_dirichlet(H_mesh,Dirichlet_bdy_H_sides,sigma,&
4720 mu_h_field,time,mode,nl,stab, la_h, vb_1, vb_2, b_ext, j_over_sigma, &
4729 INTEGER,
DIMENSION(:),
INTENT(IN) :: dirichlet_bdy_h_sides
4730 REAL(KIND=8),
INTENT(IN) :: time
4731 REAL(KIND=8),
DIMENSION(H_mesh%me),
INTENT(IN) :: sigma
4732 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field
4733 INTEGER,
INTENT(IN) :: mode
4734 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: nl
4735 REAL(KIND=8),
DIMENSION(3),
INTENT(IN) :: stab
4736 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: b_ext
4737 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: j_over_sigma
4738 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: sigma_curl_bdy_in
4739 REAL(KIND=8),
DIMENSION(H_mesh%np,6) :: src_h
4740 REAL(KIND=8),
DIMENSION(2) :: gaussp1
4741 REAL(KIND=8) :: x, ray, stab_colle_h_mu
4742 INTEGER :: i, ni, ms, k, ls, m1, count
4744 REAL(KIND=8),
DIMENSION(6) :: jsolh_anal, b_ext_l
4745 REAL(KIND=8) :: muhl1, hm1
4746 REAL(KIND=8),
DIMENSION(6,H_mesh%gauss%l_Gs) :: hloc, hlocxn
4747 REAL(KIND=8),
DIMENSION(2,H_mesh%gauss%l_Gs) :: rloc
4748 REAL(KIND=8),
DIMENSION(1) :: muloc
4756 INTEGER,
DIMENSION(H_mesh%np) :: idxn
4759 #include "petsc/finclude/petsc.h"
4760 petscerrorcode :: ierr
4774 stab_colle_h_mu = stab(3)
4777 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
4778 ms = dirichlet_bdy_h_sides(count)
4779 hm1 = stab_colle_h_mu/sum(h_mesh%gauss%rjs(:,ms))
4780 m1 = h_mesh%neighs(ms)
4781 mesh_id1 = h_mesh%i_d(m1)
4782 muloc(1) = mu_h_field(h_mesh%jj(1,m1))
4783 DO ls = 1, h_mesh%gauss%l_Gs
4784 rloc(1,ls) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
4785 rloc(2,ls) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
4789 hloc(k,:) =
hexact(h_mesh, k, rloc, mode, muloc, time)
4792 hlocxn(1,:) = hloc(3,:)*h_mesh%gauss%rnorms(2,:,ms)
4793 hlocxn(2,:) = hloc(4,:)*h_mesh%gauss%rnorms(2,:,ms)
4794 hlocxn(3,:) = hloc(5,:)*h_mesh%gauss%rnorms(1,:,ms)-hloc(1,:)*h_mesh%gauss%rnorms(2,:,ms)
4795 hlocxn(4,:) = hloc(6,:)*h_mesh%gauss%rnorms(1,:,ms)-hloc(2,:)*h_mesh%gauss%rnorms(2,:,ms)
4796 hlocxn(5,:) = -hloc(3,:)*h_mesh%gauss%rnorms(1,:,ms)
4797 hlocxn(6,:) = -hloc(4,:)*h_mesh%gauss%rnorms(1,:,ms)
4799 DO ls = 1, h_mesh%gauss%l_Gs
4804 b_ext_l(k) = sum(b_ext(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls))
4808 muhl1=sum(mu_h_field(h_mesh%jjs(:,ms))*h_mesh%gauss%wws(:,ls))
4809 gaussp1(1) = rloc(1,ls)
4810 gaussp1(2) = rloc(2,ls)
4818 IF (inputs%if_level_set.AND.inputs%variation_sigma_fluid)
THEN
4820 jsolh_anal(k) = sum(j_over_sigma(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls)) &
4821 + sigma_curl_bdy_in(index,k) &
4822 + sum(nl(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls)) &
4827 jsolh_anal(k) =
jexact_gauss(k, gaussp1, mode, one ,sigma(m1), muhl1, &
4828 time, mesh_id1, b_ext_l)/sigma(m1) &
4829 + sum(nl(h_mesh%jjs(:,ms),k)*h_mesh%gauss%wws(:,ls)) &
4834 DO ni = 1, h_mesh%gauss%n_ws
4835 i = h_mesh%jjs(ni,ms)
4836 x = h_mesh%gauss%rjs(ls,ms)*ray*h_mesh%gauss%wws(ni,ls)
4838 src_h(i,1) = src_h(i,1)+x*(-jsolh_anal(3)*h_mesh%gauss%rnorms(2,ls,ms))
4839 src_h(i,2) = src_h(i,2)+x*(-jsolh_anal(4)*h_mesh%gauss%rnorms(2,ls,ms))
4840 src_h(i,3) = src_h(i,3)+x*(jsolh_anal(1)*h_mesh%gauss%rnorms(2,ls,ms)&
4841 -jsolh_anal(5)*h_mesh%gauss%rnorms(1,ls,ms))
4842 src_h(i,4) = src_h(i,4)+x*(jsolh_anal(2)*h_mesh%gauss%rnorms(2,ls,ms)&
4843 -jsolh_anal(6)*h_mesh%gauss%rnorms(1,ls,ms))
4844 src_h(i,5) = src_h(i,5)+x*(jsolh_anal(3)*h_mesh%gauss%rnorms(1,ls,ms))
4845 src_h(i,6) = src_h(i,6)+x*(jsolh_anal(4)*h_mesh%gauss%rnorms(1,ls,ms))
4852 IF (h_mesh%np /= 0)
THEN
4854 idxn = la_h%loc_to_glob(1,:)-1
4855 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,1), add_values, ierr)
4856 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,2), add_values, ierr)
4857 idxn = la_h%loc_to_glob(2,:)-1
4858 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,4), add_values, ierr)
4859 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,3), add_values, ierr)
4860 idxn = la_h%loc_to_glob(3,:)-1
4861 CALL vecsetvalues(vb_1, h_mesh%np, idxn, src_h(:,5), add_values, ierr)
4862 CALL vecsetvalues(vb_2, h_mesh%np, idxn, src_h(:,6), add_values, ierr)
4865 CALL vecassemblybegin(vb_1,ierr)
4866 CALL vecassemblyend(vb_1,ierr)
4867 CALL vecassemblybegin(vb_2,ierr)
4868 CALL vecassemblyend(vb_2,ierr)
4880 INTEGER,
POINTER,
DIMENSION(:) :: js_d
4881 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: on_proc_loc, on_proc, not_cav_loc, not_cav
4882 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: is_ok, j_tmp
4883 INTEGER,
DIMENSION(1) :: loc
4884 INTEGER :: m, ms, i, nb_dom, idx, nb_cav, ni
4886 #include "petsc/finclude/petsc.h"
4887 mpi_comm,
INTENT(IN) :: communicator
4889 petscerrorcode :: ierr
4891 IF (inputs%nb_dom_phi==0)
RETURN
4893 CALL mpi_comm_rank(communicator, rank, ierr)
4895 nb_dom = inputs%nb_dom_phi
4896 ALLOCATE(on_proc_loc(nb_dom), on_proc(nb_dom))
4897 ALLOCATE(not_cav_loc(nb_dom), not_cav(nb_dom))
4905 IF (minval(abs(inputs%list_dom_phi-i)) /= 0)
THEN
4906 WRITE(*,*)
'error in dirichlet cavities'
4908 loc = minloc(abs(inputs%list_dom_phi-i))
4909 on_proc_loc(loc(1)) = rank
4911 IF (mesh%mes /= 0)
THEN
4912 ALLOCATE(is_ok(mesh%mes))
4913 is_ok = mesh%i_d(mesh%neighs)
4914 IF (interface_h_phi%mes /=0)
THEN
4915 is_ok(interface_h_phi%mesh2) = 0
4918 IF (sum(abs(mesh%rr(1,mesh%jjs(:,ms)))) .LT. 1.d-12)
THEN
4921 IF (inputs%my_periodic%nb_periodic_pairs /=0)
THEN
4922 IF (minval(abs(inputs%my_periodic%list_periodic-mesh%sides(ms))) == 0)
THEN
4929 IF (is_ok(ms) == 0) cycle
4931 IF (minval(abs(inputs%list_dom_phi-i)) /= 0)
THEN
4932 WRITE(*,*)
'error in dirichlet cavities'
4934 loc = minloc(abs(inputs%list_dom_phi-i))
4935 not_cav_loc(loc(1)) = rank
4938 CALL mpi_allreduce(on_proc_loc, on_proc, nb_dom, mpi_integer, mpi_max, communicator, ierr)
4939 CALL mpi_allreduce(not_cav_loc, not_cav, nb_dom, mpi_integer, mpi_max, communicator, ierr)
4941 ALLOCATE(j_tmp(
SIZE(js_d)+nb_dom))
4942 j_tmp(1:
SIZE(js_d)) = js_d
4945 IF ( (not_cav(i)==-1) .AND. (on_proc(i)==rank) )
THEN
4949 IF (mesh%i_d(m) == inputs%list_dom_phi(i))
THEN
4950 DO ni = 1, mesh%gauss%n_w
4951 IF (minval(abs(mesh%jjs-mesh%jj(ni,m))) /=0)
THEN
4952 j_tmp(idx) = mesh%jj(ni,m)
4958 WRITE(*,*)
'add ', j_tmp(idx),
'in dom ', inputs%list_dom_phi(i),
' : proc ', rank
4959 WRITE(*,*)
'add ', mesh%rr(:,j_tmp(idx)), mesh%i_d(m)
4967 nb_cav = idx -
SIZE(js_d)
4968 IF (nb_cav /= 0)
THEN
4974 DEALLOCATE(on_proc_loc, on_proc, j_tmp)
4975 DEALLOCATE(not_cav_loc, not_cav)
4976 IF (
ALLOCATED(is_ok))
DEALLOCATE(is_ok)
4978 WRITE(*,
'(a,x,i2,x,a,x,i2)')
'I have detected', nb_cav,
' cavity(ies) on proc', rank
4988 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: angles
4989 INTEGER,
INTENT(IN) :: nb_angles
4990 INTEGER,
INTENT(IN) :: nb, ne
4991 REAL(KIND=8),
INTENT(IN) :: time
4992 REAL(KIND=8),
DIMENSION(nb_angles,ne-nb+1) :: vv
4999 FUNCTION one_over_mu(H_mesh,angles,nb_angles,nb,ne,time) RESULT(vv)
5004 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: angles
5005 INTEGER,
INTENT(IN) :: nb_angles
5006 INTEGER,
INTENT(IN) :: nb, ne
5007 REAL(KIND=8),
INTENT(IN) :: time
5008 REAL(KIND=8),
DIMENSION(nb_angles,ne-nb+1) :: vv
5014 SUBROUTINE smb_sigma_prod_curl(communicator, mesh, list_mode, H_in, sigma_bar, sigma_in, V_out)
5023 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5024 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: h_in
5025 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_bar
5026 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: sigma_in
5027 REAL(KIND=8),
DIMENSION(:,:,:) :: v_out
5028 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: h_gauss, roth
5029 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: sigma_gauss
5030 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: roth_bar
5031 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
5032 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
5033 INTEGER :: m, l , i, mode, index, k
5034 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,6) :: h_in_loc
5035 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: sigma_in_loc
5037 INTEGER :: nb_procs, bloc_size, m_max_pad, code
5038 #include "petsc/finclude/petsc.h"
5039 mpi_comm :: communicator
5041 DO i = 1,
SIZE(list_mode)
5045 j_loc = mesh%jj(:,m)
5047 h_in_loc(:,k) = h_in(j_loc,k,i)
5050 sigma_in_loc(:,k) = sigma_in(j_loc,k,i)
5053 DO l = 1, mesh%gauss%l_G
5055 dw_loc = mesh%gauss%dw(:,:,l,m)
5058 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
5061 h_gauss(index,1,i) = sum(h_in_loc(:,1)*mesh%gauss%ww(:,l))
5062 h_gauss(index,3,i) = sum(h_in_loc(:,3)*mesh%gauss%ww(:,l))
5063 h_gauss(index,5,i) = sum(h_in_loc(:,5)*mesh%gauss%ww(:,l))
5065 h_gauss(index,2,i) = sum(h_in_loc(:,2)*mesh%gauss%ww(:,l))
5066 h_gauss(index,4,i) = sum(h_in_loc(:,4)*mesh%gauss%ww(:,l))
5067 h_gauss(index,6,i) = sum(h_in_loc(:,6)*mesh%gauss%ww(:,l))
5069 sigma_gauss(index,1,i) = sum(sigma_in_loc(:,1)*mesh%gauss%ww(:,l))
5070 sigma_gauss(index,2,i) = sum(sigma_in_loc(:,2)*mesh%gauss%ww(:,l))
5073 roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
5074 -sum(h_in_loc(:,3)*dw_loc(2,:))
5075 roth(index,4,i) = sum(h_in_loc(:,2)*dw_loc(2,:)) &
5076 -sum(h_in_loc(:,6)*dw_loc(1,:))
5077 roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
5078 +sum(h_in_loc(:,3)*dw_loc(1,:)) &
5079 -mode/ray*h_gauss(index,2,i)
5081 roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
5082 -sum(h_in_loc(:,4)*dw_loc(2,:))
5083 roth(index,3,i) = sum(h_in_loc(:,1)*dw_loc(2,:)) &
5084 -sum(h_in_loc(:,5)*dw_loc(1,:))
5085 roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
5086 +sum(h_in_loc(:,4)*dw_loc(1,:))&
5087 +mode/ray*h_gauss(index,1,i)
5090 roth_bar(index,k,i) = roth(index,k,i)/sum(sigma_bar(j_loc)*mesh%gauss%ww(:,l))
5096 CALL mpi_comm_size(communicator, nb_procs, code)
5097 bloc_size =
SIZE(sigma_gauss,1)/nb_procs+1
5098 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5101 bloc_size, m_max_pad)
5103 v_out = roth_bar - v_out
5107 SUBROUTINE smb_sigma_prod_curl_bdy(communicator, mesh, Dirichlet_bdy_H_sides, list_mode, H_in, sigma_bar, sigma_in, V_out)
5116 INTEGER,
DIMENSION(:),
INTENT(IN) :: dirichlet_bdy_h_sides
5117 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5118 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: h_in
5119 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: sigma_bar
5120 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: sigma_in
5121 REAL(KIND=8),
DIMENSION(:,:,:) :: v_out
5122 REAL(KIND=8),
DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),6,SIZE(list_mode)) :: h_gauss, roth
5123 REAL(KIND=8),
DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),6,SIZE(list_mode)) :: roth_bar
5124 REAL(KIND=8),
DIMENSION(mesh%gauss%l_Gs*SIZE(Dirichlet_bdy_H_sides),2,SIZE(list_mode)) :: sigma_gauss
5125 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
5126 INTEGER :: ms, ls , i, mode, index, k
5127 INTEGER,
DIMENSION(mesh%gauss%n_ws) :: j_loc
5128 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,6) :: h_in_loc
5129 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,2) :: sigma_in_loc
5131 INTEGER :: nb_procs, bloc_size, m_max_pad, code, count, m1
5132 #include "petsc/finclude/petsc.h"
5133 mpi_comm :: communicator
5135 DO i = 1,
SIZE(list_mode)
5138 DO count = 1,
SIZE(dirichlet_bdy_h_sides)
5139 ms = dirichlet_bdy_h_sides(count)
5140 m1 = mesh%neighs(ms)
5142 j_loc = mesh%jjs(:,ms)
5144 h_in_loc(:,k) = h_in(j_loc,k,i)
5147 sigma_in_loc(:,k) = sigma_in(j_loc,k,i)
5150 DO ls = 1, mesh%gauss%l_Gs
5152 dw_loc = mesh%gauss%dw_s(:,:,ls,ms)
5155 ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
5158 h_gauss(index,1,i) = sum(h_in_loc(:,1)*mesh%gauss%wws(:,ls))
5159 h_gauss(index,3,i) = sum(h_in_loc(:,3)*mesh%gauss%wws(:,ls))
5160 h_gauss(index,5,i) = sum(h_in_loc(:,5)*mesh%gauss%wws(:,ls))
5162 h_gauss(index,2,i) = sum(h_in_loc(:,2)*mesh%gauss%wws(:,ls))
5163 h_gauss(index,4,i) = sum(h_in_loc(:,4)*mesh%gauss%wws(:,ls))
5164 h_gauss(index,6,i) = sum(h_in_loc(:,6)*mesh%gauss%wws(:,ls))
5166 sigma_gauss(index,1,i) = sum(sigma_in_loc(:,1)*mesh%gauss%wws(:,ls))
5167 sigma_gauss(index,2,i) = sum(sigma_in_loc(:,2)*mesh%gauss%wws(:,ls))
5170 roth(index,1,i) = mode/ray*h_gauss(index,6,i) &
5172 -sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(2,:))
5176 roth(index,4,i) = sum(h_in(mesh%jj(:,m1),2,i)*dw_loc(2,:)) &
5177 -sum(h_in(mesh%jj(:,m1),6,i)*dw_loc(1,:))
5179 roth(index,5,i) = 1/ray*h_gauss(index,3,i) &
5181 +sum(h_in(mesh%jj(:,m1),3,i)*dw_loc(1,:)) &
5182 -mode/ray*h_gauss(index,2,i)
5185 roth(index,2,i) =-mode/ray*h_gauss(index,5,i) &
5187 -sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(2,:))
5191 roth(index,3,i) = sum(h_in(mesh%jj(:,m1),1,i)*dw_loc(2,:)) &
5192 -sum(h_in(mesh%jj(:,m1),5,i)*dw_loc(1,:))
5194 roth(index,6,i) = 1/ray*h_gauss(index,4,i) &
5196 +sum(h_in(mesh%jj(:,m1),4,i)*dw_loc(1,:))&
5197 +mode/ray*h_gauss(index,1,i)
5200 roth_bar(index,k,i) = roth(index,k,i)/(sum(sigma_bar(j_loc)*mesh%gauss%wws(:,ls)))
5206 IF (
SIZE(dirichlet_bdy_h_sides).GE.1)
THEN
5207 CALL mpi_comm_size(communicator, nb_procs, code)
5208 bloc_size =
SIZE(sigma_gauss,1)/nb_procs+1
5209 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5212 bloc_size, m_max_pad)
5214 v_out = roth_bar - v_out
5220 mu_h_field, mu_phi, sigma_tot, time, j_over_sigma)
5228 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
5229 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: sigma_tot
5230 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field
5231 REAL(KIND=8),
INTENT(IN) :: mu_phi, time
5232 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: j_over_sigma
5233 REAL(KIND=8),
DIMENSION(SIZE(mesh%rr,2),6,SIZE(list_mode)) :: j_node
5234 REAL(KIND=8),
DIMENSION(6) :: b_ext_l
5235 REAL(KIND=8) :: muhl
5236 INTEGER :: mode, mesh_id1, k, i, n
5237 INTEGER :: nb_procs, bloc_size, m_max_pad, code
5238 #include "petsc/finclude/petsc.h"
5239 mpi_comm :: communicator
5244 muhl=minval(mu_h_field)
5248 DO i = 1,
SIZE(list_mode)
5251 DO n = 1,
SIZE(mesh%rr,2)
5252 j_node(n,k,i) =
jexact_gauss(k, mesh%rr(:,n), mode, mu_phi, 1.d0, &
5253 muhl, time, mesh_id1, b_ext_l)
5258 CALL mpi_comm_size(communicator, nb_procs, code)
5259 bloc_size =
SIZE(j_node,1)/nb_procs+1
5260 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
5262 nb_procs, bloc_size, m_max_pad)
5267
subroutine surf_int(H_mesh, phi_mesh, interface_H_phi, interface_H_mu, list_dirichlet_sides_H, sigma, mu_phi, mu_H_field, time, mode, LA_H, LA_phi, vb_1, vb_2, R_fourier, index_fourier)
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, J_over_sigma)
subroutine, public create_my_ghost(mesh, LA, ifrom)
subroutine dirichlet_rhs(js_D, bs_D, b)
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, sigma_ns_in, jj_v_to_H)
real(kind=8) function, dimension(size(h_mesh%rr, 2)), public sigma_bar_in_fourier_space(H_mesh)
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
subroutine smb_sigma_prod_curl_bdy(communicator, mesh, Dirichlet_bdy_H_sides, list_mode, H_in, sigma_bar, sigma_in, V_out)
subroutine mat_dirichlet_maxwell(H_mesh, Dirichlet_bdy_H_sides, mode, stab, mu_H_field, LA_H, H_p_phi_mat1, H_p_phi_mat2, sigma_np)
subroutine dirichlet_nodes_parallel(mesh, list_dirichlet_sides, js_d)
subroutine mat_maxwell_mu(H_mesh, interface_H_mu, mode, stab, mu_H_field, sigma, LA_H, H_p_phi_mat1, H_p_phi_mat2)
real(kind=8) function, dimension(size(rr, 2), 6), public h_b_quasi_static(char_h_b, rr, m)
real(kind=8) function, public jexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t, mesh_id, opt_B_ext)
subroutine solver(my_ksp, b, x, reinit, verbose)
subroutine dirichlet_m_parallel(matrix, glob_js_D)
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(size(rr, 2)), public phiexact(TYPE, rr, m, mu_phi, t)
real(kind=8) function, dimension(2), public grad_mu_bar_in_fourier_space(pt, pt_id)
real(kind=8) function, dimension(nb_angles, ne-nb+1) one_over_mu(H_mesh, angles, nb_angles, nb, ne, time)
real(kind=8) function, dimension(nb_angles, ne-nb+1), public mu_in_real_space(H_mesh, angles, nb_angles, nb, ne, time)
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_np)
subroutine, public ref(communicator, V1_in, V2_in, V_out, temps)
subroutine, public fft_scalar_vect_dcl(communicator, V1_in, V2_in, V_out, pb, nb_procs, bloc_size, m_max_pad, temps)
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
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)
real(kind=8) function, public eexact_gauss(TYPE, rr, m, mu_phi, sigma, mu_H, t)
real(kind=8) function user_time()
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)
subroutine courant_mu(H_mesh, interface_H_mu, sigma, mu_H_field, time, mode, nl, LA_H, vb_1, vb_2, B_ext)
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, sigma_curl_bdy_in)
subroutine dirichlet_cavities(communicator, interface_H_phi, mesh, js_D)
real(kind=8) function, dimension(ne-nb+1), public mu_bar_in_fourier_space(H_mesh, nb, ne, pts, pts_ids)
subroutine error_petsc(string)
subroutine, public fft_par_cross_prod_dcl(communicator, V1_in, V2_in, V_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
subroutine scalar_with_bc_glob_js_d(pp_mesh, list_mode, pp_1_LA, pp_js_D, pp_mode_global_js_D)
subroutine, public extract(xghost, ks, ke, LA, phi)
subroutine vector_without_bc_glob_js_d(vv_mesh, list_mode, vv_3_LA, vv_mode_global_js_D)
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
real(kind=8) function, dimension(size(rr, 2)), public hexact(H_mesh, TYPE, rr, m, mu_H_field, t)
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic z
subroutine smb_current_over_sigma(communicator, mesh, list_mode, mu_H_field, mu_phi, sigma_tot, time, J_over_sigma)
subroutine smb_sigma_prod_curl(communicator, mesh, list_mode, H_in, sigma_bar, sigma_in, V_out)