Store coil fields for GS solve

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

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)

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)

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)