diff --git a/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/2d_INSAD.i b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/2d_INSAD.i new file mode 100644 index 0000000..4bd2d97 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/2d_INSAD.i @@ -0,0 +1,176 @@ +# Note: if this is run on 1 core only, ILU is not used +# And instead a method requiring more iterations is applied. + +N_X = 200 +N_Y_half = 10 +U_AVG = 1 +element_type = 'QUAD9' + +[Mesh] + [meshTop] + type = GeneratedMeshGenerator + dim = 2 + nx = ${N_X} + ny = ${N_Y_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + bias_y = 0.8 + boundary_name_prefix = 'meshTop' + elem_type = ${element_type} + [] + [meshBottom] + type = GeneratedMeshGenerator + dim = 2 + nx = ${N_X} + ny = ${N_Y_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + bias_y = 1.25 + boundary_name_prefix = 'meshBottom' + elem_type = ${element_type} + [] + [meshComplete] + type = StitchedMeshGenerator + inputs = 'meshTop meshBottom' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTop_bottom meshBottom_top' + [] + [meshRename] + type = RenameBoundaryGenerator + input = meshComplete + old_boundary = ' + meshTop_right meshBottom_right + meshTop_left meshBottom_left + meshTop_top + meshBottom_bottom + ' + new_boundary = ' + right right + left left + top + bottom + ' + [] +[] + +[Variables] + [velocity] + family = LAGRANGE_VEC + order = SECOND + [] + [pressure] + family = MONOMIAL + order = FIRST + [] +[] + +# [ICs] +# [velocity] +# type = VectorFunctionIC +# variable = velocity +# function = velocityInlet +# [] +# [] + +[Functions] + [velocityInlet] + type = ParsedVectorFunction + symbol_names = 'y_max' + symbol_values = '1' + expression_x = '(3/2) * ${U_AVG} * (1 - (y * y) / (y_max * y_max))' + expression_y = '0' + [] +[] + +[BCs] + [inlet] + type = VectorFunctionDirichletBC + variable = velocity + boundary = left + function = velocityInlet + [] + [no_slip] + type = VectorDirichletBC + variable = velocity + boundary = 'top bottom' + values = '0 0 0' + [] + [pressure_set] + type = PenaltyDirichletBC + variable = pressure + value = 0 + boundary = 'right' + penalty = 1e5 + [] + # [pressure_set] + # type = DirichletBC + # variable = pressure + # boundary = 'right' + # value = 0 + # [] +[] + +[Materials] + [const] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '1 1' + [] + [ins_mat] + type = INSADTauMaterial + velocity = velocity + pressure = pressure + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = pressure + [] + + [momentum_convection] + type = INSADMomentumAdvection + variable = velocity + [] + + [momentum_viscous] + type = INSADMomentumViscous + variable = velocity + [] + + [momentum_pressure] + type = INSADMomentumPressure + variable = velocity + pressure = pressure + integrate_p_by_parts = true + [] +[] + +[Problem] + type = FEProblem +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + l_max_its = 100 + nl_max_its = 1000 + petsc_options_iname = '-pc_type' + petsc_options_value = 'ilu' +[] + +[Outputs] + exodus = true +[] diff --git a/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INS.i b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INS.i new file mode 100644 index 0000000..1c56622 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INS.i @@ -0,0 +1,320 @@ +U_AVG = 1 + +N_X = 50 +N_Y_half = 10 +N_Z_half = 10 +ELEMENT_TYPE = HEX27 + +GRADING_R_Y = 5 +GRADING_R_Z = 5 +RATIO_Y_FWD = ${fparse GRADING_R_Y ^ (1/(N_Y_half - 1))} +RATIO_Y_INV = ${fparse 1/RATIO_Y_FWD} +RATIO_Z_FWD = ${fparse GRADING_R_Z ^ (1/(N_Z_half - 1))} +RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} + +[Mesh] + [meshTopBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopBack' + [] + [meshTopFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopFront' + [] + [meshBottomBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomBack' + [] + [meshBottomFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomFront' + [] + [meshTop] + type = StitchedMeshGenerator + inputs = 'meshTopBack meshTopFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTopBack_front meshTopFront_back' + [] + [renameMeshTop] + type = RenameBoundaryGenerator + input = meshTop + old_boundary = ' + meshTopBack_top meshTopFront_top + meshTopBack_right meshTopFront_right + meshTopBack_left meshTopFront_left + meshTopBack_bottom meshTopFront_bottom + ' + new_boundary = ' + top top + meshTop_right meshTop_right + meshTop_left meshTop_left + meshTop_bottom meshTop_bottom + ' + [] + [meshBottom] + type = StitchedMeshGenerator + inputs = 'meshBottomBack meshBottomFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshBottomBack_front meshBottomFront_back' + [] + [renameMeshBottom] + type = RenameBoundaryGenerator + input = meshBottom + old_boundary = ' + meshBottomBack_top meshBottomFront_top + meshBottomBack_right meshBottomFront_right + meshBottomBack_left meshBottomFront_left + meshBottomBack_bottom meshBottomFront_bottom + ' + new_boundary = ' + meshBottom_top meshBottom_top + meshBottom_right meshBottom_right + meshBottom_left meshBottom_left + bottom bottom + ' + [] + [mesh] + type = StitchedMeshGenerator + inputs = 'renameMeshTop renameMeshBottom' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTop_bottom meshBottom_top' + [] + [renameMesh] + type = RenameBoundaryGenerator + input = mesh + old_boundary = ' + meshBottomFront_front meshTopFront_front + meshBottomBack_back meshTopBack_back + meshBottom_left meshTop_left + meshBottom_right meshTop_right + ' + new_boundary = ' + front front + back back + left left + right right + ' + [] +[] + +[Variables] + [velocity_x] + family = LAGRANGE + order = SECOND + [] + [velocity_y] + family = LAGRANGE + order = SECOND + [] + [velocity_z] + family = LAGRANGE + order = SECOND + [] + [pressure] + family = MONOMIAL + order = FIRST + [] +[] + +[ICs] + [velocityIC_x] + type = ConstantIC + value = ${U_AVG} + variable = velocity_x + [] + [velocityIC_y] + type = ConstantIC + value = 1e-15 + variable = velocity_y + [] + [velocityIC_z] + type = ConstantIC + value = 1e-15 + variable = velocity_z + [] +[] + +[Functions] + [velocityInlet] + type = ParsedFunction + symbol_names = 'y_max z_max' + symbol_values = '1 1' + expression = '(9/4) * ${U_AVG} * (1 - (y * y) / (y_max * y_max)) * (1 - (z * z) / (z_max * z_max))' + [] +[] + +[BCs] + [inlet_x] + type = FunctionDirichletBC + variable = velocity_x + boundary = left + function = velocityInlet + [] + [inlet_y] + type = DirichletBC + variable = velocity_y + boundary = left + value = 0.0 + [] + [inlet_z] + type = DirichletBC + variable = velocity_z + boundary = left + value = 0.0 + [] + [no_slip_x] + type = DirichletBC + variable = velocity_x + boundary = 'top bottom front back' + value = 0.0 + [] + [no_slip_y] + type = DirichletBC + variable = velocity_y + boundary = 'top bottom front back' + value = 0.0 + [] + [no_slip_z] + type = DirichletBC + variable = velocity_z + boundary = 'top bottom front back' + value = 0.0 + [] + # [pressure_set] + # type = DirichletBC + # variable = pressure + # boundary = right + # value = 0.0 + # [] + [pressure_set] + type = PenaltyDirichletBC + variable = pressure + value = 0 + boundary = 'right' + penalty = 1e5 + [] +[] + +[Materials] + [const] + type = GenericConstantMaterial + prop_names = 'rho mu' + prop_values = '1 1' + [] +[] + +[Kernels] + [mass] + type = INSMass + variable = pressure + u = velocity_x + v = velocity_y + w = velocity_z + pressure = pressure + [] + [x_momentum_space] + type = INSMomentumLaplaceForm + variable = velocity_x + u = velocity_x + v = velocity_y + w = velocity_z + pressure = pressure + component = 0 + [] + [y_momentum_space] + type = INSMomentumLaplaceForm + variable = velocity_y + u = velocity_x + v = velocity_y + w = velocity_z + pressure = pressure + component = 1 + [] + [z_momentum_space] + type = INSMomentumLaplaceForm + variable = velocity_z + u = velocity_x + v = velocity_y + w = velocity_z + pressure = pressure + component = 2 + [] +[] + +[Problem] + type = FEProblem +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + automatic_scaling = true + l_max_its = 1000 + nl_max_its = 1000 + petsc_options_iname = '-pc_type' + petsc_options_value = 'bjacobi' +[] + +[Outputs] + exodus = true + execute_on = 'nonlinear' +[] diff --git a/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INSAD.i b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INSAD.i new file mode 100644 index 0000000..ffc5d0d --- /dev/null +++ b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INSAD.i @@ -0,0 +1,266 @@ +U_AVG = 1 + +N_X = 50 +N_Y_half = 10 +N_Z_half = 10 +ELEMENT_TYPE = HEX27 + +GRADING_R_Y = 5 +GRADING_R_Z = 5 +RATIO_Y_FWD = ${fparse GRADING_R_Y ^ (1/(N_Y_half - 1))} +RATIO_Y_INV = ${fparse 1/RATIO_Y_FWD} +RATIO_Z_FWD = ${fparse GRADING_R_Z ^ (1/(N_Z_half - 1))} +RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} + +[Mesh] + [meshTopBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopBack' + [] + [meshTopFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopFront' + [] + [meshBottomBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomBack' + [] + [meshBottomFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomFront' + [] + [meshTop] + type = StitchedMeshGenerator + inputs = 'meshTopBack meshTopFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTopBack_front meshTopFront_back' + [] + [renameMeshTop] + type = RenameBoundaryGenerator + input = meshTop + old_boundary = ' + meshTopBack_top meshTopFront_top + meshTopBack_right meshTopFront_right + meshTopBack_left meshTopFront_left + meshTopBack_bottom meshTopFront_bottom + ' + new_boundary = ' + top top + meshTop_right meshTop_right + meshTop_left meshTop_left + meshTop_bottom meshTop_bottom + ' + [] + [meshBottom] + type = StitchedMeshGenerator + inputs = 'meshBottomBack meshBottomFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshBottomBack_front meshBottomFront_back' + [] + [renameMeshBottom] + type = RenameBoundaryGenerator + input = meshBottom + old_boundary = ' + meshBottomBack_top meshBottomFront_top + meshBottomBack_right meshBottomFront_right + meshBottomBack_left meshBottomFront_left + meshBottomBack_bottom meshBottomFront_bottom + ' + new_boundary = ' + meshBottom_top meshBottom_top + meshBottom_right meshBottom_right + meshBottom_left meshBottom_left + bottom bottom + ' + [] + [mesh] + type = StitchedMeshGenerator + inputs = 'renameMeshTop renameMeshBottom' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTop_bottom meshBottom_top' + [] + [renameMesh] + type = RenameBoundaryGenerator + input = mesh + old_boundary = ' + meshBottomFront_front meshTopFront_front + meshBottomBack_back meshTopBack_back + meshBottom_left meshTop_left + meshBottom_right meshTop_right + ' + new_boundary = ' + front front + back back + left left + right right + ' + [] +[] + +[Variables] + [velocity] + family = LAGRANGE_VEC + order = SECOND + [] + [pressure] + family = MONOMIAL + order = FIRST + [] +[] + +[ICs] + [velocityIC] + type = VectorConstantIC + x_value = ${U_AVG} + y_value = 1e-15 + z_value = 1e-15 + variable = velocity + [] +[] + +[Functions] + [velocityInlet] + type = ParsedVectorFunction + symbol_names = 'y_max z_max' + symbol_values = '1 1' + expression_x = '(9/4) * ${U_AVG} * (1 - (y * y) / (y_max * y_max)) * (1 - (z * z) / (z_max * z_max))' + expression_y = '0' + expression_z = '0' + [] +[] + +[BCs] + [inlet] + type = VectorFunctionDirichletBC + variable = velocity + boundary = left + function = velocityInlet + [] + [no_slip] + type = VectorDirichletBC + variable = velocity + boundary = 'top bottom front back' + values = '0 0 0' + [] + [pressure_set] + type = PenaltyDirichletBC + variable = pressure + value = 0 + boundary = 'right' + penalty = 1e5 + [] +[] + +[Materials] + [const] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '1 1' + [] + [ins_mat] + type = INSADTauMaterial + velocity = velocity + pressure = pressure + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = pressure + [] + + [momentum_convection] + type = INSADMomentumAdvection + variable = velocity + [] + + [momentum_viscous] + type = INSADMomentumViscous + variable = velocity + [] + + [momentum_pressure] + type = INSADMomentumPressure + variable = velocity + pressure = pressure + integrate_p_by_parts = true + [] +[] + +[Problem] + type = FEProblem +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + automatic_scaling = true + l_max_its = 1000 + nl_max_its = 1000 + petsc_options_iname = '-pc_type' + petsc_options_value = 'bjacobi' +[] + +[Outputs] + exodus = true +[] diff --git a/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i b/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i new file mode 100644 index 0000000..5587825 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i @@ -0,0 +1,301 @@ +U_AVG = 1 + +N_X = 50 +N_Y_half = 10 +N_Z_half = 10 +ELEMENT_TYPE = HEX27 + +GRADING_R_Y = 5 +GRADING_R_Z = 5 +RATIO_Y_FWD = ${fparse GRADING_R_Y ^ (1/(N_Y_half - 1))} +RATIO_Y_INV = ${fparse 1/RATIO_Y_FWD} +RATIO_Z_FWD = ${fparse GRADING_R_Z ^ (1/(N_Z_half - 1))} +RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} + +[Mesh] + [meshTopBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopBack' + [] + [meshTopFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopFront' + [] + [meshBottomBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomBack' + [] + [meshBottomFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomFront' + [] + [meshTop] + type = StitchedMeshGenerator + inputs = 'meshTopBack meshTopFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTopBack_front meshTopFront_back' + [] + [renameMeshTop] + type = RenameBoundaryGenerator + input = meshTop + old_boundary = ' + meshTopBack_top meshTopFront_top + meshTopBack_right meshTopFront_right + meshTopBack_left meshTopFront_left + meshTopBack_bottom meshTopFront_bottom + ' + new_boundary = ' + top top + meshTop_right meshTop_right + meshTop_left meshTop_left + meshTop_bottom meshTop_bottom + ' + [] + [meshBottom] + type = StitchedMeshGenerator + inputs = 'meshBottomBack meshBottomFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshBottomBack_front meshBottomFront_back' + [] + [renameMeshBottom] + type = RenameBoundaryGenerator + input = meshBottom + old_boundary = ' + meshBottomBack_top meshBottomFront_top + meshBottomBack_right meshBottomFront_right + meshBottomBack_left meshBottomFront_left + meshBottomBack_bottom meshBottomFront_bottom + ' + new_boundary = ' + meshBottom_top meshBottom_top + meshBottom_right meshBottom_right + meshBottom_left meshBottom_left + bottom bottom + ' + [] + [mesh] + type = StitchedMeshGenerator + inputs = 'renameMeshTop renameMeshBottom' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTop_bottom meshBottom_top' + [] + [renameMesh] + type = RenameBoundaryGenerator + input = mesh + old_boundary = ' + meshBottomFront_front meshTopFront_front + meshBottomBack_back meshTopBack_back + meshBottom_left meshTop_left + meshBottom_right meshTop_right + ' + new_boundary = ' + front front + back back + left left + right right + ' + [] +[] + +[Variables] + [velocity] + family = LAGRANGE_VEC + order = SECOND + [] + [pressure] + family = MONOMIAL + order = FIRST + [] +[] + +[AuxVariables] + [currentDensity] + family = LAGRANGE_VEC + order = FIRST + [] +[] + +[ICs] + [velocityIC] + type = VectorConstantIC + x_value = ${U_AVG} + y_value = 1e-15 + z_value = 1e-15 + variable = velocity + [] +[] + +[Functions] + [velocityInlet] + type = ParsedVectorFunction + symbol_names = 'y_max z_max' + symbol_values = '1 1' + expression_x = '(9/4) * ${U_AVG} * (1 - (y * y) / (y_max * y_max)) * (1 - (z * z) / (z_max * z_max))' + expression_y = '0' + expression_z = '0' + [] + [magneticFieldFunction] + type = ParsedVectorFunction + expression_x = '0' + expression_y = '20' + expression_z = '0' + [] + [currentDensityFunction] + type = ParsedVectorFunction + expression_x = '0' + expression_y = '0' + expression_z = '1' + [] +[] + +[BCs] + [inlet] + type = VectorFunctionDirichletBC + variable = velocity + boundary = left + function = velocityInlet + [] + [no_slip] + type = VectorDirichletBC + variable = velocity + boundary = 'top bottom front back' + values = '0 0 0' + [] + [pressure_set] + type = PenaltyDirichletBC + variable = pressure + value = 0 + boundary = 'right' + penalty = 1e5 + [] +[] + +[Materials] + [const] + type = ADGenericConstantMaterial + prop_names = 'rho mu' + prop_values = '1 1' + [] + [ins_mat] + type = INSADTauMaterial + velocity = velocity + pressure = pressure + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = pressure + [] + + [momentum_convection] + type = INSADMomentumAdvection + variable = velocity + [] + + [momentum_viscous] + type = INSADMomentumViscous + variable = velocity + [] + + [momentum_pressure] + type = INSADMomentumPressure + variable = velocity + pressure = pressure + integrate_p_by_parts = true + [] + + [momentum_lorentz] + type = IMHDADMomentumLorentz + variable = velocity + currentDensity = currentDensity + magneticFieldFunction = magneticFieldFunction + [] +[] + +[AuxKernels] + [currentDensityKernel] + type = VectorFunctionAux + variable = currentDensity + function = currentDensityFunction + execute_on = INITIAL + [] +[] + +[Problem] + type = FEProblem +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + automatic_scaling = true + l_max_its = 1000 + nl_max_its = 1000 + petsc_options_iname = '-pc_type' + petsc_options_value = 'bjacobi' +[] + +[Outputs] + exodus = true +[] diff --git a/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_gradient.i b/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_gradient.i new file mode 100644 index 0000000..75563bc --- /dev/null +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_gradient.i @@ -0,0 +1,125 @@ +# Note: this case gives approximately the expected results, +# but there is a jump in electric potential from the boundary elements +# into the bulk. + +sigma = 1 + +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 50 + ny = 50 + xmax = 1 + ymax = 1 + xmin = -1 + ymin = -1 + elem_type = QUAD9 + [] +[] + +[Variables] + [currentDensity] + family = RAVIART_THOMAS + order = FIRST + [] + [electricPotential] + family = MONOMIAL + order = CONSTANT + [] +[] + +[Functions] + [potentialGradient] + type = ParsedFunction + expression = x + [] + [currentFlow] + type = ParsedVectorFunction + expression_x = -${sigma} + expression_y = 0 + div = 0 + [] +[] + +[Kernels] + [currentDensityConductivity] + type = IMHDADCurrentDensity + variable = currentDensity + [] + [electricField] + type = GradField + variable = currentDensity + coupled_scalar_variable = electricPotential + coeff = -1 + [] + + [divergenceFree] + type = DivField + variable = electricPotential + coupled_vector_variable = currentDensity + [] +[] + +[BCs] + [electricPotentialWall] + type = FunctionPenaltyDirichletBC + variable = electricPotential + function = potentialGradient + penalty = 1e5 + boundary = 'left right top bottom' + [] +[] + +[ICs] + [start] + type = FunctionIC + variable = electricPotential + function = potentialGradient + [] +[] + +[Postprocessors] + [L2Error] + type = ElementVectorL2Error + variable = currentDensity + function = currentFlow + [] + [HDivSemiError] + type = ElementHDivSemiError + variable = currentDensity + function = currentFlow + [] + [HDivError] + type = ElementHDivError + variable = currentDensity + function = currentFlow + [] +[] + +[Materials] + [const] + type = ADGenericConstantMaterial + prop_names = 'conductivity' + prop_values = '${sigma}' + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + petsc_options_iname = '-pc_type' + petsc_options_value = ' bjacobi' +[] + +[Outputs] + exodus = true + execute_on = 'NONLINEAR' +[] \ No newline at end of file diff --git a/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_solenoidal.i b/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_solenoidal.i new file mode 100644 index 0000000..41650d5 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_solenoidal.i @@ -0,0 +1,122 @@ +sigma = 1 + +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 50 + ny = 50 + xmax = 1 + ymax = 1 + xmin = -1 + ymin = -1 + elem_type = QUAD9 + [] +[] + +[Functions] + [solenoid] + type = ParsedVectorFunction + expression_x = y + expression_y = -x + div = 0 + [] +[] + +[Variables] + [currentDensity] + family = RAVIART_THOMAS + order = FIRST + [] + [electricPotential] + family = MONOMIAL + order = CONSTANT + [] +[] + +[Kernels] + [currentDensityConductivity] + type = IMHDADCurrentDensity + variable = currentDensity + [] + [electricField] + type = GradField + variable = currentDensity + coupled_scalar_variable = electricPotential + coeff = -1 + [] + [currentDensitySource] + type = VectorBodyForce + variable = currentDensity + function = solenoid + [] + + [divergenceFree] + type = DivField + variable = electricPotential + coupled_vector_variable = currentDensity + [] +[] + +[BCs] + [current_wall_value] + type = VectorPenaltyDirichletBC + variable = currentDensity + x_exact_sln = ${Functions/solenoid/expression_x} + y_exact_sln = ${Functions/solenoid/expression_y} + penalty = 1e3 + boundary = 'left right top bottom' + [] + [potential_wall_zero] # improves convergence but HDivSemiError increases? + type = PenaltyDirichletBC + variable = electricPotential + value = 0 + penalty = 1e3 + boundary = 'left right top bottom' + [] +[] + +[Postprocessors] + [L2Error] + type = ElementVectorL2Error + variable = currentDensity + function = solenoid + [] + [HDivSemiError] + type = ElementHDivSemiError + variable = currentDensity + function = solenoid + [] + [HDivError] + type = ElementHDivError + variable = currentDensity + function = solenoid + [] +[] + +[Materials] + [const] + type = ADGenericConstantMaterial + prop_names = 'conductivity' + prop_values = '${sigma}' + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + petsc_options_iname = '-pc_type' + petsc_options_value = ' bjacobi' +[] + +[Outputs] + exodus = true + execute_on = 'NONLINEAR' +[] \ No newline at end of file diff --git a/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic_solenoidal.i b/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic_solenoidal.i new file mode 100644 index 0000000..b57c230 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic_solenoidal.i @@ -0,0 +1,127 @@ +sigma = 1 + +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 3 + nx = 16 + ny = 16 + nz = 16 + xmax = 1 + ymax = 1 + zmax = 1 + xmin = -1 + ymin = -1 + zmin = -1 + elem_type = HEX27 + [] +[] + +[Functions] + [solenoid] + type = ParsedVectorFunction + expression_x = y + expression_y = -x + expression_z = 0 + div = 0 + [] +[] + +[Variables] + [currentDensity] + family = RAVIART_THOMAS + order = FIRST + [] + [electricPotential] + family = MONOMIAL + order = CONSTANT + [] +[] + +[Kernels] + [currentDensityConductivity] + type = IMHDADCurrentDensity + variable = currentDensity + [] + [electricField] + type = GradField + variable = currentDensity + coupled_scalar_variable = electricPotential + coeff = -1 + [] + [currentDensitySource] + type = VectorBodyForce + variable = currentDensity + function = solenoid + [] + + [divergenceFree] + type = DivField + variable = electricPotential + coupled_vector_variable = currentDensity + [] +[] + +[BCs] + [current_wall_value] + type = VectorPenaltyDirichletBC + variable = currentDensity + x_exact_sln = ${Functions/solenoid/expression_x} + y_exact_sln = ${Functions/solenoid/expression_y} + z_exact_sln = ${Functions/solenoid/expression_z} + penalty = 1e3 + boundary = 'front back left right top bottom' + [] + [potential_wall_zero] # improves convergence but HDivSemiError increases? + type = PenaltyDirichletBC + variable = electricPotential + value = 0 + penalty = 1e3 + boundary = 'front back left right top bottom' + [] +[] + +[Postprocessors] + [L2Error] + type = ElementVectorL2Error + variable = currentDensity + function = solenoid + [] + [HDivSemiError] + type = ElementHDivSemiError + variable = currentDensity + function = solenoid + [] + [HDivError] + type = ElementHDivError + variable = currentDensity + function = solenoid + [] +[] + +[Materials] + [const] + type = ADGenericConstantMaterial + prop_names = 'conductivity' + prop_values = '${sigma}' + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + petsc_options_iname = '-pc_type' + petsc_options_value = ' bjacobi' +[] + +[Outputs] + exodus = true + execute_on = 'NONLINEAR' +[] \ No newline at end of file diff --git a/examples/mhd/charge_conserving/development/step4_coupledAuxVelocity/3d_coupledAuxVelocity.i b/examples/mhd/charge_conserving/development/step4_coupledAuxVelocity/3d_coupledAuxVelocity.i new file mode 100644 index 0000000..b72aa3a --- /dev/null +++ b/examples/mhd/charge_conserving/development/step4_coupledAuxVelocity/3d_coupledAuxVelocity.i @@ -0,0 +1,267 @@ +sigma = 1 + +U_AVG = 1 + +N_X = 50 +N_Y_half = 10 +N_Z_half = 10 +ELEMENT_TYPE = HEX27 + +GRADING_R_Y = 5 +GRADING_R_Z = 5 +RATIO_Y_FWD = ${fparse GRADING_R_Y ^ (1/(N_Y_half - 1))} +RATIO_Y_INV = ${fparse 1/RATIO_Y_FWD} +RATIO_Z_FWD = ${fparse GRADING_R_Z ^ (1/(N_Z_half - 1))} +RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} + +[Mesh] + [meshTopBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopBack' + [] + [meshTopFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopFront' + [] + [meshBottomBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomBack' + [] + [meshBottomFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomFront' + [] + [meshTop] + type = StitchedMeshGenerator + inputs = 'meshTopBack meshTopFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTopBack_front meshTopFront_back' + [] + [renameMeshTop] + type = RenameBoundaryGenerator + input = meshTop + old_boundary = ' + meshTopBack_top meshTopFront_top + meshTopBack_right meshTopFront_right + meshTopBack_left meshTopFront_left + meshTopBack_bottom meshTopFront_bottom + ' + new_boundary = ' + top top + meshTop_right meshTop_right + meshTop_left meshTop_left + meshTop_bottom meshTop_bottom + ' + [] + [meshBottom] + type = StitchedMeshGenerator + inputs = 'meshBottomBack meshBottomFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshBottomBack_front meshBottomFront_back' + [] + [renameMeshBottom] + type = RenameBoundaryGenerator + input = meshBottom + old_boundary = ' + meshBottomBack_top meshBottomFront_top + meshBottomBack_right meshBottomFront_right + meshBottomBack_left meshBottomFront_left + meshBottomBack_bottom meshBottomFront_bottom + ' + new_boundary = ' + meshBottom_top meshBottom_top + meshBottom_right meshBottom_right + meshBottom_left meshBottom_left + bottom bottom + ' + [] + [mesh] + type = StitchedMeshGenerator + inputs = 'renameMeshTop renameMeshBottom' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTop_bottom meshBottom_top' + [] + [renameMesh] + type = RenameBoundaryGenerator + input = mesh + old_boundary = ' + meshBottomFront_front meshTopFront_front + meshBottomBack_back meshTopBack_back + meshBottom_left meshTop_left + meshBottom_right meshTop_right + ' + new_boundary = ' + front front + back back + left left + right right + ' + [] +[] + +[Functions] + [velocityFunction] + type = ParsedVectorFunction + symbol_names = 'y_max z_max' + symbol_values = '1 1' + expression_x = '(9/4) * ${U_AVG} * (1 - (y * y) / (y_max * y_max)) * (1 - (z * z) / (z_max * z_max))' + expression_y = '0' + expression_z = '0' + [] + [magneticFieldFunction] + type = ParsedVectorFunction + expression_x = '0' + expression_y = '20' + expression_z = '0' + [] +[] + +[Variables] + [currentDensity] + family = RAVIART_THOMAS + order = FIRST + [] + [electricPotential] + family = MONOMIAL + order = CONSTANT + [] +[] + +[AuxVariables] + [velocity] + family = LAGRANGE_VEC + order = FIRST + [] +[] + +[Kernels] + [currentDensityConductivity] + type = IMHDADCurrentDensity + variable = currentDensity + [] + [electricField] + type = GradField + variable = currentDensity + coupled_scalar_variable = electricPotential + coeff = -1 + [] + [currentDensityCoupling] + type = IMHDADCurrentUxB + variable = currentDensity + velocity = velocity + magneticFieldFunction = magneticFieldFunction + [] + + [divergenceFree] + type = DivField + variable = electricPotential + coupled_vector_variable = currentDensity + [] +[] + +[AuxKernels] + [velocityKernel] + type = VectorFunctionAux + variable = velocity + function = velocityFunction + execute_on = INITIAL + [] +[] + +[BCs] + [current_wall_normal] + # Trying to set J.n=0 on the walls. + # Not sure this is the correct BC to use; + # name suggests it is setting divJ=0 rather than J.n=0, + # but divJ=0 is enforced by the equations + # and having this BC gives a J field that looks right qualitatively + type = VectorDivPenaltyDirichletBC + variable = currentDensity + function_x = 0 + function_y = 0 + function_z = 0 + penalty = 1e7 + boundary = 'front back left right top bottom' + [] +[] + +[Materials] + [const] + type = ADGenericConstantMaterial + prop_names = 'conductivity' + prop_values = '${sigma}' + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + l_max_its = 1000 + nl_max_its = 1000 + petsc_options_iname = '-pc_type' + petsc_options_value = ' cholesky' +[] + +[Outputs] + exodus = true + execute_on = 'NONLINEAR' +[] \ No newline at end of file diff --git a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled_hunt.i b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled_hunt.i new file mode 100644 index 0000000..607bc69 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled_hunt.i @@ -0,0 +1,325 @@ +sigma = 1 + +U_AVG = 1 + +N_X = 50 +N_Y_half = 20 +N_Z_half = 10 +ELEMENT_TYPE = HEX27 + +GRADING_R_Y = 5 +GRADING_R_Z = 5 +RATIO_Y_FWD = ${fparse GRADING_R_Y ^ (1/(N_Y_half - 1))} +RATIO_Y_INV = ${fparse 1/RATIO_Y_FWD} +RATIO_Z_FWD = ${fparse GRADING_R_Z ^ (1/(N_Z_half - 1))} +RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} + +[Mesh] + [meshTopBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopBack' + [] + [meshTopFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopFront' + [] + [meshBottomBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomBack' + [] + [meshBottomFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomFront' + [] + [meshTop] + type = StitchedMeshGenerator + inputs = 'meshTopBack meshTopFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTopBack_front meshTopFront_back' + [] + [renameMeshTop] + type = RenameBoundaryGenerator + input = meshTop + old_boundary = ' + meshTopBack_top meshTopFront_top + meshTopBack_right meshTopFront_right + meshTopBack_left meshTopFront_left + meshTopBack_bottom meshTopFront_bottom + ' + new_boundary = ' + top top + meshTop_right meshTop_right + meshTop_left meshTop_left + meshTop_bottom meshTop_bottom + ' + [] + [meshBottom] + type = StitchedMeshGenerator + inputs = 'meshBottomBack meshBottomFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshBottomBack_front meshBottomFront_back' + [] + [renameMeshBottom] + type = RenameBoundaryGenerator + input = meshBottom + old_boundary = ' + meshBottomBack_top meshBottomFront_top + meshBottomBack_right meshBottomFront_right + meshBottomBack_left meshBottomFront_left + meshBottomBack_bottom meshBottomFront_bottom + ' + new_boundary = ' + meshBottom_top meshBottom_top + meshBottom_right meshBottom_right + meshBottom_left meshBottom_left + bottom bottom + ' + [] + [mesh] + type = StitchedMeshGenerator + inputs = 'renameMeshTop renameMeshBottom' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTop_bottom meshBottom_top' + [] + [renameMesh] + type = RenameBoundaryGenerator + input = mesh + old_boundary = ' + meshBottomFront_front meshTopFront_front + meshBottomBack_back meshTopBack_back + meshBottom_left meshTop_left + meshBottom_right meshTop_right + ' + new_boundary = ' + front front + back back + left left + right right + ' + [] +[] + +[Functions] + [velocityFunction] + type = ParsedVectorFunction + symbol_names = 'y_max z_max' + symbol_values = '1 1' + expression_x = '(9/4) * ${U_AVG} * (1 - (y * y) / (y_max * y_max)) * (1 - (z * z) / (z_max * z_max))' + expression_y = '0' + expression_z = '0' + [] + [magneticFieldFunction] + type = ParsedVectorFunction + expression_x = '0' + expression_y = '20' + expression_z = '0' + [] +[] + +[Variables] + [velocity] + family = LAGRANGE_VEC + order = SECOND + [] + [pressure] + family = LAGRANGE + order = FIRST + [] + [currentDensity] + family = RAVIART_THOMAS + order = FIRST + [] + [electricPotential] + family = MONOMIAL + order = CONSTANT + [] +[] + +[ICs] + [velocityIC] + # type = VectorConstantIC + # x_value = ${U_AVG} + # y_value = 1e-15 + # z_value = 1e-15 + type = VectorFunctionIC + function = velocityFunction + variable = velocity + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = pressure + [] + + [momentum_convection] + type = INSADMomentumAdvection + variable = velocity + [] + [momentum_viscous] + type = INSADMomentumViscous + variable = velocity + [] + [momentum_pressure] + type = INSADMomentumPressure + variable = velocity + pressure = pressure + integrate_p_by_parts = true + [] + [momentum_lorentz] + type = IMHDADMomentumLorentz + variable = velocity + currentDensity = currentDensity + magneticFieldFunction = magneticFieldFunction + [] + + [currentDensityConductivity] + type = IMHDADCurrentDensity + variable = currentDensity + [] + [electricField] + type = GradField + variable = currentDensity + coupled_scalar_variable = electricPotential + coeff = -1 + [] + [currentDensityCoupling] + type = IMHDADCurrentUxB + variable = currentDensity + velocity = velocity + magneticFieldFunction = magneticFieldFunction + [] + + [divergenceFree] + type = DivField + variable = electricPotential + coupled_vector_variable = currentDensity + [] +[] + +[BCs] + [inlet] + type = VectorFunctionDirichletBC + variable = velocity + boundary = left + function = velocityFunction + [] + [no_slip] + type = VectorDirichletBC + variable = velocity + boundary = 'top bottom front back' + values = '0 0 0' + [] + [pressure_set] + type = PenaltyDirichletBC + variable = pressure + value = 0 + boundary = 'right' + penalty = 1e5 + [] + [current_insulating] + # Setting J.n=0 on the insulating walls, inlet and outlet. + type = VectorDivPenaltyDirichletBC + variable = currentDensity + function_x = 0 + function_y = 0 + function_z = 0 + penalty = 1e7 + boundary = 'front back left right' + [] + [potential_conducting] + type = PenaltyDirichletBC + variable = electricPotential + value = 0 + boundary = 'top bottom' + penalty = 1e5 + [] +[] + +[Materials] + [const] + type = ADGenericConstantMaterial + prop_names = 'rho mu conductivity' + prop_values = '1 1 ${sigma}' + [] + [ins_mat] + type = INSADTauMaterial + velocity = velocity + pressure = pressure + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + petsc_options_iname = '-pc_type' + petsc_options_value = ' lu' + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + automatic_scaling = false + l_max_its = 1000 + nl_max_its = 1000 +[] + +[Outputs] + exodus = true + execute_on = 'NONLINEAR' +[] diff --git a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled_shercliff.i b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled_shercliff.i new file mode 100644 index 0000000..c5cacc2 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled_shercliff.i @@ -0,0 +1,331 @@ +sigma = 1 + +U_AVG = 1 + +N_X = 50 +N_Y_half = 20 +N_Z_half = 10 +ELEMENT_TYPE = HEX27 + +GRADING_R_Y = 5 +GRADING_R_Z = 5 +RATIO_Y_FWD = ${fparse GRADING_R_Y ^ (1/(N_Y_half - 1))} +RATIO_Y_INV = ${fparse 1/RATIO_Y_FWD} +RATIO_Z_FWD = ${fparse GRADING_R_Z ^ (1/(N_Z_half - 1))} +RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} + +[Mesh] + [meshTopBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopBack' + [] + [meshTopFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = 0 + ymax = 1 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_INV} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshTopFront' + [] + [meshBottomBack] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = -1 + zmax = 0 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_FWD} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomBack' + [] + [meshBottomFront] + type = GeneratedMeshGenerator + dim = 3 + nx = ${N_X} + ny = ${N_Y_half} + nz = ${N_Z_half} + xmin = 0 + xmax = 20 + ymin = -1 + ymax = 0 + zmin = 0 + zmax = 1 + bias_y = ${RATIO_Y_FWD} + bias_z = ${RATIO_Z_INV} + elem_type = ${ELEMENT_TYPE} + boundary_name_prefix = 'meshBottomFront' + [] + [meshTop] + type = StitchedMeshGenerator + inputs = 'meshTopBack meshTopFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTopBack_front meshTopFront_back' + [] + [renameMeshTop] + type = RenameBoundaryGenerator + input = meshTop + old_boundary = ' + meshTopBack_top meshTopFront_top + meshTopBack_right meshTopFront_right + meshTopBack_left meshTopFront_left + meshTopBack_bottom meshTopFront_bottom + ' + new_boundary = ' + top top + meshTop_right meshTop_right + meshTop_left meshTop_left + meshTop_bottom meshTop_bottom + ' + [] + [meshBottom] + type = StitchedMeshGenerator + inputs = 'meshBottomBack meshBottomFront' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshBottomBack_front meshBottomFront_back' + [] + [renameMeshBottom] + type = RenameBoundaryGenerator + input = meshBottom + old_boundary = ' + meshBottomBack_top meshBottomFront_top + meshBottomBack_right meshBottomFront_right + meshBottomBack_left meshBottomFront_left + meshBottomBack_bottom meshBottomFront_bottom + ' + new_boundary = ' + meshBottom_top meshBottom_top + meshBottom_right meshBottom_right + meshBottom_left meshBottom_left + bottom bottom + ' + [] + [mesh] + type = StitchedMeshGenerator + inputs = 'renameMeshTop renameMeshBottom' + clear_stitched_boundary_ids = true + stitch_boundaries_pairs = 'meshTop_bottom meshBottom_top' + [] + [renameMesh] + type = RenameBoundaryGenerator + input = mesh + old_boundary = ' + meshBottomFront_front meshTopFront_front + meshBottomBack_back meshTopBack_back + meshBottom_left meshTop_left + meshBottom_right meshTop_right + ' + new_boundary = ' + front front + back back + left left + right right + ' + [] +[] + +[Functions] + [velocityFunction] + type = ParsedVectorFunction + symbol_names = 'y_max z_max' + symbol_values = '1 1' + expression_x = '(9/4) * ${U_AVG} * (1 - (y * y) / (y_max * y_max)) * (1 - (z * z) / (z_max * z_max))' + expression_y = '0' + expression_z = '0' + [] + [magneticFieldFunction] + type = ParsedVectorFunction + expression_x = '0' + expression_y = '20' + expression_z = '0' + [] + [divergenceFreeCheck] + type = ParsedVectorFunction + div = '0' + [] +[] + +[Variables] + [velocity] + family = LAGRANGE_VEC + order = SECOND + [] + [pressure] + family = LAGRANGE + order = FIRST + [] + [currentDensity] + family = RAVIART_THOMAS + order = FIRST + [] + [electricPotential] + family = MONOMIAL + order = CONSTANT + [] +[] + +[ICs] + [velocityIC] + # type = VectorConstantIC + # x_value = ${U_AVG} + # y_value = 1e-15 + # z_value = 1e-15 + type = VectorFunctionIC + function = velocityFunction + variable = velocity + [] +[] + +[Kernels] + [mass] + type = INSADMass + variable = pressure + [] + + [momentum_convection] + type = INSADMomentumAdvection + variable = velocity + [] + [momentum_viscous] + type = INSADMomentumViscous + variable = velocity + [] + [momentum_pressure] + type = INSADMomentumPressure + variable = velocity + pressure = pressure + integrate_p_by_parts = true + [] + [momentum_lorentz] + type = IMHDADMomentumLorentz + variable = velocity + currentDensity = currentDensity + magneticFieldFunction = magneticFieldFunction + [] + + [currentDensityConductivity] + type = IMHDADCurrentDensity + variable = currentDensity + [] + [electricField] + type = GradField + variable = currentDensity + coupled_scalar_variable = electricPotential + coeff = -1 + [] + [currentDensityCoupling] + type = IMHDADCurrentUxB + variable = currentDensity + velocity = velocity + magneticFieldFunction = magneticFieldFunction + [] + + [divergenceFree] + type = DivField + variable = electricPotential + coupled_vector_variable = currentDensity + [] +[] + +[BCs] + [inlet] + type = VectorFunctionDirichletBC + variable = velocity + boundary = left + function = velocityFunction + [] + [no_slip] + type = VectorDirichletBC + variable = velocity + boundary = 'top bottom front back' + values = '0 0 0' + [] + [pressure_set] + type = PenaltyDirichletBC + variable = pressure + value = 0 + boundary = 'right' + penalty = 1e5 + [] + [current_wall_normal] + # Setting J.n=0 on the walls, inlet and outlet. + type = VectorDivPenaltyDirichletBC + variable = currentDensity + function_x = 0 + function_y = 0 + function_z = 0 + penalty = 1e7 + boundary = 'front back left right top bottom' + [] +[] + +[Materials] + [const] + type = ADGenericConstantMaterial + prop_names = 'rho mu conductivity' + prop_values = '1 1 ${sigma}' + [] + [ins_mat] + type = INSADTauMaterial + velocity = velocity + pressure = pressure + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + petsc_options_iname = '-pc_type' + petsc_options_value = ' lu' + [] +[] + +[Postprocessors] + [currentDivergenceFreeError] + type = ElementHDivSemiError + variable = currentDensity + function = divergenceFreeCheck + execute_on = 'NONLINEAR' + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + automatic_scaling = false + l_max_its = 1000 + nl_max_its = 1000 +[] + +[Outputs] + exodus = true + execute_on = 'NONLINEAR' +[] diff --git a/include/kernels/IMHDADCurrentDensity.h b/include/kernels/IMHDADCurrentDensity.h new file mode 100644 index 0000000..f408987 --- /dev/null +++ b/include/kernels/IMHDADCurrentDensity.h @@ -0,0 +1,20 @@ +#pragma once + +#include "ADKernelValue.h" + +/** + * This class computes the residual and Jacobian contributions for the + * left-hand-side of the explicit calculation of the electric current density equation. + */ +class IMHDADCurrentDensity : public ADVectorKernelValue +{ +public: + static InputParameters validParams(); + + IMHDADCurrentDensity(const InputParameters & parameters); + +protected: + virtual ADRealVectorValue precomputeQpResidual() override; + + const ADMaterialProperty & _conductivity; +}; diff --git a/include/kernels/IMHDADCurrentUxB.h b/include/kernels/IMHDADCurrentUxB.h new file mode 100644 index 0000000..e6bd4b9 --- /dev/null +++ b/include/kernels/IMHDADCurrentUxB.h @@ -0,0 +1,32 @@ +#pragma once + +#include "ADKernelValue.h" +#include "Function.h" + +class Function; + +/** + * This class computes the residual and Jacobian contributions for the + * U cross B term of the electric current density equation + * using a coupled velocity variable and an imposed magnetic field function. + */ +class IMHDADCurrentUxB : public ADVectorKernelValue +{ +public: + static InputParameters validParams(); + + IMHDADCurrentUxB(const InputParameters & parameters); + +protected: + virtual ADRealVectorValue precomputeQpResidual() override; + + const ADVectorVariableValue & _velocity; + + /// Optional vectorValue function + const Function * const _magnetic_field; + + /// Optional component function value + const Function & _magnetic_field_x; + const Function & _magnetic_field_y; + const Function & _magnetic_field_z; +}; diff --git a/include/kernels/IMHDADMomentumLorentz.h b/include/kernels/IMHDADMomentumLorentz.h new file mode 100644 index 0000000..a04a7fe --- /dev/null +++ b/include/kernels/IMHDADMomentumLorentz.h @@ -0,0 +1,32 @@ +#pragma once + +#include "ADKernelValue.h" +#include "Function.h" + +class Function; + +/** + * This class computes the residual and Jacobian contributions for the + * Lorentz force term of the incompressible inductionless MHD momentum equation + * using a coupled current density variable and an imposed magnetic field function. + */ +class IMHDADMomentumLorentz : public ADVectorKernelValue +{ +public: + static InputParameters validParams(); + + IMHDADMomentumLorentz(const InputParameters & parameters); + +protected: + virtual ADRealVectorValue precomputeQpResidual() override; + + const ADVectorVariableValue & _current_density; + + /// Optional vectorValue function + const Function * const _magnetic_field; + + /// Optional component function value + const Function & _magnetic_field_x; + const Function & _magnetic_field_y; + const Function & _magnetic_field_z; +}; diff --git a/scripts/NSFVBase.patch b/scripts/NSFVBase.patch deleted file mode 100644 index 3c47c88..0000000 --- a/scripts/NSFVBase.patch +++ /dev/null @@ -1,22 +0,0 @@ -diff --git a/modules/navier_stokes/include/base/NSFVBase.h b/modules/navier_stokes/include/base/NSFVBase.h -index e90b4a85..7a9ac856 100644 ---- a/modules/navier_stokes/include/base/NSFVBase.h -+++ b/modules/navier_stokes/include/base/NSFVBase.h -@@ -3537,12 +3537,13 @@ NSFVBase::checkBlockwiseConsistency(const std::string block_param_name - "Block '" + block + - "' is not present in the block restriction of the fluid flow action!"); - -- for (const auto & param_name : parameter_names) -+ for (unsigned int param_i = 0; param_i < parameter_names.size(); ++param_i) - { -- const std::vector & param_vector = parameters().template get>(param_name); -+ const std::vector & param_vector = -+ parameters().template get>(parameter_names[param_i]); - if (block_names.size() != param_vector.size()) -- paramError(param_name, -- "The number of entries in '" + param_name + "' (" + -+ paramError(parameter_names[param_i], -+ "The number of entries in '" + parameter_names[param_i] + "' (" + - std::to_string(param_vector.size()) + - ") is not the same as the number of blocks" - " (" + diff --git a/scripts/csd3_python_prep.sh b/scripts/csd3_python_prep.sh index 675af34..e4aed63 100755 --- a/scripts/csd3_python_prep.sh +++ b/scripts/csd3_python_prep.sh @@ -6,3 +6,4 @@ module load python/3.11.0-icl pip3.11 install --user --upgrade packaging pip3.11 install --user --upgrade pyyaml +pip3.11 install --user --upgrade jinja2 diff --git a/scripts/install_csd3_cclake.sh b/scripts/install_csd3_cclake.sh index c3adc8b..70c092d 100755 --- a/scripts/install_csd3_cclake.sh +++ b/scripts/install_csd3_cclake.sh @@ -60,7 +60,7 @@ METHODS="opt" ./scripts/update_and_rebuild_libmesh.sh --with-mpi # 8 for each first order variable # 27 for each second order variable -./configure --with-derivative-size=81 +./configure --with-derivative-size=89 cd $PROTEUS_DIR make -j $MOOSE_JOBS diff --git a/scripts/install_csd3_sapphire.sh b/scripts/install_csd3_sapphire.sh index 78ef597..de963b4 100755 --- a/scripts/install_csd3_sapphire.sh +++ b/scripts/install_csd3_sapphire.sh @@ -41,7 +41,7 @@ unset PETSC_DIR PETSC_ARCH # Build libMesh -./scripts/update_and_rebuild_libmesh.sh --with-mpi +./scripts/update_and_rebuild_libmesh.sh --with-mpi --disable-netgen # Build WASP @@ -52,11 +52,7 @@ unset PETSC_DIR PETSC_ARCH # 8 for each first order variable # 27 for each second order variable -./configure --with-derivative-size=81 - -# Apply NSFV patch - -git apply $PROTEUS_DIR/scripts/NSFVBase.patch +./configure --with-derivative-size=89 cd $PROTEUS_DIR make -j $MOOSE_JOBS diff --git a/scripts/install_ubuntu.sh b/scripts/install_ubuntu.sh index cb81fcc..0dadab8 100755 --- a/scripts/install_ubuntu.sh +++ b/scripts/install_ubuntu.sh @@ -61,7 +61,7 @@ unset PETSC_DIR PETSC_ARCH # 8 for each first order variable # 27 for each second order variable -./configure --with-derivative-size=81 +./configure --with-derivative-size=89 cd $PROTEUS_DIR make -j $MOOSE_JOBS diff --git a/src/kernels/IMHDADCurrentDensity.C b/src/kernels/IMHDADCurrentDensity.C new file mode 100644 index 0000000..b0b2fbc --- /dev/null +++ b/src/kernels/IMHDADCurrentDensity.C @@ -0,0 +1,28 @@ +#include "IMHDADCurrentDensity.h" + +registerMooseObject("ProteusApp", IMHDADCurrentDensity); + +InputParameters +IMHDADCurrentDensity::validParams() +{ + InputParameters params = ADVectorKernelValue::validParams(); + params.addClassDescription( + "The electric current term ($\\sigma^{-1}\\vec{J}$), " + "with the weak form of ($\\phi_u, \\sigma^{-1}\\vec{J}$). " + "The Jacobian is computed using automatic differentiation."); + params.addParam("conductivity", "conductivity", "The name of the conductivity"); + + return params; +} + +IMHDADCurrentDensity::IMHDADCurrentDensity(const InputParameters & parameters) + : ADVectorKernelValue(parameters), + _conductivity(getADMaterialProperty("conductivity")) +{ +} + +ADRealVectorValue +IMHDADCurrentDensity::precomputeQpResidual() +{ + return _u[_qp]/_conductivity[_qp]; +} diff --git a/src/kernels/IMHDADCurrentUxB.C b/src/kernels/IMHDADCurrentUxB.C new file mode 100644 index 0000000..253cd43 --- /dev/null +++ b/src/kernels/IMHDADCurrentUxB.C @@ -0,0 +1,51 @@ +#include "IMHDADCurrentUxB.h" + +registerMooseObject("ProteusApp", IMHDADCurrentUxB); + +InputParameters +IMHDADCurrentUxB::validParams() +{ + InputParameters params = ADVectorKernelValue::validParams(); + params.addClassDescription( + "The electric current velocity source term ($-\\vec{U} \\times \\vec{B}$), " + "with the weak form of ($\\phi_u, -\\vec{U} \\times \\vec{B}$). " + "The Jacobian is computed using automatic differentiation."); + params.addRequiredCoupledVar("velocity", "The variable representing the velocity"); + params.addParam("magneticFieldFunction", + "A function defining a vectorValue method that describes the imposed magnetic field"); + params.addParam( + "magneticFieldFunction_x", "0", "A function that describes the x-component of the imposed magnetic field"); + params.addParam( + "magneticFieldFunction_y", "0", "A function that describes the y-component of the imposed magnetic field"); + params.addParam( + "magneticFieldFunction_z", "0", "A function that describes the z-component of the imposed magnetic field"); + + return params; +} + +IMHDADCurrentUxB::IMHDADCurrentUxB(const InputParameters & parameters) + : ADVectorKernelValue(parameters), + _velocity(adCoupledVectorValue("velocity")), + _magnetic_field(isParamValid("magneticFieldFunction") ? &getFunction("magneticFieldFunction") : nullptr), + _magnetic_field_x(getFunction("magneticFieldFunction_x")), + _magnetic_field_y(getFunction("magneticFieldFunction_y")), + _magnetic_field_z(getFunction("magneticFieldFunction_z")) +{ + if (_magnetic_field && parameters.isParamSetByUser("_magnetic_field_x")) + paramError("magneticFieldFunction_x", "The 'magneticFieldFunction' and 'magneticFieldFunction_x' parameters cannot both be set."); + if (_magnetic_field && parameters.isParamSetByUser("_magnetic_field_y")) + paramError("magneticFieldFunction_y", "The 'magneticFieldFunction' and 'magneticFieldFunction_y' parameters cannot both be set."); + if (_magnetic_field && parameters.isParamSetByUser("_magnetic_field_z")) + paramError("magneticFieldFunction_z", "The 'magneticFieldFunction' and 'magneticFieldFunction_z' parameters cannot both be set."); +} + +ADRealVectorValue +IMHDADCurrentUxB::precomputeQpResidual() +{ + if (_magnetic_field) + return -_velocity[_qp].cross(_magnetic_field->vectorValue(_t, _q_point[_qp])); + else + return -_velocity[_qp].cross(RealVectorValue(_magnetic_field_x.value(_t, _q_point[_qp]), + _magnetic_field_y.value(_t, _q_point[_qp]), + _magnetic_field_z.value(_t, _q_point[_qp]))); +} diff --git a/src/kernels/IMHDADMomentumLorentz.C b/src/kernels/IMHDADMomentumLorentz.C new file mode 100644 index 0000000..d4e46f2 --- /dev/null +++ b/src/kernels/IMHDADMomentumLorentz.C @@ -0,0 +1,51 @@ +#include "IMHDADMomentumLorentz.h" + +registerMooseObject("ProteusApp", IMHDADMomentumLorentz); + +InputParameters +IMHDADMomentumLorentz::validParams() +{ + InputParameters params = ADVectorKernelValue::validParams(); + params.addClassDescription( + "The Lorentz force term ($\\vec{J} \\times \\vec{B}_0$), " + "with the weak form of ($\\phi_u, \\vec{J} \\times \\vec{B}_0$)." + "The Jacobian is computed using automatic differentiation."); + params.addRequiredCoupledVar("currentDensity", "The variable representing the electric current density"); + params.addParam("magneticFieldFunction", + "A function defining a vectorValue method that describes the imposed magnetic field"); + params.addParam( + "magneticFieldFunction_x", "0", "A function that describes the x-component of the imposed magnetic field"); + params.addParam( + "magneticFieldFunction_y", "0", "A function that describes the y-component of the imposed magnetic field"); + params.addParam( + "magneticFieldFunction_z", "0", "A function that describes the z-component of the imposed magnetic field"); + + return params; +} + +IMHDADMomentumLorentz::IMHDADMomentumLorentz(const InputParameters & parameters) + : ADVectorKernelValue(parameters), + _current_density(adCoupledVectorValue("currentDensity")), + _magnetic_field(isParamValid("magneticFieldFunction") ? &getFunction("magneticFieldFunction") : nullptr), + _magnetic_field_x(getFunction("magneticFieldFunction_x")), + _magnetic_field_y(getFunction("magneticFieldFunction_y")), + _magnetic_field_z(getFunction("magneticFieldFunction_z")) +{ + if (_magnetic_field && parameters.isParamSetByUser("_magnetic_field_x")) + paramError("magneticFieldFunction_x", "The 'magneticFieldFunction' and 'magneticFieldFunction_x' parameters cannot both be set."); + if (_magnetic_field && parameters.isParamSetByUser("_magnetic_field_y")) + paramError("magneticFieldFunction_y", "The 'magneticFieldFunction' and 'magneticFieldFunction_y' parameters cannot both be set."); + if (_magnetic_field && parameters.isParamSetByUser("_magnetic_field_z")) + paramError("magneticFieldFunction_z", "The 'magneticFieldFunction' and 'magneticFieldFunction_z' parameters cannot both be set."); +} + +ADRealVectorValue +IMHDADMomentumLorentz::precomputeQpResidual() +{ + if (_magnetic_field) + return -_current_density[_qp].cross(_magnetic_field->vectorValue(_t, _q_point[_qp])); + else + return -_current_density[_qp].cross(RealVectorValue(_magnetic_field_x.value(_t, _q_point[_qp]), + _magnetic_field_y.value(_t, _q_point[_qp]), + _magnetic_field_z.value(_t, _q_point[_qp]))); +}