test/modules/LinearBoltzmannSolvers/Transport_Steady_Cyl/Transport2DCyl_1Monoenergetic.lua

-- 2D transport test in axialsymmetric cylindrical geometry with
-- vacuum boundary condition - monoenergetic.
-- SDM: PWLD
-- Test: Max-value=1.00000
num_procs = 4
--Structured mesh
 
 
 
 
--############################################### Check num_procs
if (check_num_procs==nil and chi_number_of_processes ~= num_procs) then
    chiLog(LOG_0ERROR,"Incorrect amount of processors. " ..
                      "Expected "..tostring(num_procs)..
                      ". Pass check_num_procs=false to override if possible.")
    os.exit(false)
end
 
--############################################### Setup mesh
dim = 2
length = {1.0, 2.0, }
ncells = {50, 100, }
nodes = {}
for d = 1, dim do
  delta = length[d]/ncells[d]
  nodes[d] = {}
  for i = 0, ncells[d] do
    nodes[d][i+1] = i*delta
  end
end
 
meshgen1 = chi_mesh.OrthogonalMeshGenerator.Create({ node_sets = {nodes[1],nodes[2]} })
chi_mesh.MeshGenerator.Execute(meshgen1)
 
--############################################### Set Material IDs
vol0 = chi_mesh.RPPLogicalVolume.Create
({ xmin=0.0,xmax=length[1],ymin=0.0,ymax=length[2], infz=true })
chiVolumeMesherSetProperty(MATID_FROMLOGICAL,vol0,0)
 
--############################################### Add materials
ngrp = 1
sigmat = 25.0
ratioc = 0.1
source = sigmat * (1 - ratioc)
 
material0 = chiPhysicsAddMaterial("Material_0");
chiPhysicsMaterialAddProperty(material0,TRANSPORT_XSECTIONS)
chiPhysicsMaterialAddProperty(material0,ISOTROPIC_MG_SOURCE)
chiPhysicsMaterialSetProperty(material0, TRANSPORT_XSECTIONS,
                              SIMPLEXS1, ngrp, sigmat, ratioc)
chiPhysicsMaterialSetProperty(material0, ISOTROPIC_MG_SOURCE,
                              SINGLE_VALUE, source)
 
--############################################### Setup Physics
pquad0 = chiCreateCylindricalProductQuadrature(GAUSS_LEGENDRE_CHEBYSHEV, 4, 8)
 
lbs_block =
{
  coord_system = 2,
  num_groups = ngrp,
  groupsets =
  {
    {
      groups_from_to = {0, ngrp-1},
      angular_quadrature_handle = pquad0,
      angle_aggregation_type = "azimuthal",
      inner_linear_method = "gmres",
      l_max_its = 100,
      l_abs_tol = 1.0e-12,
    }
  }
}
 
lbs_options =
{
  boundary_conditions = { { name = "xmin", type = "reflecting"} },
  scattering_order = 0,
}
phys1 = lbs.DiscreteOrdinatesCurvilinearSolver.Create(lbs_block)
lbs.SetOptions(phys1, lbs_options)
 
--############################################### Initialize and Execute Solver
ss_solver = lbs.SteadyStateSolver.Create({lbs_solver_handle = phys1})
 
chiSolverInitialize(ss_solver)
chiSolverExecute(ss_solver)
 
 
--phys0 = chiLBSCurvilinearCreateSolver(LBSCurvilinear.CYLINDRICAL)
--
----  angular quadrature
--pquad = chiCreateCylindricalProductQuadrature(GAUSS_LEGENDRE_CHEBYSHEV, 4, 8)
--
----  groups
--groups = {}
--for g = 1, ngrp do
--  groups[g] = chiLBSCreateGroup(phys0)
--end
--
----  groupsets
--gs0 = chiLBSCreateGroupset(phys0)
--chiLBSGroupsetAddGroups(phys0, gs0, 0, ngrp-1)
--chiLBSGroupsetSetQuadrature(phys0, gs0, pquad)
--chiLBSGroupsetSetAngleAggregationType(phys0, gs0, LBSGroupset.ANGLE_AGG_AZIMUTHAL)
--chiLBSGroupsetSetIterativeMethod(phys0, gs0, KRYLOV_GMRES_CYCLES)
--chiLBSGroupsetSetResidualTolerance(phys0, gs0, 1.0e-12)
--chiLBSGroupsetSetMaxIterations(phys0, gs0, 100)
--chiLBSGroupsetSetGMRESRestartIntvl(phys0, gs0, 30)
--
----  spatial discretisation
--chiLBSSetProperty(phys0, DISCRETIZATION_METHOD, PWLD)
--
----  scattering order
--chiLBSSetProperty(phys0, SCATTERING_ORDER, 0)
--
----------------------------------------------------------------------------------
----  boundary conditions
----------------------------------------------------------------------------------
--dirichlet_value = {}
--for g = 1, ngrp do
--  dirichlet_value[g] = 0
--end
--chiLBSSetProperty(phys0, BOUNDARY_CONDITION,
--                  XMIN, LBSBoundaryTypes.REFLECTING)
--chiLBSSetProperty(phys0, BOUNDARY_CONDITION,
--                  XMAX, LBSBoundaryTypes.INCIDENT_ISOTROPIC, dirichlet_value)
--chiLBSSetProperty(phys0, BOUNDARY_CONDITION,
--                  YMIN, LBSBoundaryTypes.INCIDENT_ISOTROPIC, dirichlet_value)
--chiLBSSetProperty(phys0, BOUNDARY_CONDITION,
--                  YMAX, LBSBoundaryTypes.INCIDENT_ISOTROPIC, dirichlet_value)
--
----------------------------------------------------------------------------------
----  solvers
----------------------------------------------------------------------------------
--chiSolverInitialize(phys0)
--chiSolverExecute(phys0)
 
--############################################### Exports
fflist, count = chiLBSGetScalarFieldFunctionList(phys1)
if master_export == nil then
  chiExportMultiFieldFunctionToVTK(fflist, "ZRZPhi")
end
 
--############################################### Volume integrations
ffi1 = chiFFInterpolationCreate(VOLUME)
curffi = ffi1
chiFFInterpolationSetProperty(curffi, OPERATION, OP_MAX)
chiFFInterpolationSetProperty(curffi, LOGICAL_VOLUME, vol0)
chiFFInterpolationSetProperty(curffi, ADD_FIELDFUNCTION, fflist[1])
 
chiFFInterpolationInitialize(curffi)
chiFFInterpolationExecute(curffi)
maxval = chiFFInterpolationGetValue(curffi)
 
chiLog(LOG_0, string.format("Max-value=%.5f", maxval))