TokaMaker: Coupled TokaMaker/MUG time-dependent evolution (OpenFUSIONToolkit#320)

Author: sguizzo

src/examples/TokaMaker/AdvancedWorkflows/Blanket/blanket_mesh.h5

@@ -304,20 +304,22 @@ subroutine fem_vec_load(self,source,filename,path,err_flag)
 NULLIFY(lge)
 DEBUG_STACK_POP
 end subroutine fem_vec_load
+
 !------------------------------------------------------------------------------
-!> Needs docs
+!> Generate graph structure for composite FE representation 
 !------------------------------------------------------------------------------
-subroutine fem_mat_create(self,new,mask)
+subroutine fem_graph_create(self, graphs, known_graphs, nknown_graphs, mask)
 CLASS(oft_fem_comp_type), INTENT(inout) :: self
-CLASS(oft_matrix), POINTER, INTENT(out) :: new
+TYPE(oft_graph_ptr), ALLOCATABLE, INTENT(out) :: graphs(:,:)
+TYPE(oft_graph_ptr), ALLOCATABLE, INTENT(out) :: known_graphs(:)
+INTEGER(i4), INTENT(out) :: nknown_graphs
 INTEGER(i4), OPTIONAL, INTENT(in) :: mask(:,:)
-INTEGER(i4) :: i,j,k,nknown_graphs
+INTEGER(i4) :: i,j,k
 INTEGER(i4), ALLOCATABLE, DIMENSION(:,:) :: mat_mask,graph_ids
 CLASS(oft_vector), POINTER :: tmp_vec
-TYPE(oft_graph_ptr), ALLOCATABLE :: graphs(:,:),known_graphs(:)
 DEBUG_STACK_PUSH
 !---
-IF(oft_debug_print(2))WRITE(*,'(2X,A)')'Building composite FE matrix'
+IF(oft_debug_print(2))WRITE(*,'(2X,A)')'Building composite FE graph'
 ALLOCATE(mat_mask(self%nfields,self%nfields))
 mat_mask=1
 IF(PRESENT(mask))mat_mask=mask
@@ -351,7 +353,7 @@ subroutine fem_mat_create(self,new,mask)
     IF(mat_mask(i,j)==0)CYCLE
     IF(mat_mask(i,j)==2)THEN
       IF(i/=j)CALL oft_abort('Identity only valid on diagonal.', &
-      'fem_mat_create',__FILE__)
+      'fem_graph_create',__FILE__)
       !---Setup identity graph
       CALL self%fields(i)%fe%vec_create(tmp_vec)
       CALL create_identity_graph(graphs(i,j)%g,tmp_vec)
@@ -383,14 +385,41 @@ subroutine fem_mat_create(self,new,mask)
     END IF
   END DO
 END DO
+DEALLOCATE(mat_mask,graph_ids)
+DEBUG_STACK_POP
+end subroutine fem_graph_create
+
+!------------------------------------------------------------------------------
+!> Generate composite FE matrix from user-provided (or default) graph structure
+!------------------------------------------------------------------------------
+subroutine fem_mat_create(self,new,mask, graphs_in)
+CLASS(oft_fem_comp_type), INTENT(inout) :: self
+CLASS(oft_matrix), POINTER, INTENT(out) :: new
+INTEGER(i4), OPTIONAL, INTENT(in) :: mask(:,:)
+TYPE(oft_graph_ptr), OPTIONAL, INTENT(in) :: graphs_in(:,:)
+INTEGER(i4) :: i,nknown_graphs
+CLASS(oft_vector), POINTER :: tmp_vec
+TYPE(oft_graph_ptr), ALLOCATABLE :: graphs(:,:), known_graphs(:)
+DEBUG_STACK_PUSH
+!---
+IF (PRESENT(graphs_in)) THEN
+  IF(oft_debug_print(2))WRITE(*,'(2X,A)')'Using user-provided graphs for FE matrix'
+  graphs = graphs_in
+ELSE
+  IF(oft_debug_print(2))WRITE(*,'(2X,A)')'Building composite FE matrix'
+  CALL fem_graph_create(self, graphs, known_graphs, nknown_graphs, mask)
+END IF
 !---
 CALL self%vec_create(tmp_vec)
 CALL create_matrix(new,graphs,tmp_vec,tmp_vec)
 CALL tmp_vec%delete
-DO i=1,nknown_graphs
-  DEALLOCATE(known_graphs(i)%g)
-END DO
-DEALLOCATE(graphs,known_graphs,mat_mask,graph_ids,tmp_vec)
+IF (.NOT. PRESENT(graphs_in)) THEN
+  DO i=1,nknown_graphs
+    DEALLOCATE(known_graphs(i)%g)
+  END DO
+  DEALLOCATE(known_graphs)
+END IF
+DEALLOCATE(graphs,tmp_vec)
 DEBUG_STACK_POP
 end subroutine fem_mat_create
 !------------------------------------------------------------------------------

src/examples/TokaMaker/AdvancedWorkflows/Blanket/vde_blanket.ipynb

@@ -2590,7 +2590,7 @@ subroutine native_matrix_setup_full(self,u)
     jp=1
     last=0
     do mm=1,netmp ! Loop over input points
-      IF(letmp(1,mm)==last.AND.periodic)THEN
+      IF(letmp(1,mm)==last)THEN
         m=jp
       ELSE
         last=letmp(1,mm)
@@ -2631,7 +2631,7 @@ subroutine native_matrix_setup_full(self,u)
   jp=1
   last=0
   do mm=1,neout ! Loop over input points
-    IF(leout(1,mm)==last.AND.periodic)THEN
+    IF(leout(1,mm)==last)THEN
       m=jp
     ELSE
       last=leout(1,mm)

src/fem/fem_composite.F90

@@ -20,6 +20,7 @@ set( OBJS
   thin_wall.F90
   thin_wall_solvers.F90
   thin_wall_hodlr.F90
+  mugtok_td.F90
 )
 
 add_library( oftphysics ${OBJS} )

src/lin_alg/native_la.F90

@@ -188,6 +188,7 @@ MODULE oft_gs
   LOGICAL, POINTER, DIMENSION(:) :: saddle_pmask => NULL() !< Point mask for saddle search
   LOGICAL, POINTER, DIMENSION(:) :: saddle_cmask => NULL() !< Cell mask for saddle search
   LOGICAL, POINTER, DIMENSION(:) :: saddle_rmask => NULL() !< Region mask for saddle search
+  LOGICAL, POINTER, DIMENSION(:) :: ignore_rmask => NULL() !< Mask for regions to ignore when populating physics (defaults to false)
   INTEGER(i4), POINTER, DIMENSION(:) :: limiter_nds => NULL() !< List of limiter nodes
   INTEGER(i4), POINTER, DIMENSION(:) :: bc_rhs_list => NULL() !< List of terms interacting with free-boundary BC
   INTEGER(i4), POINTER, DIMENSION(:) :: olbp => NULL() !< Oriented list of boundary points
@@ -809,6 +810,10 @@ subroutine gs_init(self)
   END DO
 END DO
 !---Build operators
+IF(.NOT.ASSOCIATED(self%ignore_rmask))THEN
+  ALLOCATE(self%ignore_rmask(smesh%nreg))
+  self%ignore_rmask=.FALSE.
+END IF
 IF(.NOT.ASSOCIATED(self%dels))THEN
   IF(self%free)THEN
     CALL compute_bcmat(self)
@@ -5462,6 +5467,8 @@ subroutine build_dels(mat,self,bc,dt,scale)
 IF(nnonaxi>0)allocate(nonaxi_tmp(self%fe_rep%nce))
 !$omp do schedule(dynamic,1) ordered
 do i=1,self%fe_rep%mesh%nc
+  !---Skip cell if in 'ignore' region
+  IF(self%ignore_rmask(self%fe_rep%mesh%reg(i))) CYCLE
   !---Get local reconstructed operators
   lop=0.d0
   IF(nnonaxi>0)nonaxi_tmp=0.d0
@@ -5668,6 +5675,8 @@ subroutine dels_coil_part(self,mat,iCoil,dt_in,bc,main_scale,nonaxi_vals)
 allocate(j_lag(self%fe_rep%nce),col_tmp(1,self%fe_rep%nce))
 !!$omp do schedule(static)
 DO j=1,smesh%nc
+  !---Skip cell if in 'ignore' region
+  IF(self%ignore_rmask(self%fe_rep%mesh%reg(j))) CYCLE
   nturns=self%coil_nturns(smesh%reg(j),iCoil)
   eta_wt=0.d0
   IF(eta_reg(smesh%reg(j))>0.d0)eta_wt=1.d0/(dt_in*eta_reg(smesh%reg(j)))

src/physics/CMakeLists.txt

@@ -538,6 +538,8 @@ subroutine apply_rhs(self,a,b)
 allocate(rop(lag_rep%nce),gop(3,lag_rep%nce)) ! Reconstructed gradient operator
 !$omp do schedule(static)
 do i=1,mesh%nc
+    !---Skip cell if in 'ignore' region
+    IF(self%gs_device%ignore_rmask(mesh%reg(i))) CYCLE
     IF(mesh%reg(i)==1)CYCLE
     !---Get local to global DOF mapping
     call lag_rep%ncdofs(i,j)