diff --git a/unstructured/M3Dmodules.f90 b/unstructured/M3Dmodules.f90 index 3873c4b6a..270fa39f4 100644 --- a/unstructured/M3Dmodules.f90 +++ b/unstructured/M3Dmodules.f90 @@ -489,6 +489,7 @@ module arrays type(field_type) :: temporary_field type(field_type) :: psi_coil_field + type(field_type), allocatable :: psi_coil_fields(:) ! the indicies of the named fields within the field vector integer, parameter :: u_g = 1 diff --git a/unstructured/coils.f90 b/unstructured/coils.f90 index 89d1cfb22..227080626 100644 --- a/unstructured/coils.f90 +++ b/unstructured/coils.f90 @@ -13,7 +13,7 @@ module coils ! reads coil and current data from file !====================================================== subroutine load_coils(xc, zc, ic, numcoils, coil_filename, current_filename, & - coil_mask, filaments) + coil_mask, filaments, numcoils_base) use math use read_ascii @@ -27,6 +27,7 @@ subroutine load_coils(xc, zc, ic, numcoils, coil_filename, current_filename, & character*(*) :: coil_filename, current_filename ! input files integer, intent(out), dimension(maxfilaments), optional :: coil_mask integer, intent(out), dimension(maxfilaments), optional :: filaments + integer, intent(out), optional :: numcoils_base integer, parameter :: fcoil = 34 @@ -88,6 +89,7 @@ subroutine load_coils(xc, zc, ic, numcoils, coil_filename, current_filename, & deallocate(turnsv) end if + if(rank.eq.0) print *, coil_filename, ' ', current_filename call read_ascii_column(current_filename, cv, ncur, 0, 1) call read_ascii_column(current_filename, phasev, ncur, 0, 2) if(ncur.ne.ncoil) then @@ -104,6 +106,8 @@ subroutine load_coils(xc, zc, ic, numcoils, coil_filename, current_filename, & if(syv(i).lt.1) syv(i) = 1 end do + if(present(numcoils_base)) numcoils_base = ncoil + cv = amu0 * 1000. * cv / (sxv*syv) / twopi ! divide coils into sub-coils @@ -157,7 +161,7 @@ subroutine load_coils(xc, zc, ic, numcoils, coil_filename, current_filename, & end subroutine load_coils - subroutine field_from_coils(xc, zc, ic, nc, f, ipole, ierr) + subroutine field_from_coils(xc, zc, ic, nc, f, ipole, ierr, ncs, fs, mask) use basic use field use mesh_mod @@ -168,26 +172,122 @@ subroutine field_from_coils(xc, zc, ic, nc, f, ipole, ierr) include 'mpif.h' - real, intent(in), dimension(nc) :: xc, zc ! array of coil positions - complex, intent(in), dimension(nc) :: ic ! array of coil currents - integer, intent(in) :: nc ! number of coils - type(field_type), intent(inout) :: f ! poloidal flux field - integer, intent(in) :: ipole ! type of field to add + real, intent(in), dimension(nc) :: xc, zc ! array of coil positions + complex, intent(in), dimension(nc) :: ic ! array of coil currents + integer, intent(in) :: nc ! number of coils + type(field_type), intent(inout) :: f ! poloidal flux field + integer, intent(in) :: ipole ! type of field to add integer, intent(out) :: ierr - + ! poloidal flux field for each coil + integer, intent(in), optional :: ncs ! number of real coils + type(field_type), intent(in), dimension(ncs), optional :: fs(:) ! fields + integer, intent(in), dimension(nc), optional :: mask !filaments to coils + integer :: numelms, k, itri type(field_type) :: psi_vec vectype, dimension(dofs_per_element) :: dofs real, allocatable :: g(:,:) + real :: cur - if(nc.le.0) return + integer :: inode, numnodes, icounter_t + vectype, dimension(dofs_per_node) :: datain, dataout - allocate(g(MAX_PTS,nc)) + if(nc.le.0) return ierr = 0 + if(present(ncs)) then + if(.not.present(fs) .or. .not.present(mask)) then + if(myrank.eq.0) then + print *, "ERROR: field_from_coils missing some optional arguments" + end if + call safestop(1) + end if + numnodes = owned_nodes() + do icounter_t=1,numnodes + inode = nodes_owned(icounter_t) + call get_node_data(f, inode, dataout) + do k=1, ncs + cur = sum(ic, mask=(mask.eq.k)) + call get_node_data(fs(k), inode, datain) + datain = datain*cur + dataout = dataout + datain + end do + call set_node_data(f, inode, dataout) + enddo + call finalize(f%vec) + + return + end if + + allocate(g(MAX_PTS,nc)) call create_field(psi_vec) psi_vec = 0. + + numelms = local_elements() + do itri=1,numelms + temp79a = 0. + + ! calculate field due to coil currents + call define_element_quadrature(itri,int_pts_main,int_pts_tor) + call define_fields(itri,0,1,0) + call gvect0(x_79,z_79,npoints,xc,zc,nc,g(1:npoints,:),ipole,ierr) + do k=1, nc + temp79a = temp79a - g(:,k)*ic(k) + end do + + dofs = intx2(mu79(:,:,OP_1),temp79a) + + call vector_insert_block(psi_vec%vec, itri, 1, dofs, MAT_ADD) + enddo + + if(allocated(g)) deallocate(g) + + call sum_shared(psi_vec%vec) + call newsolve(mass_mat_lhs%mat,psi_vec%vec,ierr) + call add(f, psi_vec) + call destroy_field(psi_vec) + + end subroutine field_from_coils + + subroutine store_field_from_coils(xc, zc, nc, nc_base, fs, mask, ipole, ierr) + use basic + use field + use mesh_mod + use m3dc1_nint + use newvar_mod + + implicit none + + include 'mpif.h' + + real, intent(in), dimension(nc) :: xc, zc ! array of coil positions + integer, intent(in) :: nc ! number of filaments + integer, intent(in) :: nc_base ! number of coils + integer, dimension(maxfilaments), intent(in) :: mask ! relate filaments to coils + type(field_type), allocatable, intent(inout) :: fs(:) ! poloidal flux field for each coil + integer, intent(in) :: ipole ! type of field to add + integer, intent(out) :: ierr + + integer :: numelms, k, j, itri, nfil + vectype, dimension(dofs_per_element) :: dofs + real, allocatable :: g(:,:) + + if(allocated(fs)) then + if (myrank.eq.0) print *,"ERROR: store_field_from_coils called twice somehow" + call safestop(1) + end if + + if(nc.le.0) return + + ierr = 0 + + allocate(fs(nc_base)) + allocate(g(MAX_PTS,nc)) + + do j=1,nc_base + call create_field(fs(j)) + enddo numelms = local_elements() do itri=1,numelms @@ -196,25 +296,30 @@ subroutine field_from_coils(xc, zc, ic, nc, f, ipole, ierr) call define_fields(itri,0,1,0) ! Field due to coil currents - temp79a = 0. call gvect0(x_79,z_79,npoints,xc,zc,nc,g(1:npoints,:),ipole,ierr) - do k=1, nc - temp79a = temp79a - g(:,k)*ic(k) - end do - - dofs = intx2(mu79(:,:,OP_1),temp79a) - - call vector_insert_block(psi_vec%vec, itri, 1, dofs, MAT_ADD) + do j=1, nc_base + temp79a = 0. + nfil = 0 + do k=1, nc + if (mask(k).eq.j) then + temp79a = temp79a - g(:,k) + nfil = nfil + 1 + end if + end do + temp79a = temp79a/nfil + dofs = intx2(mu79(:,:,OP_1), temp79a) + call vector_insert_block(fs(j)%vec, itri, 1, dofs, MAT_ADD) + enddo enddo deallocate(g) - call sum_shared(psi_vec%vec) - call newsolve(mass_mat_lhs%mat,psi_vec%vec,ierr) - call add(f, psi_vec) + do j=1,nc_base + call sum_shared(fs(j)%vec) + call newsolve(mass_mat_lhs%mat,fs(j)%vec,ierr) + enddo - call destroy_field(psi_vec) - end subroutine field_from_coils + end subroutine store_field_from_coils !============================================================ subroutine gvect(r,z,npt,xi,zi,n,g,nmult,ierr) diff --git a/unstructured/gradshafranov.f90 b/unstructured/gradshafranov.f90 index 8a74a6a8a..75d484054 100644 --- a/unstructured/gradshafranov.f90 +++ b/unstructured/gradshafranov.f90 @@ -38,10 +38,11 @@ module gradshafranov real, dimension(maxfilaments) :: xc_vac, zc_vac complex, dimension(maxfilaments) :: ic_vac, ic_out - integer :: numcoils_vac + integer :: numcoils_vac, numcoils_pf integer, dimension(maxfilaments) :: coil_mask integer, dimension(maxfilaments) :: filaments + integer :: igs_store_coils real, dimension(maxcoils) :: gs_vertical_feedback real, dimension(maxcoils) :: gs_radial_feedback real, dimension(maxcoils) :: gs_vertical_feedback_i @@ -184,8 +185,18 @@ subroutine coil_feedback(itnum) if(myrank.eq.0 .and. iprint.ge.1) & print *, "Calculating fields due to coils in feedback loop" psi_coil_field = 0. - call field_from_coils(xc_vac,zc_vac,ic_out,numcoils_vac, & - psi_coil_field,0,ierr) + + if(igs_store_coils.eq.0) then + call field_from_coils(xc_vac,zc_vac,ic_out,numcoils_vac,psi_coil_field,0,ierr) + elseif(igs_store_coils.eq.1) then + if(.not.allocated(psi_coil_fields)) then + call store_field_from_coils(xc_vac,zc_vac,numcoils_vac,numcoils_pf,& + psi_coil_fields,coil_mask,0,ierr) + end if + call field_from_coils(xc_vac,zc_vac,ic_out,numcoils_vac,psi_coil_field,0,ierr,& + numcoils_pf,psi_coil_fields,coil_mask) + end if + if(myrank.eq.0 .and. iprint.ge.2) & print *, "Done calculating fields due to coils" if(ierr.ne.0) call safestop(5) @@ -260,7 +271,7 @@ subroutine pf_coil_field(ierr) select case(idevice) case(-1) call load_coils(xc,zc,ic,numcoils,'coil.dat','current.dat',coil_mask,& - filaments) + filaments,numcoils_pf) xc_vac = xc zc_vac = zc ic_vac = ic @@ -284,7 +295,16 @@ subroutine pf_coil_field(ierr) ! Field due to coil currents if(myrank.eq.0 .and. iprint.ge.1) & print *, "Calculating fields due to coils" - call field_from_coils(xc,zc,ic,numcoils,psi_coil_field,ipole,ierr) + if(igs_store_coils.eq.0) then + call field_from_coils(xc,zc,ic,numcoils,psi_coil_field,ipole,ierr) + elseif(igs_store_coils.eq.1) then + if(.not.allocated(psi_coil_fields)) then + call store_field_from_coils(xc,zc,numcoils,numcoils_pf,& + psi_coil_fields,coil_mask,ipole,ierr) + end if + call field_from_coils(xc,zc,ic,numcoils,psi_coil_field,ipole,ierr,& + numcoils_pf,psi_coil_fields,coil_mask) + end if if(myrank.eq.0 .and. iprint.ge.1) & print *, "Done calculating fields due to coils" if(ierr.ne.0) call safestop(5) diff --git a/unstructured/input.f90 b/unstructured/input.f90 index 729b3a398..0cf827a98 100644 --- a/unstructured/input.f90 +++ b/unstructured/input.f90 @@ -656,6 +656,8 @@ subroutine set_defaults "1: read profile_pscale for factor to scale p and p'", gs_grp) call add_var_double("vscale", vscale, 1., & "Factor multiplying toroidal rotation profile", gs_grp) + call add_var_int("igs_store_coils", igs_store_coils, 0, & + "Store coil fields in GS solve", gs_grp) call add_var_double_array("gs_vertical_feedback", gs_vertical_feedback, & maxcoils, 0., & "Proportional feedback of each coil to vertical displacements", gs_grp)