From 05b4c3491c891a2015b449ae7246e5f08d5a26e4 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Tue, 14 May 2024 09:53:57 +0100 Subject: [PATCH 01/43] Add 2D Navier-Stokes 2nd order lagrange velocity 1st order monomial pressure example --- .../development/step1_2d_NS_elementTypes.i | 166 ++++++++++++++++++ 1 file changed, 166 insertions(+) create mode 100644 examples/mhd/charge_conserving/development/step1_2d_NS_elementTypes.i diff --git a/examples/mhd/charge_conserving/development/step1_2d_NS_elementTypes.i b/examples/mhd/charge_conserving/development/step1_2d_NS_elementTypes.i new file mode 100644 index 0000000..6950fcf --- /dev/null +++ b/examples/mhd/charge_conserving/development/step1_2d_NS_elementTypes.i @@ -0,0 +1,166 @@ +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 = 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 +[] From eb28f796b03ca1a06fef814eac55f3ceee60978d Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Tue, 14 May 2024 11:52:28 +0100 Subject: [PATCH 02/43] Add penalty Dirichlet BC --- .../development/step1_2d_NS_elementTypes.i | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/examples/mhd/charge_conserving/development/step1_2d_NS_elementTypes.i b/examples/mhd/charge_conserving/development/step1_2d_NS_elementTypes.i index 6950fcf..9aae768 100644 --- a/examples/mhd/charge_conserving/development/step1_2d_NS_elementTypes.i +++ b/examples/mhd/charge_conserving/development/step1_2d_NS_elementTypes.i @@ -95,7 +95,14 @@ element_type = 'QUAD9' 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 From 7c744bdfb718c7fa7e6d5684e21086c1fa73dfd9 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Tue, 14 May 2024 11:52:51 +0100 Subject: [PATCH 03/43] Add initial 3D NS element type test case --- .../development/step1_3d_NS_elementTypes.i | 354 ++++++++++++++++++ 1 file changed, 354 insertions(+) create mode 100644 examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i diff --git a/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i b/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i new file mode 100644 index 0000000..597d99f --- /dev/null +++ b/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i @@ -0,0 +1,354 @@ +U_AVG = 1 + +N_X = 10 +N_Y_half = 4 +N_Z_half = 4 +ELEMENT_TYPE = HEX20 + +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 + ' + [] +[] + +# 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 +# ' +# [] +# [] + +# # N_X = 200 +# # N_Y_half = 10 +# U_AVG = 1 +# element_type = 'QUAD9' + +# [Mesh] +# [mesh] +# type = GeneratedMeshGenerator +# dim = 2 +# nx = 200 +# ny = 20 +# # nz = 20 +# xmin = 0 +# xmax = 20 +# ymin = -1 +# ymax = 1 +# # zmin = -1 +# # zmax = 1 +# elem_type = ${element_type} +# [] +# [] + +[Variables] + [velocity] + family = LAGRANGE_VEC + order = SECOND + [] + [pressure] + family = L2_LAGRANGE + 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' + [] + # [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' + 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 +[] From 1528f98d08d6cd283b33098554f17bae4f909134 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Tue, 14 May 2024 16:17:03 +0100 Subject: [PATCH 04/43] Fix missing no-slip BCs, change element types, change executioner options, add ICs --- .../development/step1_3d_NS_elementTypes.i | 133 ++++-------------- 1 file changed, 29 insertions(+), 104 deletions(-) diff --git a/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i b/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i index 597d99f..aeb1155 100644 --- a/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i +++ b/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i @@ -1,9 +1,9 @@ U_AVG = 1 -N_X = 10 -N_Y_half = 4 -N_Z_half = 4 -ELEMENT_TYPE = HEX20 +N_X = 200 +N_Y_half = 10 +N_Z_half = 10 +ELEMENT_TYPE = HEX27 GRADING_R_Y = 5 GRADING_R_Z = 5 @@ -149,119 +149,43 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [] [] -# 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 -# ' -# [] -# [] - -# # N_X = 200 -# # N_Y_half = 10 -# U_AVG = 1 -# element_type = 'QUAD9' - -# [Mesh] -# [mesh] -# type = GeneratedMeshGenerator -# dim = 2 -# nx = 200 -# ny = 20 -# # nz = 20 -# xmin = 0 -# xmax = 20 -# ymin = -1 -# ymax = 1 -# # zmin = -1 -# # zmax = 1 -# elem_type = ${element_type} -# [] -# [] - [Variables] [velocity] family = LAGRANGE_VEC order = SECOND [] [pressure] - family = L2_LAGRANGE + family = MONOMIAL order = FIRST [] [] -# [ICs] -# [velocity] -# type = VectorFunctionIC -# variable = velocity -# function = velocityInlet -# [] -# [] +[ICs] + [velocity] + type = VectorFunctionIC + variable = velocity + function = velocityInlet + [] + [pressure] + type = FunctionIC + variable = pressure + function = pressureGradient + [] +[] [Functions] [velocityInlet] type = ParsedVectorFunction - symbol_names = 'y_max' - symbol_values = '1' - expression_x = '(3/2) * ${U_AVG} * (1 - (y * y) / (y_max * y_max))' + 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' + [] + [pressureGradient] + type = ParsedFunction + expression = '-0.01*x + 0.2' [] - # [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] @@ -274,7 +198,7 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [no_slip] type = VectorDirichletBC variable = velocity - boundary = 'top bottom' + boundary = 'top bottom front back' values = '0 0 0' [] [pressure_set] @@ -343,12 +267,13 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [Executioner] type = Steady solve_type = NEWTON - l_max_its = 100 + l_max_its = 1000 nl_max_its = 1000 petsc_options_iname = '-pc_type' - petsc_options_value = 'ilu' + petsc_options_value = 'bjacobi' [] [Outputs] exodus = true + execute_on = 'nonlinear' [] From 9ad4044ea8f72e7b774a8dd7d8066f7c5b3cbc1e Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Tue, 14 May 2024 16:17:40 +0100 Subject: [PATCH 05/43] Initial Lorentz force kernel implementation --- .../development/step2_2d_NS_LorentzForce.i | 194 ++++++++++++++++++ include/kernels/IMHDADMomentumLorentz.h | 22 ++ src/kernels/IMHDADMomentumLorentz.C | 30 +++ 3 files changed, 246 insertions(+) create mode 100644 examples/mhd/charge_conserving/development/step2_2d_NS_LorentzForce.i create mode 100644 include/kernels/IMHDADMomentumLorentz.h create mode 100644 src/kernels/IMHDADMomentumLorentz.C diff --git a/examples/mhd/charge_conserving/development/step2_2d_NS_LorentzForce.i b/examples/mhd/charge_conserving/development/step2_2d_NS_LorentzForce.i new file mode 100644 index 0000000..20005f1 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step2_2d_NS_LorentzForce.i @@ -0,0 +1,194 @@ +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 + [] +[] + +[AuxVariables] + [magneticField] + family = LAGRANGE_VEC + 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' + [] + [magneticFieldFunction] + type = ParsedVectorFunction + expression_x = '0' + expression_y = '20' + [] +[] + +[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 + [] +[] + +[AuxKernels] + [magneticFieldKernel] + type = VectorFunctionAux + variable = magneticField + function = magneticFieldFunction + execute_on = INITIAL + [] +[] + +[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/include/kernels/IMHDADMomentumLorentz.h b/include/kernels/IMHDADMomentumLorentz.h new file mode 100644 index 0000000..0b2eff5 --- /dev/null +++ b/include/kernels/IMHDADMomentumLorentz.h @@ -0,0 +1,22 @@ +#pragma once + +#include "ADKernelValue.h" + +/** + * This class computes the residual and Jacobian contributions for the + * Lorentz force term of the incompressible inductionless MHD momentum equation. + */ +class IMHDADMomentumLorentz : public ADVectorKernelValue +{ +public: + static InputParameters validParams(); + + IMHDADMomentumLorentz(const InputParameters & parameters); + +protected: + virtual ADRealVectorValue precomputeQpResidual() override; + + const ADVectorVariableValue & _current_density; + + const ADVectorVariableValue & _magnetic_field; +}; diff --git a/src/kernels/IMHDADMomentumLorentz.C b/src/kernels/IMHDADMomentumLorentz.C new file mode 100644 index 0000000..b6c59ff --- /dev/null +++ b/src/kernels/IMHDADMomentumLorentz.C @@ -0,0 +1,30 @@ +#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.addRequiredCoupledVar("magneticField", "The variable representing the magnetic field"); + + return params; +} + +IMHDADMomentumLorentz::IMHDADMomentumLorentz(const InputParameters & parameters) + : ADVectorKernelValue(parameters), + _current_density(adCoupledVectorValue("currentDensity")), + _magnetic_field(adCoupledVectorValue("magneticField")) +{ +} + +ADRealVectorValue +IMHDADMomentumLorentz::precomputeQpResidual() +{ + return _current_density[_qp].cross(_magnetic_field[_qp]); +} From a4c4fff09a1188053709c1d4b5d86e03bf2c2197 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Tue, 14 May 2024 16:45:15 +0100 Subject: [PATCH 06/43] Add current and Lorentz force kernel --- .../development/step2_2d_NS_LorentzForce.i | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/examples/mhd/charge_conserving/development/step2_2d_NS_LorentzForce.i b/examples/mhd/charge_conserving/development/step2_2d_NS_LorentzForce.i index 20005f1..331c329 100644 --- a/examples/mhd/charge_conserving/development/step2_2d_NS_LorentzForce.i +++ b/examples/mhd/charge_conserving/development/step2_2d_NS_LorentzForce.i @@ -70,6 +70,10 @@ element_type = 'QUAD9' family = LAGRANGE_VEC order = FIRST [] + [currentDensity] + family = LAGRANGE_VEC + order = FIRST + [] [] [ICs] @@ -93,6 +97,12 @@ element_type = 'QUAD9' expression_x = '0' expression_y = '20' [] + [currentDensityFunction] + type = ParsedVectorFunction + expression_x = '0' + expression_y = '0' + expression_z = '1' + [] [] [BCs] @@ -158,6 +168,13 @@ element_type = 'QUAD9' pressure = pressure integrate_p_by_parts = true [] + + [momentum_lorentz] + type = IMHDADMomentumLorentz + variable = velocity + currentDensity = currentDensity + magneticField = magneticField + [] [] [AuxKernels] @@ -167,6 +184,12 @@ element_type = 'QUAD9' function = magneticFieldFunction execute_on = INITIAL [] + [currentDensityKernel] + type = VectorFunctionAux + variable = currentDensity + function = currentDensityFunction + execute_on = INITIAL + [] [] [Problem] From 632805736f71e8599bff5079befcb26e76f2f509 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 15 May 2024 16:17:36 +0100 Subject: [PATCH 07/43] Add automatic scaling, remove function ICs, reduce x resolution --- .../development/step1_3d_NS_elementTypes.i | 26 +++++-------------- 1 file changed, 7 insertions(+), 19 deletions(-) diff --git a/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i b/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i index aeb1155..2c4f831 100644 --- a/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i +++ b/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i @@ -1,6 +1,6 @@ U_AVG = 1 -N_X = 200 +N_X = 50 N_Y_half = 10 N_Z_half = 10 ELEMENT_TYPE = HEX27 @@ -161,15 +161,12 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [] [ICs] - [velocity] - type = VectorFunctionIC + [velocityIC] + type = VectorConstantIC + x_value = ${U_AVG} + y_value = 1e-15 + z_value = 1e-15 variable = velocity - function = velocityInlet - [] - [pressure] - type = FunctionIC - variable = pressure - function = pressureGradient [] [] @@ -182,10 +179,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} expression_y = '0' expression_z = '0' [] - [pressureGradient] - type = ParsedFunction - expression = '-0.01*x + 0.2' - [] [] [BCs] @@ -208,12 +201,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} boundary = 'right' penalty = 1e5 [] - # [pressure_set] - # type = DirichletBC - # variable = pressure - # boundary = 'right' - # value = 0 - # [] [] [Materials] @@ -267,6 +254,7 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [Executioner] type = Steady solve_type = NEWTON + automatic_scaling = true l_max_its = 1000 nl_max_its = 1000 petsc_options_iname = '-pc_type' From e6a1ced8225dcf0dd564c80ee97092afae63b3e1 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 16 May 2024 09:16:08 +0100 Subject: [PATCH 08/43] Rename files --- .../{step1_2d_NS_elementTypes.i => step1_2d_INSAD_elementTypes.i} | 0 .../{step1_3d_NS_elementTypes.i => step1_3d_INSAD_elementTypes.i} | 0 .../{step2_2d_NS_LorentzForce.i => step2_2d_INSAD_LorentzForce.i} | 0 3 files changed, 0 insertions(+), 0 deletions(-) rename examples/mhd/charge_conserving/development/{step1_2d_NS_elementTypes.i => step1_2d_INSAD_elementTypes.i} (100%) rename examples/mhd/charge_conserving/development/{step1_3d_NS_elementTypes.i => step1_3d_INSAD_elementTypes.i} (100%) rename examples/mhd/charge_conserving/development/{step2_2d_NS_LorentzForce.i => step2_2d_INSAD_LorentzForce.i} (100%) diff --git a/examples/mhd/charge_conserving/development/step1_2d_NS_elementTypes.i b/examples/mhd/charge_conserving/development/step1_2d_INSAD_elementTypes.i similarity index 100% rename from examples/mhd/charge_conserving/development/step1_2d_NS_elementTypes.i rename to examples/mhd/charge_conserving/development/step1_2d_INSAD_elementTypes.i diff --git a/examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i b/examples/mhd/charge_conserving/development/step1_3d_INSAD_elementTypes.i similarity index 100% rename from examples/mhd/charge_conserving/development/step1_3d_NS_elementTypes.i rename to examples/mhd/charge_conserving/development/step1_3d_INSAD_elementTypes.i diff --git a/examples/mhd/charge_conserving/development/step2_2d_NS_LorentzForce.i b/examples/mhd/charge_conserving/development/step2_2d_INSAD_LorentzForce.i similarity index 100% rename from examples/mhd/charge_conserving/development/step2_2d_NS_LorentzForce.i rename to examples/mhd/charge_conserving/development/step2_2d_INSAD_LorentzForce.i From b8d57955be284de1c931d63041de8300798c6bf3 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 16 May 2024 11:07:36 +0100 Subject: [PATCH 09/43] Add INS (non-AD) version --- .../development/step1_3d_INS_elementTypes.i | 320 ++++++++++++++++++ 1 file changed, 320 insertions(+) create mode 100644 examples/mhd/charge_conserving/development/step1_3d_INS_elementTypes.i diff --git a/examples/mhd/charge_conserving/development/step1_3d_INS_elementTypes.i b/examples/mhd/charge_conserving/development/step1_3d_INS_elementTypes.i new file mode 100644 index 0000000..1c56622 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step1_3d_INS_elementTypes.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' +[] From 226506fed6036aff97e6d20af513fe2b1386edd5 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 16 May 2024 11:15:19 +0100 Subject: [PATCH 10/43] Tidy up case file directory structure --- .../2d_INSAD.i} | 0 .../3d_INS.i} | 0 .../3d_INSAD.i} | 1 - .../2d_INSAD.i} | 0 4 files changed, 1 deletion(-) rename examples/mhd/charge_conserving/development/{step1_2d_INSAD_elementTypes.i => step1_NavierStokes_elementTypes/2d_INSAD.i} (100%) rename examples/mhd/charge_conserving/development/{step1_3d_INS_elementTypes.i => step1_NavierStokes_elementTypes/3d_INS.i} (100%) rename examples/mhd/charge_conserving/development/{step1_3d_INSAD_elementTypes.i => step1_NavierStokes_elementTypes/3d_INSAD.i} (99%) rename examples/mhd/charge_conserving/development/{step2_2d_INSAD_LorentzForce.i => step2_NavierStokes_LorentzForce/2d_INSAD.i} (100%) diff --git a/examples/mhd/charge_conserving/development/step1_2d_INSAD_elementTypes.i b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/2d_INSAD.i similarity index 100% rename from examples/mhd/charge_conserving/development/step1_2d_INSAD_elementTypes.i rename to examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/2d_INSAD.i diff --git a/examples/mhd/charge_conserving/development/step1_3d_INS_elementTypes.i b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INS.i similarity index 100% rename from examples/mhd/charge_conserving/development/step1_3d_INS_elementTypes.i rename to examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INS.i diff --git a/examples/mhd/charge_conserving/development/step1_3d_INSAD_elementTypes.i b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INSAD.i similarity index 99% rename from examples/mhd/charge_conserving/development/step1_3d_INSAD_elementTypes.i rename to examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INSAD.i index 2c4f831..ffc5d0d 100644 --- a/examples/mhd/charge_conserving/development/step1_3d_INSAD_elementTypes.i +++ b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/3d_INSAD.i @@ -263,5 +263,4 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [Outputs] exodus = true - execute_on = 'nonlinear' [] diff --git a/examples/mhd/charge_conserving/development/step2_2d_INSAD_LorentzForce.i b/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/2d_INSAD.i similarity index 100% rename from examples/mhd/charge_conserving/development/step2_2d_INSAD_LorentzForce.i rename to examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/2d_INSAD.i From 0982349d8ac14f00f5e83c635428a6fb4eec69eb Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 16 May 2024 15:11:00 +0100 Subject: [PATCH 11/43] Delete non-functional 2d Lorentz force case, add 3d case --- .../2d_INSAD.i | 217 ------------ .../3d_INSAD.i | 311 ++++++++++++++++++ 2 files changed, 311 insertions(+), 217 deletions(-) delete mode 100644 examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/2d_INSAD.i create mode 100644 examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i diff --git a/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/2d_INSAD.i b/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/2d_INSAD.i deleted file mode 100644 index 331c329..0000000 --- a/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/2d_INSAD.i +++ /dev/null @@ -1,217 +0,0 @@ -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 - [] -[] - -[AuxVariables] - [magneticField] - family = LAGRANGE_VEC - order = FIRST - [] - [currentDensity] - family = LAGRANGE_VEC - 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' - [] - [magneticFieldFunction] - type = ParsedVectorFunction - expression_x = '0' - expression_y = '20' - [] - [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' - 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 - [] - - [momentum_lorentz] - type = IMHDADMomentumLorentz - variable = velocity - currentDensity = currentDensity - magneticField = magneticField - [] -[] - -[AuxKernels] - [magneticFieldKernel] - type = VectorFunctionAux - variable = magneticField - function = magneticFieldFunction - execute_on = INITIAL - [] - [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 - 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/step2_NavierStokes_LorentzForce/3d_INSAD.i b/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i new file mode 100644 index 0000000..14028b3 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i @@ -0,0 +1,311 @@ +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] + [magneticField] + family = LAGRANGE_VEC + order = FIRST + [] + [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 + magneticField = magneticField + [] +[] + +[AuxKernels] + [magneticFieldKernel] + type = VectorFunctionAux + variable = magneticField + function = magneticFieldFunction + execute_on = INITIAL + [] + [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 +[] From efbead440a4da96ee9bf54d0817766c7c76774fa Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 16 May 2024 15:12:03 +0100 Subject: [PATCH 12/43] Fix incorrect sign of Lorentz force term (opposite orientation) --- src/kernels/IMHDADMomentumLorentz.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kernels/IMHDADMomentumLorentz.C b/src/kernels/IMHDADMomentumLorentz.C index b6c59ff..b9d2244 100644 --- a/src/kernels/IMHDADMomentumLorentz.C +++ b/src/kernels/IMHDADMomentumLorentz.C @@ -26,5 +26,5 @@ IMHDADMomentumLorentz::IMHDADMomentumLorentz(const InputParameters & parameters) ADRealVectorValue IMHDADMomentumLorentz::precomputeQpResidual() { - return _current_density[_qp].cross(_magnetic_field[_qp]); + return -_current_density[_qp].cross(_magnetic_field[_qp]); } From 1a921679ea0ea63b758b2ff83c47356a3319f339 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 16 May 2024 15:12:39 +0100 Subject: [PATCH 13/43] Testing traction form of viscous momentum term --- .../development/step2_NavierStokes_LorentzForce/3d_INSAD.i | 1 + 1 file changed, 1 insertion(+) 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 index 14028b3..a95e01d 100644 --- a/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i +++ b/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i @@ -253,6 +253,7 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [momentum_viscous] type = INSADMomentumViscous variable = velocity + viscous_form = traction [] [momentum_pressure] From 084f7f39599ade63124b91d3e214cba684027592 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 16 May 2024 15:13:48 +0100 Subject: [PATCH 14/43] Switch back to default laplace form of viscous equation due to concerns about y and z velocity near outlet in traction form --- .../development/step2_NavierStokes_LorentzForce/3d_INSAD.i | 1 - 1 file changed, 1 deletion(-) 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 index a95e01d..14028b3 100644 --- a/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i +++ b/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i @@ -253,7 +253,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [momentum_viscous] type = INSADMomentumViscous variable = velocity - viscous_form = traction [] [momentum_pressure] From 82efbfd3abc095f1eae9f72652fc3a047e537985 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 16 May 2024 17:15:35 +0100 Subject: [PATCH 15/43] Add initial electrostatic implementation --- .../step3_electrostatic/3d_electrostatic.i | 121 ++++++++++++++++++ include/kernels/IMHDADCurrentDensity.h | 20 +++ src/kernels/IMHDADCurrentDensity.C | 28 ++++ 3 files changed, 169 insertions(+) create mode 100644 examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic.i create mode 100644 include/kernels/IMHDADCurrentDensity.h create mode 100644 src/kernels/IMHDADCurrentDensity.C diff --git a/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic.i b/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic.i new file mode 100644 index 0000000..d20a229 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic.i @@ -0,0 +1,121 @@ +[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 + [] +[] + +[Variables] + [currentDensity] + family = RAVIART_THOMAS + order = FIRST + [] + [electricPotential] + family = MONOMIAL + order = CONSTANT + # order = FIRST + [] +[] + +[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] + # [epot_live] + # type = PenaltyDirichletBC + # variable = electricPotential + # value = 100 + # boundary = 'left' + # penalty = 1e5 + # [] + # [epot_ground] + # type = PenaltyDirichletBC + # variable = electricPotential + # value = 0 + # boundary = 'right' + # penalty = 1e5 + # [] + [current_wall_in] + type = VectorDivPenaltyDirichletBC + variable = currentDensity + function_x = cos(y*3.14159/2)*cos(z*3.14159/2) + penalty = 1e8 + boundary = 'left' + [] + [current_wall_out] + type = VectorDivPenaltyDirichletBC + variable = currentDensity + function_x = 0 + penalty = 1e8 + boundary = 'right' + [] +[] + +# [Postprocessors] +# [L2Error] +# type = ElementVectorL2Error +# variable = u +# function = f +# [] +# [HDivSemiError] +# type = ElementHDivSemiError +# variable = u +# function = f +# [] +# [HDivError] +# type = ElementHDivError +# variable = u +# function = f +# [] +# [] + +[Materials] + [const] + type = ADGenericConstantMaterial + prop_names = 'conductivity' + prop_values = '1' + [] +[] + +[Preconditioning] + [SMP] + type = SMP + full = true + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + petsc_options_iname = '-pc_type -ksp_rtol -ksp_norm_type' + petsc_options_value = ' jacobi 1e-12 preconditioned' +[] + +[Outputs] + exodus = true +[] \ No newline at end of file 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/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]; +} From 013121a78832d4a90b0eb9067302597ef06bce05 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Mon, 20 May 2024 18:02:52 +0100 Subject: [PATCH 16/43] Add new 2d electrostatic test cases and update 3d electrostatic case --- .../2d_electrostatic_gradient.i | 121 ++++++++++++++++++ .../2d_electrostatic_solenoidal.i | 115 +++++++++++++++++ ...static.i => 3d_electrostatic_solenoidal.i} | 90 +++++++------ 3 files changed, 280 insertions(+), 46 deletions(-) create mode 100644 examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_gradient.i create mode 100644 examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_solenoidal.i rename examples/mhd/charge_conserving/development/step3_electrostatic/{3d_electrostatic.i => 3d_electrostatic_solenoidal.i} (50%) 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..6067232 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_gradient.i @@ -0,0 +1,121 @@ +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..4fbf562 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_solenoidal.i @@ -0,0 +1,115 @@ +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' + [] +[] + +[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.i b/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic_solenoidal.i similarity index 50% rename from examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic.i rename to examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic_solenoidal.i index d20a229..4dfc767 100644 --- a/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic.i +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic_solenoidal.i @@ -1,3 +1,5 @@ +sigma = 1 + [Mesh] [gmg] type = GeneratedMeshGenerator @@ -15,6 +17,16 @@ [] [] +[Functions] + [solenoid] + type = ParsedVectorFunction + expression_x = y + expression_y = -x + expression_z = 0 + div = 0 + [] +[] + [Variables] [currentDensity] family = RAVIART_THOMAS @@ -23,7 +35,6 @@ [electricPotential] family = MONOMIAL order = CONSTANT - # order = FIRST [] [] @@ -38,6 +49,11 @@ coupled_scalar_variable = electricPotential coeff = -1 [] + [currentDensitySource] + type = VectorBodyForce + variable = currentDensity + [] + [divergenceFree] type = DivField variable = electricPotential @@ -46,59 +62,40 @@ [] [BCs] - # [epot_live] - # type = PenaltyDirichletBC - # variable = electricPotential - # value = 100 - # boundary = 'left' - # penalty = 1e5 - # [] - # [epot_ground] - # type = PenaltyDirichletBC - # variable = electricPotential - # value = 0 - # boundary = 'right' - # penalty = 1e5 - # [] - [current_wall_in] - type = VectorDivPenaltyDirichletBC + [current_wall_value] + type = VectorPenaltyDirichletBC variable = currentDensity - function_x = cos(y*3.14159/2)*cos(z*3.14159/2) - penalty = 1e8 - boundary = 'left' + 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' [] - [current_wall_out] - type = VectorDivPenaltyDirichletBC +[] + +[Postprocessors] + [L2Error] + type = ElementVectorL2Error + variable = currentDensity + function = solenoid + [] + [HDivSemiError] + type = ElementHDivSemiError variable = currentDensity - function_x = 0 - penalty = 1e8 - boundary = 'right' + function = solenoid + [] + [HDivError] + type = ElementHDivError + variable = currentDensity + function = solenoid [] [] -# [Postprocessors] -# [L2Error] -# type = ElementVectorL2Error -# variable = u -# function = f -# [] -# [HDivSemiError] -# type = ElementHDivSemiError -# variable = u -# function = f -# [] -# [HDivError] -# type = ElementHDivError -# variable = u -# function = f -# [] -# [] - [Materials] [const] type = ADGenericConstantMaterial prop_names = 'conductivity' - prop_values = '1' + prop_values = '${sigma}' [] [] @@ -112,10 +109,11 @@ [Executioner] type = Steady solve_type = NEWTON - petsc_options_iname = '-pc_type -ksp_rtol -ksp_norm_type' - petsc_options_value = ' jacobi 1e-12 preconditioned' + petsc_options_iname = '-pc_type' + petsc_options_value = ' bjacobi' [] [Outputs] exodus = true + execute_on = 'NONLINEAR' [] \ No newline at end of file From 528e1a3efdb192c189bbea0c86a44f57594153ff Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Tue, 21 May 2024 11:42:27 +0100 Subject: [PATCH 17/43] Add missing function argument to VectorBodyForce kernel --- .../step3_electrostatic/3d_electrostatic_solenoidal.i | 1 + 1 file changed, 1 insertion(+) 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 index 4dfc767..575fc5d 100644 --- a/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic_solenoidal.i +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic_solenoidal.i @@ -52,6 +52,7 @@ sigma = 1 [currentDensitySource] type = VectorBodyForce variable = currentDensity + function = solenoid [] [divergenceFree] From 174831168fa8860339329df45e5fc65a21c8ef69 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Tue, 21 May 2024 17:38:11 +0100 Subject: [PATCH 18/43] Add UxB electric potential coupling term and test case --- .../3d_coupledAuxVelocity.i | 277 ++++++++++++++++++ include/kernels/IMHDADCurrentUxB.h | 21 ++ src/kernels/IMHDADCurrentUxB.C | 30 ++ 3 files changed, 328 insertions(+) create mode 100644 examples/mhd/charge_conserving/development/step4_coupledAuxVelocity/3d_coupledAuxVelocity.i create mode 100644 include/kernels/IMHDADCurrentUxB.h create mode 100644 src/kernels/IMHDADCurrentUxB.C 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..9e95bf9 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step4_coupledAuxVelocity/3d_coupledAuxVelocity.i @@ -0,0 +1,277 @@ +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 + [] + [magneticField] + 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 + magneticField = magneticField + [] + + [divergenceFree] + type = DivField + variable = electricPotential + coupled_vector_variable = currentDensity + [] +[] + +[AuxKernels] + [magneticFieldKernel] + type = VectorFunctionAux + variable = magneticField + function = magneticFieldFunction + execute_on = INITIAL + [] + [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/include/kernels/IMHDADCurrentUxB.h b/include/kernels/IMHDADCurrentUxB.h new file mode 100644 index 0000000..6bea65b --- /dev/null +++ b/include/kernels/IMHDADCurrentUxB.h @@ -0,0 +1,21 @@ +#pragma once + +#include "ADKernelValue.h" + +/** + * This class computes the residual and Jacobian contributions for the + * U cross B term of the electric current density equation. + */ +class IMHDADCurrentUxB : public ADVectorKernelValue +{ +public: + static InputParameters validParams(); + + IMHDADCurrentUxB(const InputParameters & parameters); + +protected: + virtual ADRealVectorValue precomputeQpResidual() override; + + const ADVectorVariableValue & _velocity; + const ADVectorVariableValue & _magnetic_field; +}; diff --git a/src/kernels/IMHDADCurrentUxB.C b/src/kernels/IMHDADCurrentUxB.C new file mode 100644 index 0000000..382f481 --- /dev/null +++ b/src/kernels/IMHDADCurrentUxB.C @@ -0,0 +1,30 @@ +#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.addRequiredCoupledVar("magneticField", "The variable representing the magnetic field"); + + return params; +} + +IMHDADCurrentUxB::IMHDADCurrentUxB(const InputParameters & parameters) + : ADVectorKernelValue(parameters), + _velocity(adCoupledVectorValue("velocity")), + _magnetic_field(adCoupledVectorValue("magneticField")) +{ +} + +ADRealVectorValue +IMHDADCurrentUxB::precomputeQpResidual() +{ + return _velocity[_qp].cross(_magnetic_field[_qp]); +} From 6c40b3770dadf6a365f446b58364d4f03d9286b3 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 22 May 2024 11:11:17 +0100 Subject: [PATCH 19/43] Add WIP fully coupled example --- .../3d_coupledAuxVelocity.i | 336 ++++++++++++++++++ 1 file changed, 336 insertions(+) create mode 100644 examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupledAuxVelocity.i diff --git a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupledAuxVelocity.i b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupledAuxVelocity.i new file mode 100644 index 0000000..b6fa7c5 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupledAuxVelocity.i @@ -0,0 +1,336 @@ +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] + [velocity] + family = LAGRANGE_VEC + order = SECOND + [] + [pressure] + family = MONOMIAL + order = FIRST + [] + [currentDensity] + family = RAVIART_THOMAS + order = FIRST + [] + [electricPotential] + family = MONOMIAL + order = CONSTANT + [] +[] + +[AuxVariables] + [magneticField] + family = LAGRANGE_VEC + order = FIRST + [] +[] + +[ICs] + [velocityIC] + type = VectorConstantIC + x_value = ${U_AVG} + y_value = 1e-15 + z_value = 1e-15 + 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 + magneticField = magneticField + [] + + [currentDensityConductivity] + type = IMHDADCurrentDensity + variable = currentDensity + [] + [electricField] + type = GradField + variable = currentDensity + coupled_scalar_variable = electricPotential + coeff = -1 + [] + [currentDensityCoupling] + type = IMHDADCurrentUxB + variable = currentDensity + velocity = velocity + magneticField = magneticField + [] + + [divergenceFree] + type = DivField + variable = electricPotential + coupled_vector_variable = currentDensity + [] +[] + +[AuxKernels] + [magneticFieldKernel] + type = VectorFunctionAux + variable = magneticField + function = magneticFieldFunction + execute_on = INITIAL + [] +[] + +[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] + # 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 = 'rho mu conductivity' + prop_values = '1 1 ${sigma}' + [] + [ins_mat] + type = INSADTauMaterial + velocity = velocity + pressure = pressure + [] +[] + +[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' +[] \ No newline at end of file From 7bdf6c86ea1ee34c68ccbc4a2ea388003b8264f7 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 23 May 2024 14:07:21 +0100 Subject: [PATCH 20/43] Rename fully coupled case (previous name incorrectly included AuxVelocity) --- .../{3d_coupledAuxVelocity.i => 3d_coupled.i} | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) rename examples/mhd/charge_conserving/development/step5_MHD_coupling/{3d_coupledAuxVelocity.i => 3d_coupled.i} (97%) diff --git a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupledAuxVelocity.i b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i similarity index 97% rename from examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupledAuxVelocity.i rename to examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i index b6fa7c5..95d0d99 100644 --- a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupledAuxVelocity.i +++ b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i @@ -196,10 +196,12 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [ICs] [velocityIC] - type = VectorConstantIC - x_value = ${U_AVG} - y_value = 1e-15 - z_value = 1e-15 + # type = VectorConstantIC + # x_value = ${U_AVG} + # y_value = 1e-15 + # z_value = 1e-15 + type = VectorFunctionIC + function = velocityFunction variable = velocity [] [] @@ -327,10 +329,10 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} l_max_its = 1000 nl_max_its = 1000 petsc_options_iname = '-pc_type' - petsc_options_value = ' bjacobi' + petsc_options_value = ' ksp' [] [Outputs] exodus = true execute_on = 'NONLINEAR' -[] \ No newline at end of file +[] From 0a14084bd6c3e5cb9505d6d8c724159c6418cce6 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Fri, 24 May 2024 16:51:48 +0100 Subject: [PATCH 21/43] Change sign of IMHDADCurrentUxB residual --- src/kernels/IMHDADCurrentUxB.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kernels/IMHDADCurrentUxB.C b/src/kernels/IMHDADCurrentUxB.C index 382f481..d0b609d 100644 --- a/src/kernels/IMHDADCurrentUxB.C +++ b/src/kernels/IMHDADCurrentUxB.C @@ -26,5 +26,5 @@ IMHDADCurrentUxB::IMHDADCurrentUxB(const InputParameters & parameters) ADRealVectorValue IMHDADCurrentUxB::precomputeQpResidual() { - return _velocity[_qp].cross(_magnetic_field[_qp]); + return -_velocity[_qp].cross(_magnetic_field[_qp]); } From 0637ea640235e7b0c218a6f178c4ddfecbaca260 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Fri, 24 May 2024 16:52:14 +0100 Subject: [PATCH 22/43] Add WIP INSADAugmentedLagrangian kernel --- include/kernels/INSADAugmentedLagrangian.h | 28 +++++++++++++++++++ src/kernels/INSADAugmentedLagrangian.C | 31 ++++++++++++++++++++++ 2 files changed, 59 insertions(+) create mode 100644 include/kernels/INSADAugmentedLagrangian.h create mode 100644 src/kernels/INSADAugmentedLagrangian.C diff --git a/include/kernels/INSADAugmentedLagrangian.h b/include/kernels/INSADAugmentedLagrangian.h new file mode 100644 index 0000000..bc1bcff --- /dev/null +++ b/include/kernels/INSADAugmentedLagrangian.h @@ -0,0 +1,28 @@ +#pragma once + +#include "ADKernel.h" + +/** + * This class computes the residual and Jacobian contributions for the + * augmented Lagrangian stabilisation term. + */ +class INSADAugmentedLagrangian : public ADVectorKernel +{ +public: + static InputParameters validParams(); + + INSADAugmentedLagrangian(const InputParameters & parameters); + +protected: + // virtual ADRealVectorValue computeQpResidual() override; + virtual ADReal computeQpResidual() override; + + /// div of the vector variable + const VectorVariableDivergence & _div_u; + + /// div of the vector test function + const VectorVariableTestDivergence & _div_test; + + /// scalar coefficient + ADReal _coeff; +}; diff --git a/src/kernels/INSADAugmentedLagrangian.C b/src/kernels/INSADAugmentedLagrangian.C new file mode 100644 index 0000000..46c2f5b --- /dev/null +++ b/src/kernels/INSADAugmentedLagrangian.C @@ -0,0 +1,31 @@ +#include "INSADAugmentedLagrangian.h" + +registerMooseObject("ProteusApp", INSADAugmentedLagrangian); + +InputParameters +INSADAugmentedLagrangian::validParams() +{ + InputParameters params = ADVectorKernel::validParams(); + params.addClassDescription( + "The augmented Lagrangian term, " + "with the weak form of ($\\alpha(\\nabla \\cdot \\phi_u, \\nabla \\cdot u)$), " + "where $\\alpha>0$ is the scalar AL-stabilisation parameter. " + "The Jacobian is computed using automatic differentiation."); + params.addParam("coeff", 1.0, "The AL-stabilisation parameter") + + return params; +} + +INSADAugmentedLagrangian::INSADAugmentedLagrangian(const InputParameters & parameters) + : ADVectorKernel(parameters), + + _div_test(_var.divPhi()), + _coeff(getParam("coeff")) +{ +} + +ADReal +INSADAugmentedLagrangian::precomputeQpResidual() +{ + return _coeff * _div_u[_qp] * _div_test[_qp]; +} From b4cc3f94c2c2015a15414cfca9139e13e30362b6 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Fri, 24 May 2024 16:53:17 +0100 Subject: [PATCH 23/43] Delete INSADAugmentedLagrangian kernels; found INSADMomentumGradDiv already implements the same functionality --- include/kernels/INSADAugmentedLagrangian.h | 28 ------------------- src/kernels/INSADAugmentedLagrangian.C | 31 ---------------------- 2 files changed, 59 deletions(-) delete mode 100644 include/kernels/INSADAugmentedLagrangian.h delete mode 100644 src/kernels/INSADAugmentedLagrangian.C diff --git a/include/kernels/INSADAugmentedLagrangian.h b/include/kernels/INSADAugmentedLagrangian.h deleted file mode 100644 index bc1bcff..0000000 --- a/include/kernels/INSADAugmentedLagrangian.h +++ /dev/null @@ -1,28 +0,0 @@ -#pragma once - -#include "ADKernel.h" - -/** - * This class computes the residual and Jacobian contributions for the - * augmented Lagrangian stabilisation term. - */ -class INSADAugmentedLagrangian : public ADVectorKernel -{ -public: - static InputParameters validParams(); - - INSADAugmentedLagrangian(const InputParameters & parameters); - -protected: - // virtual ADRealVectorValue computeQpResidual() override; - virtual ADReal computeQpResidual() override; - - /// div of the vector variable - const VectorVariableDivergence & _div_u; - - /// div of the vector test function - const VectorVariableTestDivergence & _div_test; - - /// scalar coefficient - ADReal _coeff; -}; diff --git a/src/kernels/INSADAugmentedLagrangian.C b/src/kernels/INSADAugmentedLagrangian.C deleted file mode 100644 index 46c2f5b..0000000 --- a/src/kernels/INSADAugmentedLagrangian.C +++ /dev/null @@ -1,31 +0,0 @@ -#include "INSADAugmentedLagrangian.h" - -registerMooseObject("ProteusApp", INSADAugmentedLagrangian); - -InputParameters -INSADAugmentedLagrangian::validParams() -{ - InputParameters params = ADVectorKernel::validParams(); - params.addClassDescription( - "The augmented Lagrangian term, " - "with the weak form of ($\\alpha(\\nabla \\cdot \\phi_u, \\nabla \\cdot u)$), " - "where $\\alpha>0$ is the scalar AL-stabilisation parameter. " - "The Jacobian is computed using automatic differentiation."); - params.addParam("coeff", 1.0, "The AL-stabilisation parameter") - - return params; -} - -INSADAugmentedLagrangian::INSADAugmentedLagrangian(const InputParameters & parameters) - : ADVectorKernel(parameters), - - _div_test(_var.divPhi()), - _coeff(getParam("coeff")) -{ -} - -ADReal -INSADAugmentedLagrangian::precomputeQpResidual() -{ - return _coeff * _div_u[_qp] * _div_test[_qp]; -} From 6a94965c33b71b89076ffd6383b33b39a202a138 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Fri, 24 May 2024 16:55:04 +0100 Subject: [PATCH 24/43] Increase y resolution for higher Ha cases, add optional FSP preconditioning framework (doesn't seem to help so far) --- .../step5_MHD_coupling/3d_coupled.i | 36 +++++++++++++++++-- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i index 95d0d99..c63cbb7 100644 --- a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i +++ b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i @@ -3,7 +3,7 @@ sigma = 1 U_AVG = 1 N_X = 50 -N_Y_half = 10 +N_Y_half = 20 N_Z_half = 10 ELEMENT_TYPE = HEX27 @@ -316,9 +316,41 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [] [Preconditioning] + active = 'SMP' + # active = 'FSP' [SMP] type = SMP full = true + petsc_options_iname = '-pc_type' + petsc_options_value = ' cholesky' + [] + [FSP] + type = FSP + topsplit = 'splitting' + [splitting] + splitting = 'velocity pressure currentDensity electricPotential' + splitting_type = additive + [] + [velocity] + vars = 'velocity' + petsc_options_iname = '-pc_type' + petsc_options_value = ' cholesky' + [] + [pressure] + vars = 'pressure' + petsc_options_iname = '-pc_type' + petsc_options_value = ' cholesky' + [] + [currentDensity] + vars = 'currentDensity' + petsc_options_iname = '-pc_type' + petsc_options_value = ' cholesky' + [] + [electricPotential] + vars = 'electricPotential' + petsc_options_iname = '-pc_type' + petsc_options_value = ' cholesky' + [] [] [] @@ -328,8 +360,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} automatic_scaling = true l_max_its = 1000 nl_max_its = 1000 - petsc_options_iname = '-pc_type' - petsc_options_value = ' ksp' [] [Outputs] From c2e8b1fb6593df2c89bf01bfc42329215a6bb374 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Fri, 24 May 2024 17:19:58 +0100 Subject: [PATCH 25/43] Add augmented Lagrangian terms for velocity and current density to demonstrate implementation; neither seems to help --- .../development/step5_MHD_coupling/3d_coupled.i | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i index c63cbb7..156817c 100644 --- a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i +++ b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i @@ -232,6 +232,11 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} currentDensity = currentDensity magneticField = magneticField [] + [momentum_augmented_lagrangian] + type = INSADMomentumGradDiv + variable = velocity + gamma = 1 + [] [currentDensityConductivity] type = IMHDADCurrentDensity @@ -249,6 +254,11 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} velocity = velocity magneticField = magneticField [] + [current_augmented_lagrangian] + type = INSADMomentumGradDiv + variable = currentDensity + gamma = 1 + [] [divergenceFree] type = DivField From 3845e22cf55107a36b3ceaa00d209b5609d5df45 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Fri, 24 May 2024 17:23:09 +0100 Subject: [PATCH 26/43] Remove AL terms --- .../development/step5_MHD_coupling/3d_coupled.i | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i index 156817c..c63cbb7 100644 --- a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i +++ b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i @@ -232,11 +232,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} currentDensity = currentDensity magneticField = magneticField [] - [momentum_augmented_lagrangian] - type = INSADMomentumGradDiv - variable = velocity - gamma = 1 - [] [currentDensityConductivity] type = IMHDADCurrentDensity @@ -254,11 +249,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} velocity = velocity magneticField = magneticField [] - [current_augmented_lagrangian] - type = INSADMomentumGradDiv - variable = currentDensity - gamma = 1 - [] [divergenceFree] type = DivField From e46eed78739f6ae24ba319ee738e5530d5cb04d3 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 30 May 2024 16:28:45 +0100 Subject: [PATCH 27/43] Switch Lorentz force kernel to use a vector function for B rather than an AuxVariable --- .../3d_INSAD.i | 12 +------- include/kernels/IMHDADMomentumLorentz.h | 14 ++++++++-- src/kernels/IMHDADMomentumLorentz.C | 28 +++++++++++++++++-- 3 files changed, 38 insertions(+), 16 deletions(-) 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 index 14028b3..5587825 100644 --- a/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i +++ b/examples/mhd/charge_conserving/development/step2_NavierStokes_LorentzForce/3d_INSAD.i @@ -161,10 +161,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [] [AuxVariables] - [magneticField] - family = LAGRANGE_VEC - order = FIRST - [] [currentDensity] family = LAGRANGE_VEC order = FIRST @@ -266,17 +262,11 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} type = IMHDADMomentumLorentz variable = velocity currentDensity = currentDensity - magneticField = magneticField + magneticFieldFunction = magneticFieldFunction [] [] [AuxKernels] - [magneticFieldKernel] - type = VectorFunctionAux - variable = magneticField - function = magneticFieldFunction - execute_on = INITIAL - [] [currentDensityKernel] type = VectorFunctionAux variable = currentDensity diff --git a/include/kernels/IMHDADMomentumLorentz.h b/include/kernels/IMHDADMomentumLorentz.h index 0b2eff5..a04a7fe 100644 --- a/include/kernels/IMHDADMomentumLorentz.h +++ b/include/kernels/IMHDADMomentumLorentz.h @@ -1,10 +1,14 @@ #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. + * 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 { @@ -18,5 +22,11 @@ class IMHDADMomentumLorentz : public ADVectorKernelValue const ADVectorVariableValue & _current_density; - const ADVectorVariableValue & _magnetic_field; + /// 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/src/kernels/IMHDADMomentumLorentz.C b/src/kernels/IMHDADMomentumLorentz.C index b9d2244..3bcb05f 100644 --- a/src/kernels/IMHDADMomentumLorentz.C +++ b/src/kernels/IMHDADMomentumLorentz.C @@ -11,7 +11,15 @@ IMHDADMomentumLorentz::validParams() "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.addRequiredCoupledVar("magneticField", "The variable representing the magnetic field"); + params.addParam("magneticFieldFunction", + "A function for the imposed magnetic field that defines a vectorValue method. " + "This cannot be supplied with the component parameters."); + 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; } @@ -19,12 +27,26 @@ IMHDADMomentumLorentz::validParams() IMHDADMomentumLorentz::IMHDADMomentumLorentz(const InputParameters & parameters) : ADVectorKernelValue(parameters), _current_density(adCoupledVectorValue("currentDensity")), - _magnetic_field(adCoupledVectorValue("magneticField")) + _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() { - return -_current_density[_qp].cross(_magnetic_field[_qp]); + 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]))); } From b1cdd52444c0af7c04a4c1a1b5b1bf45b2af572a Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 30 May 2024 16:39:19 +0100 Subject: [PATCH 28/43] Update to use new version of Lorentz force kernel --- .../step5_MHD_coupling/3d_coupled.i | 18 +----------------- 1 file changed, 1 insertion(+), 17 deletions(-) diff --git a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i index c63cbb7..6a0e0a8 100644 --- a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i +++ b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i @@ -187,13 +187,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [] [] -[AuxVariables] - [magneticField] - family = LAGRANGE_VEC - order = FIRST - [] -[] - [ICs] [velocityIC] # type = VectorConstantIC @@ -230,7 +223,7 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} type = IMHDADMomentumLorentz variable = velocity currentDensity = currentDensity - magneticField = magneticField + magneticFieldFunction = magneticFieldFunction [] [currentDensityConductivity] @@ -257,15 +250,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [] [] -[AuxKernels] - [magneticFieldKernel] - type = VectorFunctionAux - variable = magneticField - function = magneticFieldFunction - execute_on = INITIAL - [] -[] - [BCs] [inlet] type = VectorFunctionDirichletBC From 24c56b7f97a57b4988a76209a8686f672a95aee0 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Tue, 4 Jun 2024 13:16:08 +0100 Subject: [PATCH 29/43] Switch UxB kernel to use vector function for B, update coupled example --- .../step5_MHD_coupling/3d_coupled.i | 4 +-- include/kernels/IMHDADCurrentUxB.h | 15 ++++++++-- src/kernels/IMHDADCurrentUxB.C | 28 +++++++++++++++++-- 3 files changed, 40 insertions(+), 7 deletions(-) diff --git a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i index 6a0e0a8..271f6cc 100644 --- a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i +++ b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i @@ -240,7 +240,7 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} type = IMHDADCurrentUxB variable = currentDensity velocity = velocity - magneticField = magneticField + magneticFieldFunction = magneticFieldFunction [] [divergenceFree] @@ -306,7 +306,7 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} type = SMP full = true petsc_options_iname = '-pc_type' - petsc_options_value = ' cholesky' + petsc_options_value = ' lu' [] [FSP] type = FSP diff --git a/include/kernels/IMHDADCurrentUxB.h b/include/kernels/IMHDADCurrentUxB.h index 6bea65b..e6bd4b9 100644 --- a/include/kernels/IMHDADCurrentUxB.h +++ b/include/kernels/IMHDADCurrentUxB.h @@ -1,10 +1,14 @@ #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. + * 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 { @@ -17,5 +21,12 @@ class IMHDADCurrentUxB : public ADVectorKernelValue virtual ADRealVectorValue precomputeQpResidual() override; const ADVectorVariableValue & _velocity; - const ADVectorVariableValue & _magnetic_field; + + /// 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/src/kernels/IMHDADCurrentUxB.C b/src/kernels/IMHDADCurrentUxB.C index d0b609d..85bc03a 100644 --- a/src/kernels/IMHDADCurrentUxB.C +++ b/src/kernels/IMHDADCurrentUxB.C @@ -11,7 +11,15 @@ IMHDADCurrentUxB::validParams() "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.addRequiredCoupledVar("magneticField", "The variable representing the magnetic field"); + params.addParam("magneticFieldFunction", + "A function for the imposed magnetic field that defines a vectorValue method. " + "This cannot be supplied with the component parameters."); + 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; } @@ -19,12 +27,26 @@ IMHDADCurrentUxB::validParams() IMHDADCurrentUxB::IMHDADCurrentUxB(const InputParameters & parameters) : ADVectorKernelValue(parameters), _velocity(adCoupledVectorValue("velocity")), - _magnetic_field(adCoupledVectorValue("magneticField")) + _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() { - return -_velocity[_qp].cross(_magnetic_field[_qp]); + 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]))); } From d0670e15c55def1e3ce20968fe63c76058289fd1 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 6 Jun 2024 11:29:06 +0100 Subject: [PATCH 30/43] Specify shercliff case in coupled example filename --- .../step5_MHD_coupling/{3d_coupled.i => 3d_coupled_shercliff.i} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/mhd/charge_conserving/development/step5_MHD_coupling/{3d_coupled.i => 3d_coupled_shercliff.i} (100%) diff --git a/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled_shercliff.i similarity index 100% rename from examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled.i rename to examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled_shercliff.i From 2346dfe45f4f4acb40d403933b5c47b89f36c278 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 6 Jun 2024 11:29:56 +0100 Subject: [PATCH 31/43] Switch to LAGRANGE FIRST pressure --- .../development/step5_MHD_coupling/3d_coupled_shercliff.i | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 index 271f6cc..78b2ca8 100644 --- 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 @@ -174,7 +174,7 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} order = SECOND [] [pressure] - family = MONOMIAL + family = LAGRANGE order = FIRST [] [currentDensity] From 563253972429cd6f91a8f629bc77d727fbcc636d Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 6 Jun 2024 11:47:19 +0100 Subject: [PATCH 32/43] Disable automatic scaling, remove resolved comments about uncertainty in J.n BC --- .../development/step5_MHD_coupling/3d_coupled_shercliff.i | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) 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 index 78b2ca8..b8c59c1 100644 --- 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 @@ -271,11 +271,7 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} penalty = 1e5 [] [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 + # Setting J.n=0 on the walls, inlet and outlet. type = VectorDivPenaltyDirichletBC variable = currentDensity function_x = 0 @@ -341,7 +337,7 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [Executioner] type = Steady solve_type = NEWTON - automatic_scaling = true + automatic_scaling = false l_max_its = 1000 nl_max_its = 1000 [] From 56b21fc448dfae8d3b9cfe076fc62a8f497aa2f7 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Thu, 6 Jun 2024 15:03:48 +0100 Subject: [PATCH 33/43] Add Hunt case --- .../step5_MHD_coupling/3d_coupled_hunt.i | 355 ++++++++++++++++++ 1 file changed, 355 insertions(+) create mode 100644 examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled_hunt.i 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..cee09f1 --- /dev/null +++ b/examples/mhd/charge_conserving/development/step5_MHD_coupling/3d_coupled_hunt.i @@ -0,0 +1,355 @@ +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] + active = 'SMP' + # active = 'FSP' + [SMP] + type = SMP + full = true + petsc_options_iname = '-pc_type' + petsc_options_value = ' lu' + [] + [FSP] + type = FSP + topsplit = 'splitting' + [splitting] + splitting = 'velocity pressure currentDensity electricPotential' + splitting_type = additive + [] + [velocity] + vars = 'velocity' + petsc_options_iname = '-pc_type' + petsc_options_value = ' cholesky' + [] + [pressure] + vars = 'pressure' + petsc_options_iname = '-pc_type' + petsc_options_value = ' cholesky' + [] + [currentDensity] + vars = 'currentDensity' + petsc_options_iname = '-pc_type' + petsc_options_value = ' cholesky' + [] + [electricPotential] + vars = 'electricPotential' + petsc_options_iname = '-pc_type' + petsc_options_value = ' cholesky' + [] + [] +[] + +[Executioner] + type = Steady + solve_type = NEWTON + automatic_scaling = false + l_max_its = 1000 + nl_max_its = 1000 +[] + +[Outputs] + exodus = true + execute_on = 'NONLINEAR' +[] From f865e2ae9f9aafcf2d25a7ded44bf82d737133b1 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 21 Aug 2024 15:08:12 +0100 Subject: [PATCH 34/43] Add zero value condition for electricPotential --- .../step3_electrostatic/2d_electrostatic_solenoidal.i | 7 +++++++ .../step3_electrostatic/3d_electrostatic_solenoidal.i | 7 +++++++ 2 files changed, 14 insertions(+) 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 index 4fbf562..41650d5 100644 --- a/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_solenoidal.i +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_solenoidal.i @@ -67,6 +67,13 @@ sigma = 1 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] 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 index 575fc5d..b57c230 100644 --- a/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic_solenoidal.i +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/3d_electrostatic_solenoidal.i @@ -72,6 +72,13 @@ sigma = 1 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] From 9fb60b0f9d5bf0cac989c43e136c0b5c17ce87c1 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 21 Aug 2024 15:08:57 +0100 Subject: [PATCH 35/43] Add zero-divergence check --- .../step5_MHD_coupling/3d_coupled_shercliff.i | 13 +++++++++++++ 1 file changed, 13 insertions(+) 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 index b8c59c1..bcbd9e8 100644 --- 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 @@ -166,6 +166,10 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} expression_y = '20' expression_z = '0' [] + [divergenceFreeCheck] + type = ParsedVectorFunction + div = '0' + [] [] [Variables] @@ -334,6 +338,15 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [] [] +[Postprocessors] + [currentDivergenceFreeError] + type = ElementHDivSemiError + variable = currentDensity + function = divergenceFreeCheck + execute_on = 'NONLINEAR' + [] +[] + [Executioner] type = Steady solve_type = NEWTON From d8ce7cd06e877864340664249b05a87f74c729b6 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Tue, 22 Oct 2024 17:01:11 +0100 Subject: [PATCH 36/43] Increase AD derivative size to 89 (required for charge conservative MHD formulation) --- scripts/install_csd3_cclake.sh | 2 +- scripts/install_csd3_sapphire.sh | 2 +- scripts/install_ubuntu.sh | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) 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..2d36bad 100755 --- a/scripts/install_csd3_sapphire.sh +++ b/scripts/install_csd3_sapphire.sh @@ -52,7 +52,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 # Apply NSFV patch 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 From 4f5d3dae6a1a03e83e94b57cd5e0b4b8aa144391 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 23 Oct 2024 15:32:02 +0100 Subject: [PATCH 37/43] Add jinja2 dependency --- scripts/csd3_python_prep.sh | 1 + 1 file changed, 1 insertion(+) 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 From e91ede44af65771ec919aa9645f7948856651af0 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 23 Oct 2024 15:35:29 +0100 Subject: [PATCH 38/43] Remove NSFVBase.patch (no longer required) --- scripts/NSFVBase.patch | 22 ---------------------- scripts/install_csd3_sapphire.sh | 4 ---- 2 files changed, 26 deletions(-) delete mode 100644 scripts/NSFVBase.patch 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/install_csd3_sapphire.sh b/scripts/install_csd3_sapphire.sh index 2d36bad..50ec23b 100755 --- a/scripts/install_csd3_sapphire.sh +++ b/scripts/install_csd3_sapphire.sh @@ -54,10 +54,6 @@ unset PETSC_DIR PETSC_ARCH ./configure --with-derivative-size=89 -# Apply NSFV patch - -git apply $PROTEUS_DIR/scripts/NSFVBase.patch - cd $PROTEUS_DIR make -j $MOOSE_JOBS From 34fff534d945dc483dfb4952c2a23558127aacd8 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 23 Oct 2024 15:38:57 +0100 Subject: [PATCH 39/43] Add --disable-netgen flag to libMesh --- scripts/install_csd3_sapphire.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/install_csd3_sapphire.sh b/scripts/install_csd3_sapphire.sh index 50ec23b..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 From 725dca25a67052a01395e0331d6facdd11b76137 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 23 Oct 2024 16:21:09 +0100 Subject: [PATCH 40/43] Clarify magneticFieldFunction argument descriptions --- src/kernels/IMHDADCurrentUxB.C | 3 +-- src/kernels/IMHDADMomentumLorentz.C | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/kernels/IMHDADCurrentUxB.C b/src/kernels/IMHDADCurrentUxB.C index 85bc03a..253cd43 100644 --- a/src/kernels/IMHDADCurrentUxB.C +++ b/src/kernels/IMHDADCurrentUxB.C @@ -12,8 +12,7 @@ IMHDADCurrentUxB::validParams() "The Jacobian is computed using automatic differentiation."); params.addRequiredCoupledVar("velocity", "The variable representing the velocity"); params.addParam("magneticFieldFunction", - "A function for the imposed magnetic field that defines a vectorValue method. " - "This cannot be supplied with the component parameters."); + "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( diff --git a/src/kernels/IMHDADMomentumLorentz.C b/src/kernels/IMHDADMomentumLorentz.C index 3bcb05f..d4e46f2 100644 --- a/src/kernels/IMHDADMomentumLorentz.C +++ b/src/kernels/IMHDADMomentumLorentz.C @@ -12,8 +12,7 @@ IMHDADMomentumLorentz::validParams() "The Jacobian is computed using automatic differentiation."); params.addRequiredCoupledVar("currentDensity", "The variable representing the electric current density"); params.addParam("magneticFieldFunction", - "A function for the imposed magnetic field that defines a vectorValue method. " - "This cannot be supplied with the component parameters."); + "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( From 845d449d87a18bba85f1ad201ac23a15374ce4f0 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 23 Oct 2024 18:19:19 +0100 Subject: [PATCH 41/43] Remove FSP preconditioning option; currently not working --- .../step5_MHD_coupling/3d_coupled_hunt.i | 30 ------------------- .../step5_MHD_coupling/3d_coupled_shercliff.i | 30 ------------------- 2 files changed, 60 deletions(-) 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 index cee09f1..607bc69 100644 --- 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 @@ -303,42 +303,12 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [] [Preconditioning] - active = 'SMP' - # active = 'FSP' [SMP] type = SMP full = true petsc_options_iname = '-pc_type' petsc_options_value = ' lu' [] - [FSP] - type = FSP - topsplit = 'splitting' - [splitting] - splitting = 'velocity pressure currentDensity electricPotential' - splitting_type = additive - [] - [velocity] - vars = 'velocity' - petsc_options_iname = '-pc_type' - petsc_options_value = ' cholesky' - [] - [pressure] - vars = 'pressure' - petsc_options_iname = '-pc_type' - petsc_options_value = ' cholesky' - [] - [currentDensity] - vars = 'currentDensity' - petsc_options_iname = '-pc_type' - petsc_options_value = ' cholesky' - [] - [electricPotential] - vars = 'electricPotential' - petsc_options_iname = '-pc_type' - petsc_options_value = ' cholesky' - [] - [] [] [Executioner] 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 index bcbd9e8..c5cacc2 100644 --- 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 @@ -300,42 +300,12 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [] [Preconditioning] - active = 'SMP' - # active = 'FSP' [SMP] type = SMP full = true petsc_options_iname = '-pc_type' petsc_options_value = ' lu' [] - [FSP] - type = FSP - topsplit = 'splitting' - [splitting] - splitting = 'velocity pressure currentDensity electricPotential' - splitting_type = additive - [] - [velocity] - vars = 'velocity' - petsc_options_iname = '-pc_type' - petsc_options_value = ' cholesky' - [] - [pressure] - vars = 'pressure' - petsc_options_iname = '-pc_type' - petsc_options_value = ' cholesky' - [] - [currentDensity] - vars = 'currentDensity' - petsc_options_iname = '-pc_type' - petsc_options_value = ' cholesky' - [] - [electricPotential] - vars = 'electricPotential' - petsc_options_iname = '-pc_type' - petsc_options_value = ' cholesky' - [] - [] [] [Postprocessors] From 9999b1b55a6b2c753eec3b4292c44cc68f5c398f Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 23 Oct 2024 18:20:22 +0100 Subject: [PATCH 42/43] Update case to work with magnetic field changes (AuxVariable to function) --- .../step4_coupledAuxVelocity/3d_coupledAuxVelocity.i | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/examples/mhd/charge_conserving/development/step4_coupledAuxVelocity/3d_coupledAuxVelocity.i b/examples/mhd/charge_conserving/development/step4_coupledAuxVelocity/3d_coupledAuxVelocity.i index 9e95bf9..b72aa3a 100644 --- a/examples/mhd/charge_conserving/development/step4_coupledAuxVelocity/3d_coupledAuxVelocity.i +++ b/examples/mhd/charge_conserving/development/step4_coupledAuxVelocity/3d_coupledAuxVelocity.i @@ -184,10 +184,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} family = LAGRANGE_VEC order = FIRST [] - [magneticField] - family = LAGRANGE_VEC - order = FIRST - [] [] [Kernels] @@ -205,7 +201,7 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} type = IMHDADCurrentUxB variable = currentDensity velocity = velocity - magneticField = magneticField + magneticFieldFunction = magneticFieldFunction [] [divergenceFree] @@ -216,12 +212,6 @@ RATIO_Z_INV = ${fparse 1/RATIO_Z_FWD} [] [AuxKernels] - [magneticFieldKernel] - type = VectorFunctionAux - variable = magneticField - function = magneticFieldFunction - execute_on = INITIAL - [] [velocityKernel] type = VectorFunctionAux variable = velocity From b8bac3981bbed36380c5e96506c041570df08d63 Mon Sep 17 00:00:00 2001 From: Rupert Eardley Date: Wed, 23 Oct 2024 18:21:12 +0100 Subject: [PATCH 43/43] Add comments on issues with running these cases --- .../development/step1_NavierStokes_elementTypes/2d_INSAD.i | 3 +++ .../step3_electrostatic/2d_electrostatic_gradient.i | 4 ++++ 2 files changed, 7 insertions(+) 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 index 9aae768..4bd2d97 100644 --- a/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/2d_INSAD.i +++ b/examples/mhd/charge_conserving/development/step1_NavierStokes_elementTypes/2d_INSAD.i @@ -1,3 +1,6 @@ +# 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 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 index 6067232..75563bc 100644 --- a/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_gradient.i +++ b/examples/mhd/charge_conserving/development/step3_electrostatic/2d_electrostatic_gradient.i @@ -1,3 +1,7 @@ +# 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]