From a177081c1b7906f02d47277c2c5f3e8a5100abfe Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Wed, 7 Jul 2021 01:00:36 -0400 Subject: [PATCH 01/12] Not working: Trying to setup solenoidal projection with full schur. Fails because the preconditioner matrix is assumed to exist, when it doesn't. Also depends on an options path in Full_Projection that isn't compatible with being called from solenoidal projection instead of momentum equation. --- assemble/Full_Projection.F90 | 55 ++- assemble/Makefile.dependencies | 11 +- assemble/Momentum_Equation.F90 | 2 +- assemble/Solenoidal_interpolation.F90 | 299 ++++++++---- femtools/FEFields.F90 | 67 ++- femtools/Makefile.dependencies | 3 +- schemas/adaptivity_options.rnc | 10 +- schemas/adaptivity_options.rng | 10 +- schemas/diagnostic_algorithms.rnc | 22 +- schemas/diagnostic_algorithms.rng | 22 +- schemas/fluidity_options.rnc | 107 +++-- schemas/fluidity_options.rng | 117 +++-- schemas/mesh_options.rnc | 2 +- schemas/mesh_options.rng | 2 +- schemas/prescribed_field_options.rnc | 3 +- schemas/prescribed_field_options.rng | 5 +- schemas/prognostic_field_options.rnc | 16 +- schemas/prognostic_field_options.rng | 16 +- schemas/solvers.rnc | 47 +- schemas/solvers.rng | 43 +- schemas/test_advection_diffusion_options.rnc | 40 +- schemas/test_advection_diffusion_options.rng | 457 +------------------ 22 files changed, 657 insertions(+), 699 deletions(-) diff --git a/assemble/Full_Projection.F90 b/assemble/Full_Projection.F90 index 5eae4e6da9..fd91fab365 100644 --- a/assemble/Full_Projection.F90 +++ b/assemble/Full_Projection.F90 @@ -38,6 +38,7 @@ module Full_Projection use data_structures use sparse_tools use fields + use field_options use petsc_tools use signal_vars use sparse_tools_petsc @@ -62,7 +63,7 @@ module Full_Projection !-------------------------------------------------------------------------------------------------------------------- subroutine petsc_solve_full_projection(x,ctp_m,inner_m,ct_m,rhs,pmat, velocity, & - state, inner_mesh, auxiliary_matrix) + state, inner_mesh, option_path, inner_option_path, auxiliary_matrix) !-------------------------------------------------------------------------------------------------------------------- ! Solve Schur complement problem the nice way (using petsc) !! @@ -81,6 +82,7 @@ subroutine petsc_solve_full_projection(x,ctp_m,inner_m,ct_m,rhs,pmat, velocity, ! state, and inner_mesh are used to setup mg preconditioner of inner solve type(state_type), intent(in):: state type(mesh_type), intent(in):: inner_mesh + character(len=*), optional, intent(in) :: option_path, inner_option_path ! p1-p1 stabilization matrix or free surface terms: type(csr_matrix), optional, intent(in) :: auxiliary_matrix @@ -89,6 +91,7 @@ subroutine petsc_solve_full_projection(x,ctp_m,inner_m,ct_m,rhs,pmat, velocity, Vec y, b ! PETSc solution vector (y), PETSc RHS (b) character(len=OPTION_PATH_LEN) solver_option_path, name + character(len=OPTION_PATH_LEN) l_option_path, l_inner_option_path integer literations type(petsc_numbering_type) petsc_numbering logical lstartfromzero @@ -101,11 +104,25 @@ subroutine petsc_solve_full_projection(x,ctp_m,inner_m,ct_m,rhs,pmat, velocity, ewrite(2,*) 'Ignoring setting from solver option.' lstartfromzero=.true. + if(present(option_path)) then + l_option_path = option_path + else + l_option_path = x%option_path + end if + + if(present(inner_option_path)) then + l_inner_option_path = inner_option_path + else + l_inner_option_path = trim(complete_field_path(l_option_path))//& + "/scheme/use_projection_method"//& + "/full_schur_complement/inner_matrix[0]" + end if + ! Convert matrices to PETSc format, setup PETSc object and PETSc numbering, then ! Build Schur complement and set KSP. ewrite(2,*) 'Entering PETSc setup for Full Projection Solve' call petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering,name,solver_option_path, & - lstartfromzero,inner_m,ctp_m,ct_m,x%option_path,pmat, & + lstartfromzero,inner_m,ctp_m,ct_m,l_option_path,l_inner_option_path,pmat, & rhs, velocity, state, inner_mesh, auxiliary_matrix) ewrite(2,*) 'Create RHS and solution Vectors in PETSc Format' @@ -120,7 +137,7 @@ subroutine petsc_solve_full_projection(x,ctp_m,inner_m,ct_m,rhs,pmat, velocity, end if ewrite(2,*) 'Entering Core PETSc Solve' - ! Solve Ay = b using KSP and PC. Also check convergence. We call this the inner solve. + ! Solve Ay = b using KSP and PC. Also check convergence. call petsc_solve_core(y, A, b, ksp, petsc_numbering, solver_option_path, lstartfromzero, & literations, sfield=x, x0=x%val, nomatrixdump=.true.) @@ -138,7 +155,7 @@ end subroutine petsc_solve_full_projection !-------------------------------------------------------------------------------------------------------- subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,solver_option_path, & - lstartfromzero,inner_m,div_matrix_comp, div_matrix_incomp,option_path,preconditioner_matrix,rhs, & + lstartfromzero,inner_m,div_matrix_comp, div_matrix_incomp,option_path,inner_option_path,preconditioner_matrix,rhs, & velocity, state, inner_mesh, auxiliary_matrix) !-------------------------------------------------------------------------------------------------------- @@ -175,8 +192,8 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so type(csr_matrix), intent(inout) :: preconditioner_matrix ! Stabilization matrix: type(csr_matrix), optional, intent(in) :: auxiliary_matrix - ! Option path: - character(len=*), intent(in):: option_path + ! Option paths: + character(len=*), intent(in):: option_path, inner_option_path type(vector_field), intent(in) :: velocity ! used to retrieve strong diricihlet bcs ! state, and inner_mesh are used to setup mg preconditioner of inner solve type(state_type), intent(in):: state @@ -203,39 +220,36 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so Mat S ! PETSc Stabilization matrix (auxiliary_matrix) Mat pmat ! PETSc preconditioning matrix - character(len=OPTION_PATH_LEN) :: inner_option_path, inner_solver_option_path + character(len=OPTION_PATH_LEN) :: inner_solver_option_path integer, dimension(:,:), pointer :: save_gnn2unn type(integer_set), dimension(velocity%dim):: boundary_row_set - integer reference_node, i, rotation_stat + integer reference_node, i, rotation_stat, ref_stat logical parallel, have_auxiliary_matrix, have_preconditioner_matrix logical :: apply_reference_node, apply_reference_node_from_coordinates, reference_node_owned ! Sort option paths etc... solver_option_path=complete_solver_option_path(option_path) - inner_option_path= trim(option_path)//& - "/prognostic/scheme/use_projection_method& - &/full_schur_complement/inner_matrix[0]" if (have_option(trim(option_path)//'/name')) then call get_option(trim(option_path)//'/name', name) - ewrite(1,*) 'Inside petsc_solve_(block_)csr, solving for: ', trim(name) + ewrite(1,*) 'Inside petsc_solve_full_projection, solving for: ', trim(name) else - ewrite(1,*) 'Inside petsc_solve_(block_)csr, solving using option_path: ', trim(option_path) + ewrite(1,*) 'Inside petsc_solve_full_projection, solving using option_path: ', trim(option_path) name=option_path end if - ! Are we applying a reference pressure node? - apply_reference_node = have_option(trim(option_path)//& - '/prognostic/reference_node') - apply_reference_node_from_coordinates = have_option(trim(option_path)//& - '/prognostic/reference_coordinates') + ! Are we applying a reference node? + apply_reference_node = have_option(trim(complete_field_path(option_path, stat=ref_stat))//& + &"/reference_node") + apply_reference_node_from_coordinates = have_option(trim(complete_field_path(option_path, stat=ref_stat))//& + &"/reference_coordinates") ! If so, impose reference pressure node: if(apply_reference_node) then - call get_option(trim(option_path)//& - '/prognostic/reference_node', reference_node) + call get_option(trim(complete_field_path(option_path, stat=ref_stat))//& + '/reference_node', reference_node) if (GetProcNo()==1) then ewrite(2,*) 'Imposing_reference_pressure_node' allocate(ghost_nodes(1:1)) @@ -425,6 +439,7 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so inner_solver_option_path, petsc_numbering=petsc_numbering_u, startfromzero_in=.true.) ! leaving out petsc_numbering and mesh, so "iteration_vtus" monitor won't work! + ! FIXME: broken logic here, shouldn't depend on option path ! Assemble preconditioner matrix in petsc format (if required): have_preconditioner_matrix=.not.(have_option(trim(option_path)//& "/prognostic/scheme/use_projection_method& diff --git a/assemble/Makefile.dependencies b/assemble/Makefile.dependencies index d974dc9f71..5a12ff33f3 100644 --- a/assemble/Makefile.dependencies +++ b/assemble/Makefile.dependencies @@ -732,11 +732,12 @@ Solenoidal_interpolation.o ../include/solenoidal_interpolation_module.mod: \ ../include/divergence_matrix_cg.mod ../include/divergence_matrix_cv.mod \ ../include/element_numbering.mod ../include/fdebug.h \ ../include/fefields.mod ../include/fetools.mod ../include/field_options.mod \ - ../include/fields.mod ../include/fldebug.mod ../include/futils.mod \ - ../include/global_parameters.mod ../include/interpolation_module.mod \ - ../include/linked_lists.mod ../include/momentum_cg.mod \ - ../include/momentum_dg.mod ../include/solvers.mod \ - ../include/sparse_matrices_fields.mod ../include/sparse_tools.mod \ + ../include/fields.mod ../include/fldebug.mod ../include/full_projection.mod \ + ../include/futils.mod ../include/global_parameters.mod \ + ../include/interpolation_module.mod ../include/linked_lists.mod \ + ../include/momentum_cg.mod ../include/momentum_dg.mod \ + ../include/solvers.mod ../include/sparse_matrices_fields.mod \ + ../include/sparse_tools.mod ../include/sparse_tools_petsc.mod \ ../include/sparsity_patterns.mod ../include/state_module.mod \ ../include/supermesh_construction.mod ../include/tensors.mod \ ../include/transform_elements.mod ../include/vector_tools.mod diff --git a/assemble/Momentum_Equation.F90 b/assemble/Momentum_Equation.F90 index ce5689a698..71a68e4f0e 100644 --- a/assemble/Momentum_Equation.F90 +++ b/assemble/Momentum_Equation.F90 @@ -493,7 +493,7 @@ subroutine solve_momentum(state, at_first_timestep, timestep) full_projection_preconditioner => cmc_m case default ! Developer error... out of sync options input and code - FLAbort("Unknown Matrix Type for Full_Projection") + FLAbort("Unknown preconditioner matrix for full schur complement in momentum equation.") end select ! Decide on configuration of inner_m for full_projection solve: diff --git a/assemble/Solenoidal_interpolation.F90 b/assemble/Solenoidal_interpolation.F90 index bbbf4dfcdd..0c11003535 100644 --- a/assemble/Solenoidal_interpolation.F90 +++ b/assemble/Solenoidal_interpolation.F90 @@ -9,6 +9,7 @@ module solenoidal_interpolation_module use futils use spud use sparse_tools + use sparse_tools_petsc use vector_tools use tensors use element_numbering, only: FAMILY_SIMPLEX @@ -24,14 +25,16 @@ module solenoidal_interpolation_module use interpolation_module use sparse_matrices_fields use solvers + use full_projection use fefields use dgtools use assemble_cmc, only: assemble_cmc_dg, repair_stiff_nodes,& - zero_stiff_nodes, assemble_masslumped_cmc + zero_stiff_nodes, assemble_masslumped_cmc, assemble_diagonal_schur use boundary_conditions_from_options use divergence_matrix_cv, only: assemble_divergence_matrix_cv use divergence_matrix_cg, only: assemble_divergence_matrix_cg - use momentum_cg, only: correct_masslumped_velocity, add_kmk_matrix, add_kmk_rhs, assemble_kmk_matrix + use momentum_cg, only: correct_velocity_cg, correct_masslumped_velocity, & + add_kmk_matrix, add_kmk_rhs, assemble_kmk_matrix use momentum_dg, only: correct_velocity_dg implicit none @@ -95,15 +98,24 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & type(scalar_field), pointer :: s_field - logical :: dg, lump_mass, lump_on_submesh, div_cv, div_cg, apply_kmk + logical :: dg, discontinuous, div_cv, div_cg, apply_kmk + logical :: assemble_cmc, assemble_lagrange_mass, assemble_cdiagmc + logical :: assemble_lumped_mass, lump_mass_on_submesh + logical :: full_schur, assemble_schur_aux integer :: dim, j real :: dt + character(len=FIELD_NAME_LEN) :: pressure_pmat type(block_csr_matrix), pointer :: ct_m type(block_csr_matrix), pointer :: ctp_m type(scalar_field) :: ct_rhs, kmk_rhs type(csr_sparsity) :: ct_m_sparsity, cmc_m_sparsity - type(csr_matrix) :: cmc_m + type(csr_sparsity) :: field_mass_sparsity, lagrange_mass_sparsity + type(csr_matrix), target :: cmc_m, cdiagmc_m + type(csr_matrix) :: schur_aux + type(csr_matrix), target :: lagrange_mass + type(petsc_csr_matrix), target :: field_mass + type(csr_matrix), pointer :: pschur type(block_csr_matrix) :: inverse_field_mass type(scalar_field) :: field_lumped_mass type(vector_field) :: inverse_field_lumped_mass_vector @@ -137,14 +149,29 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & end if dg = (continuity(v_field) < 0) + discontinuous = have_option(trim(l_option_path)//& + "/interpolated_field/discontinuous") + + if (dg) then + assert(discontinuous) + end if - lump_mass = have_option(trim(l_option_path)//& + assemble_lumped_mass = have_option(trim(l_option_path)//& "/interpolated_field/discontinuous/lump_mass_matrix") & - .or. .not.dg + .or. & + have_option(trim(l_option_path)//& + "/interpolated_field/continuous/lump_mass_matrix") & + .or. & + have_option(trim(l_option_path)//& + "/interpolated_field/continuous/full_schur_complement"//& + "/preconditioner_matrix::LumpedSchurComplement") - lump_on_submesh = have_option(trim(l_option_path)//& + lump_mass_on_submesh = have_option(trim(l_option_path)//& "/interpolated_field/continuous/lump_mass_matrix/use_submesh") & - .and. .not.dg + .or. & + have_option(trim(l_option_path)//& + "/interpolated_field/continuous/full_schur_complement"//& + "/preconditioner_matrix[0]/lump_on_submesh") div_cv = have_option(trim(l_option_path) //& &"/lagrange_multiplier/spatial_discretisation/control_volumes") @@ -164,6 +191,48 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & & "/lagrange_multiplier/spatial_discretisation/continuous_galerkin/remove_stabilisation_term") .and. & & .not. div_cv) + full_schur = have_option(trim(l_option_path)//& + "/interpolated_field/continuous/full_schur_complement") + + assemble_cmc = .true. ! true, unless we're using full_schur and even then, it's complicated... + assemble_cdiagmc = .false. ! only used as a possible preconditioner for full_schur + assemble_lagrange_mass = .false. ! only used as a possible preconditioner for full_schur + assemble_schur_aux = .false. + if(full_schur) then + assemble_schur_aux = apply_kmk + ! Check to see whether pressure cmc_m preconditioning matrix is needed: + call get_option(trim(l_option_path)//& + "/interpolated_field/continuous/full_schur_complement"//& + "/preconditioner_matrix[0]/name", pressure_pmat) + + select case(pressure_pmat) + case("LumpedSchurComplement") + pschur => cmc_m + case("DiagonalSchurComplement") + assemble_cmc = .false. + assemble_cdiagmc = .true. + pschur => cdiagmc_m + case("LagrangeMassMatrix") + assemble_cmc = .false. + assemble_lagrange_mass = .true. + pschur => lagrange_mass + case("NoPreconditionerMatrix") + assemble_cmc = .false. + nullify(pschur) + case default + ! Developer error... out of sync options input and code + FLAbort("Unknown preconditioner matrix for full schur complement in solenoidal interpolation.") + end select + + ! always need a field_mass matrix if doing full_schur + field_mass_sparsity = make_sparsity(v_field%mesh, v_field%mesh, "FieldMassSparsity") + call allocate(field_mass, field_mass_sparsity, (/v_field%dim, v_field%dim/), & + group_size=(/v_field%dim, v_field%dim/), & + diagonal=.true., name="FieldMassMatrix") + call zero(field_mass) + + end if + ct_m_sparsity = make_sparsity(lagrange_mesh, v_field%mesh, "DivergenceSparsity") allocate(ct_m) call allocate(ct_m, ct_m_sparsity, blocks=(/1, dim/), name="DivergenceMatrix") @@ -181,11 +250,32 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & call allocate(ct_rhs, lagrange_mesh, "DivergenceRHS") call zero(ct_rhs) + + if (assemble_cmc.or.assemble_schur_aux.or.assemble_cdiagmc) then + cmc_m_sparsity = make_sparsity_transpose(lagrange_mesh, v_field%mesh, "LagrangeProjectionSparsity") + end if - cmc_m_sparsity = make_sparsity_transpose(lagrange_mesh, v_field%mesh, "LagrangeProjectionSparsity") - call allocate(cmc_m, cmc_m_sparsity, name="LagrangeProjectionMatrix") - call zero(cmc_m) + if (assemble_cmc) then + call allocate(cmc_m, cmc_m_sparsity, name="LagrangeProjectionMatrix") + call zero(cmc_m) + end if + + if (assemble_schur_aux) then + call allocate(schur_aux, cmc_m_sparsity, name="SchurAuxilliaryMatrix") + call zero(schur_aux) + end if + if (assemble_cdiagmc) then + call allocate(cdiagmc_m, cmc_m_sparsity, name="DiagonalSchurMatrix") + call zero(cdiagmc_m) + end if + + if (assemble_lagrange_mass) then + lagrange_mass_sparsity = make_sparsity(lagrange_mesh, lagrange_mesh, "LagrangeMassSparsity") + call allocate(lagrange_mass, lagrange_mass_sparsity, name="LagrangeMassMatrix") + call zero(lagrange_mass) + end if + call allocate(projec_rhs, lagrange_mesh, "LagrangeProjectionRHS") call zero(projec_rhs) @@ -198,64 +288,49 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & call zero(kmk_rhs) end if - if(lump_mass) then + if(assemble_lumped_mass) then call allocate(field_lumped_mass, v_field%mesh, "FieldLumpedMass") call zero(field_lumped_mass) call allocate(inverse_field_lumped_mass_vector, dim, v_field%mesh, & "InverseFieldLumpedMassVector") - else if (.not. dg) then - FLExit("Not possible to not lump the mass if not dg.") end if - if(div_cg) then - if(lump_mass) then - if(lump_on_submesh) then - call assemble_divergence_matrix_cg(ct_m, local_state, ct_rhs=ct_rhs, & - test_mesh=lagrange_mesh, field=v_field, & - option_path = trim(l_option_path)//"/lagrange_multiplier") - ! now get the mass matrix lumped on the submesh - call compute_lumped_mass_on_submesh(local_state, field_lumped_mass) - - else - call assemble_divergence_matrix_cg(ct_m, local_state, ct_rhs=ct_rhs, & - test_mesh=lagrange_mesh, field=v_field, & - option_path = trim(l_option_path)//"/lagrange_multiplier", & - grad_mass_lumped = field_lumped_mass) - end if - elseif(dg) then - call assemble_divergence_matrix_cg(ct_m, local_state, ct_rhs=ct_rhs, & - test_mesh=lagrange_mesh, field=v_field, & - option_path = trim(l_option_path)//"/lagrange_multiplier") - - ! now get the dg inverse mass matrix - call construct_inverse_mass_matrix_dg(inverse_field_mass, v_field, coordinate) - + if(full_schur) then + if(assemble_lagrange_mass) then + call compute_mass(coordinate, lagrange_mesh, lagrange_mass) + end if + if (assemble_lumped_mass.and.(.not.lump_mass_on_submesh)) then + call compute_mass(coordinate, v_field%mesh, field_mass, field_lumped_mass) else - FLExit("Not possible to not lump the mass if not dg.") + call compute_mass(coordinate, v_field%mesh, field_mass) + end if + else + if(assemble_lumped_mass.and.(.not.lump_mass_on_submesh)) then + call compute_lumped_mass(coordinate, field_lumped_mass) end if + end if + + if(assemble_lumped_mass.and.lump_mass_on_submesh) then + ! get the mass matrix lumped on the submesh + call compute_lumped_mass_on_submesh(local_state, field_lumped_mass) + else if (discontinuous) then + ! get the inverse discontinuous mass matrix + call construct_inverse_mass_matrix_dg(inverse_field_mass, v_field, coordinate) + end if + if(div_cg) then + call assemble_divergence_matrix_cg(ct_m, local_state, ct_rhs=ct_rhs, & + test_mesh=lagrange_mesh, field=v_field, & + option_path = trim(l_option_path)//"/lagrange_multiplier") + ! If CG lagrange with CV tested divergence then form the other C matrix. ! This will overwrite the ct_rhs formed above. if (cg_lagrange_cv_test_divergence) then - call assemble_divergence_matrix_cv(ctp_m, local_state, ct_rhs=ct_rhs, & test_mesh=lagrange_mesh, field=v_field) - end if - elseif(div_cv) then - if(lump_mass) then - if(lump_on_submesh) then - call compute_lumped_mass_on_submesh(local_state, field_lumped_mass) - else - call compute_lumped_mass(coordinate, field_lumped_mass) - end if - else if(dg) then - call construct_inverse_mass_matrix_dg(inverse_field_mass, v_field, coordinate) - else - FLExit("Not possible to not lump the mass if not dg.") - end if - + else if(div_cv) then call assemble_divergence_matrix_cv(ct_m, local_state, ct_rhs=ct_rhs, & test_mesh=lagrange_mesh, field=v_field) @@ -264,24 +339,30 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & FLAbort("Unknown spatial discretisation option for the lagrange multiplier.") end if - if(lump_mass) then - call invert(field_lumped_mass) - - do j=1, inverse_field_lumped_mass_vector%dim - call set(inverse_field_lumped_mass_vector, j, field_lumped_mass) - end do - - call apply_dirichlet_conditions_inverse_mass(inverse_field_lumped_mass_vector, v_field) + if(assemble_cmc) then + if(assemble_lumped_mass) then + call invert(field_lumped_mass) + + do j=1, inverse_field_lumped_mass_vector%dim + call set(inverse_field_lumped_mass_vector, j, field_lumped_mass) + end do + + call apply_dirichlet_conditions_inverse_mass(inverse_field_lumped_mass_vector, v_field) - call assemble_masslumped_cmc(cmc_m, ctp_m, inverse_field_lumped_mass_vector, ct_m) - else if(dg) then - call assemble_cmc_dg(cmc_m, ctp_m, ct_m, inverse_field_mass) - else - FLExit("Not possible to not lump the mass if not dg.") + call assemble_masslumped_cmc(cmc_m, ctp_m, inverse_field_lumped_mass_vector, ct_m) + else if(dg) then + call assemble_cmc_dg(cmc_m, ctp_m, ct_m, inverse_field_mass) + else + FLExit("Not possible to not lump the mass if not dg or full_schur.") + end if + + if(stiff_nodes_repair) then + call repair_stiff_nodes(cmc_m, stiff_nodes_list) + end if end if - if(stiff_nodes_repair) then - call repair_stiff_nodes(cmc_m, stiff_nodes_list) + if(assemble_cdiagmc) then + call assemble_diagonal_schur(cdiagmc_m, v_field, field_mass, ctp_m, ct_m) end if call mult(projec_rhs, ctp_m, v_field) @@ -291,13 +372,19 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & ! construction of the sparsity call insert(local_state, lagrange, name="Pressure") call insert(local_state, v_field, name="Velocity") - call insert(local_state, lagrange, name="Pressure") - call insert(local_state, v_field, name="Velocity") ! end of hack... you can look again now ewrite(2,*) "Assembling P1-P1 stabilisation" call assemble_kmk_matrix(local_state, lagrange_mesh, coordinate, theta_pg=1.0) - call add_kmk_matrix(local_state, cmc_m) + if(assemble_cmc) then + call add_kmk_matrix(local_state, cmc_m) + end if + if(assemble_schur_aux) then + call add_kmk_matrix(local_state, schur_aux) + end if + if(assemble_cdiagmc) then + call add_kmk_matrix(local_state, cdiagmc_m) + end if if(associated(s_field)) then ! Should the timestep be passed in here? not sure @@ -315,48 +402,84 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & call scale(projec_rhs, -1.0) call addto(projec_rhs, ct_rhs) - call impose_reference_pressure_node(cmc_m, projec_rhs, coordinate, trim(l_option_path)//"/lagrange_multiplier") + if (assemble_cmc) then + call impose_reference_pressure_node(cmc_m, projec_rhs, coordinate, trim(l_option_path)//"/lagrange_multiplier") - if(stiff_nodes_repair) then - call zero_stiff_nodes(projec_rhs, stiff_nodes_list) + ! only have a stiff_nodes_list if we're assembling cmc + if(stiff_nodes_repair) then + call zero_stiff_nodes(projec_rhs, stiff_nodes_list) + end if end if - call petsc_solve(lagrange, cmc_m, projec_rhs) + if (full_schur) then + if (apply_kmk) then + call petsc_solve_full_projection(lagrange, ctp_m, field_mass, ct_m, & + projec_rhs, pschur, v_field, & + local_state, v_field%mesh, & + option_path=trim(l_option_path)//"/lagrange_multiplier",& + inner_option_path=trim(l_option_path)//& + "/interpolated_field/continuous/full_schur_complement/inner_matrix[0]",& + auxiliary_matrix=schur_aux) + else + call petsc_solve_full_projection(lagrange, ctp_m, field_mass, ct_m, & + projec_rhs, pschur, v_field, & + local_state, v_field%mesh, & + option_path=trim(l_option_path)//"/lagrange_multiplier",& + inner_option_path=trim(l_option_path)//& + "/interpolated_field/continuous/full_schur_complement/inner_matrix[0]") + end if + else + call petsc_solve(lagrange, cmc_m, projec_rhs) + end if if(associated(s_field)) then call addto(s_field, lagrange) end if - if(lump_mass) then + if(full_schur) then + call correct_velocity_cg(v_field, field_mass, ct_m, lagrange, local_state) + else if(assemble_lumped_mass) then call correct_masslumped_velocity(v_field, inverse_field_lumped_mass_vector, ct_m, lagrange) - else if(dg) then + else if(discontinuous) then call correct_velocity_dg(v_field, inverse_field_mass, ct_m, lagrange) - else - FLExit("Not possible to not lump the mass if not dg.") end if + if (full_schur) then + call deallocate(field_mass) + end if call deallocate(ct_m_sparsity) call deallocate(ct_m) deallocate(ct_m) - call deallocate(ct_rhs) - call deallocate(cmc_m_sparsity) - call deallocate(cmc_m) - call deallocate(projec_rhs) - call deallocate(lagrange) if (cg_lagrange_cv_test_divergence) then call deallocate(ctp_m) deallocate(ctp_m) end if + call deallocate(ct_rhs) + if (assemble_cmc.or.assemble_schur_aux.or.assemble_cdiagmc) then + call deallocate(cmc_m_sparsity) + end if + if (assemble_cmc) then + call deallocate(cmc_m) + end if + if (assemble_cdiagmc) then + call deallocate(cdiagmc_m) + end if + if (assemble_schur_aux) then + call deallocate(schur_aux) + end if + if (assemble_lagrange_mass) then + call deallocate(lagrange_mass) + end if + call deallocate(projec_rhs) + call deallocate(lagrange) if(apply_kmk) then call deallocate(kmk_rhs) end if - if(lump_mass) then + if(assemble_lumped_mass) then call deallocate(field_lumped_mass) call deallocate(inverse_field_lumped_mass_vector) - else if(dg) then + else if(discontinuous) then call deallocate(inverse_field_mass) - else - FLExit("Not possible to not lump the mass if not dg.") end if call deallocate(local_state) diff --git a/femtools/FEFields.F90 b/femtools/FEFields.F90 index 36fc824f5f..8a8c5913a8 100644 --- a/femtools/FEFields.F90 +++ b/femtools/FEFields.F90 @@ -8,6 +8,7 @@ module fefields use elements, only: element_type use parallel_tools use sparse_tools + use sparse_tools_petsc use transform_elements, only: transform_to_physical, element_volume use fetools, only: shape_shape, shape_rhs, shape_vector_rhs use fields @@ -25,6 +26,10 @@ module fefields module procedure project_scalar_field, project_vector_field end interface + interface compute_mass + module procedure compute_mass_csr, compute_mass_petsc_csr + end interface + private public :: compute_lumped_mass, compute_mass, compute_projection_matrix, add_source_to_rhs, & compute_lumped_mass_on_submesh, compute_cv_mass, project_field @@ -211,7 +216,7 @@ subroutine compute_lumped_mass(positions, lumped_mass, density, vfrac) end subroutine compute_lumped_mass - subroutine compute_mass(positions, mesh, mass, lumped_mass, density) + subroutine compute_mass_csr(positions, mesh, mass, lumped_mass, density) type(vector_field), intent(in) :: positions type(mesh_type), intent(in) :: mesh type(csr_matrix), intent(inout) :: mass @@ -226,7 +231,7 @@ subroutine compute_mass(positions, mesh, mass, lumped_mass, density) real, dimension(ele_ngi(mesh, 1)) :: density_gi - ewrite(1,*) 'In compute_mass' + ewrite(1,*) 'In compute_mass_csr' if(present(density)) then l_density => density @@ -258,7 +263,63 @@ subroutine compute_mass(positions, mesh, mass, lumped_mass, density) deallocate(l_density) end if - end subroutine compute_mass + end subroutine compute_mass_csr + + subroutine compute_mass_petsc_csr(positions, mesh, mass, lumped_mass, density) + type(vector_field), intent(in) :: positions + type(mesh_type), intent(in) :: mesh + type(petsc_csr_matrix), intent(inout) :: mass + type(scalar_field), intent(inout), optional :: lumped_mass + type(scalar_field), intent(inout), target, optional :: density + + integer :: ele, dimi + real, dimension(ele_ngi(mesh, 1)) :: detwei + type(element_type), pointer :: t_shape + real, dimension(ele_loc(mesh, 1), ele_loc(mesh, 1)) :: mass_matrix + type(scalar_field), pointer :: l_density + + real, dimension(ele_ngi(mesh, 1)) :: density_gi + + ewrite(1,*) 'In compute_mass_petsc_csr' + + ! check our assumption below of a square matrix is valid + if(blocks(mass,1)/=blocks(mass,2)) then + FLAbort("compute_mass_petsc_csr assumes block square matrices") + end if + + if(present(density)) then + l_density => density + else + allocate(l_density) + call allocate(l_density, mesh, name="LocalDensity", field_type=FIELD_TYPE_CONSTANT) + call set(l_density, 1.0) + end if + + call zero(mass) + if(present(lumped_mass)) then + assert(lumped_mass%mesh==mesh) + call zero(lumped_mass) + end if + + do ele=1,ele_count(mesh) + t_shape => ele_shape(mesh, ele) + density_gi = ele_val_at_quad(l_density, ele) + call transform_to_physical(positions, ele, detwei=detwei) + mass_matrix = shape_shape(t_shape, t_shape, detwei*density_gi) + do dimi=1,blocks(mass,1) + call addto(mass, dimi, dimi, ele_nodes(mesh, ele), ele_nodes(mesh,ele), mass_matrix) + end do + if(present(lumped_mass)) then + call addto(lumped_mass, ele_nodes(lumped_mass, ele), sum(mass_matrix, 2)) + end if + end do + + if(.not.present(density)) then + call deallocate(l_density) + deallocate(l_density) + end if + + end subroutine compute_mass_petsc_csr subroutine compute_lumped_mass_on_submesh(state, lumped_mass, density, vfrac) diff --git a/femtools/Makefile.dependencies b/femtools/Makefile.dependencies index a47d302e80..e30cb40bbc 100644 --- a/femtools/Makefile.dependencies +++ b/femtools/Makefile.dependencies @@ -373,7 +373,8 @@ FEFields.o ../include/fefields.mod: FEFields.F90 \ ../include/field_options.mod ../include/fields.mod ../include/fldebug.mod \ ../include/halos.mod ../include/parallel_tools.mod \ ../include/sparse_matrices_fields.mod ../include/sparse_tools.mod \ - ../include/state_module.mod ../include/transform_elements.mod + ../include/sparse_tools_petsc.mod ../include/state_module.mod \ + ../include/transform_elements.mod ../include/fetools.mod: FETools.o @true diff --git a/schemas/adaptivity_options.rnc b/schemas/adaptivity_options.rnc index 3a94eb7d80..b90d5ad027 100644 --- a/schemas/adaptivity_options.rnc +++ b/schemas/adaptivity_options.rnc @@ -15,7 +15,7 @@ adaptivity_preprocessing = real }?, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -1450,7 +1450,7 @@ interpolation_algorithm_vector = }?, ## Solver options for the conservative potential calculation element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, ( ## Galerkin projection. By default, conservative, non-dissipative and @@ -1520,7 +1520,7 @@ interpolation_algorithm_vector = ), ## Solver options for the decomposition element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, comment }| @@ -1583,7 +1583,7 @@ interpolation_algorithm_vector = ## Solver options for the geopressure conservative potential ## calculation element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, ( ## Galerkin projection. By default, conservative, non-dissipative and @@ -1753,7 +1753,7 @@ continuous_projection = ## conservation properties are affected by the tolerance of the ## linear solve. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }| ## Lump the mass matrix on the left hand side of the galerkin projection. ## Hence solver options aren't necessary. diff --git a/schemas/adaptivity_options.rng b/schemas/adaptivity_options.rng index cc7fa3be3a..c4305c2645 100644 --- a/schemas/adaptivity_options.rng +++ b/schemas/adaptivity_options.rng @@ -18,7 +18,7 @@ The length scale is alpha*local element width (a scalar) - + @@ -1717,7 +1717,7 @@ Automatic when not using P1P1. Solver options for the conservative potential calculation - + @@ -1801,7 +1801,7 @@ residual Solver options for the decomposition - + @@ -1872,7 +1872,7 @@ functions. Solver options for the geopressure conservative potential calculation - + @@ -2078,7 +2078,7 @@ Defaults to computer precision if unspecified. This method requires the inversion of a mass matrix. Note that conservation properties are affected by the tolerance of the linear solve. - + Lump the mass matrix on the left hand side of the galerkin projection. diff --git a/schemas/diagnostic_algorithms.rnc b/schemas/diagnostic_algorithms.rnc index fd5eb0d198..3437581b3c 100644 --- a/schemas/diagnostic_algorithms.rnc +++ b/schemas/diagnostic_algorithms.rnc @@ -279,7 +279,7 @@ scalar_galerkin_projection_algorithm = scalar_source_field, ## Solver options. Required if projecting onto a continuous mesh. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }? } ) @@ -292,7 +292,7 @@ vector_galerkin_projection_algorithm = vector_source_field, ## Solver options. Required if projecting onto a continuous mesh. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }? } ) @@ -312,7 +312,7 @@ helmholtz_smoothed_scalar_algorithm = }, ## Solver options. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -332,7 +332,7 @@ helmholtz_smoothed_vector_algorithm = }, ## Solver options. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -352,7 +352,7 @@ helmholtz_smoothed_tensor_algorithm = }, ## Solver options. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -372,7 +372,7 @@ helmholtz_anisotropic_smoothed_scalar_algorithm = }, ## Solver options. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -392,7 +392,7 @@ helmholtz_anisotropic_smoothed_vector_algorithm = }, ## Solver options. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -412,7 +412,7 @@ helmholtz_anisotropic_smoothed_tensor_algorithm = }, ## Solver options. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -589,7 +589,7 @@ finite_element_divergence_algorithm = ## Consistent mass Galerkin projection of divergence. Requires ## solver options. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }| ## Lumped mass Galerkin projection of divergence element lump_mass { @@ -832,7 +832,7 @@ scalar_potential_algorithm = }?, ## Linear solver options. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -886,7 +886,7 @@ projection_scalar_potential_algorithm = }?, ## Linear solver options. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) diff --git a/schemas/diagnostic_algorithms.rng b/schemas/diagnostic_algorithms.rng index c8446e2de8..0f83edd825 100644 --- a/schemas/diagnostic_algorithms.rng +++ b/schemas/diagnostic_algorithms.rng @@ -390,7 +390,7 @@ this option to these these internal algorithms. Solver options. Required if projecting onto a continuous mesh. - + @@ -408,7 +408,7 @@ this option to these these internal algorithms. Solver options. Required if projecting onto a continuous mesh. - + @@ -432,7 +432,7 @@ The length scale is alpha*local element width (a scalar) Solver options. - + @@ -455,7 +455,7 @@ The length scale is alpha*local element width (a scalar) Solver options. - + @@ -478,7 +478,7 @@ The length scale is alpha*local element width (a scalar) Solver options. - + @@ -501,7 +501,7 @@ The length scale is alpha*local element size (a tensor). Solver options. - + @@ -524,7 +524,7 @@ The length scale is alpha*local element size (a tensor). Solver options. - + @@ -547,7 +547,7 @@ The length scale is alpha*local element size (a tensor). Solver options. - + @@ -788,7 +788,7 @@ the finite element C^T matrix. Consistent mass Galerkin projection of divergence. Requires solver options. - + Lumped mass Galerkin projection of divergence @@ -1117,7 +1117,7 @@ all boundaries. Linear solver options. - + @@ -1190,7 +1190,7 @@ Automatic when not using P1P1. Linear solver options. - + diff --git a/schemas/fluidity_options.rnc b/schemas/fluidity_options.rnc index 829225691b..44c1c77451 100644 --- a/schemas/fluidity_options.rnc +++ b/schemas/fluidity_options.rnc @@ -2970,7 +2970,7 @@ diagnostic_cv_gradient_vector_field = }?, ## Solver options are necessary if you're not lumping your mass or if you're field isn't dg element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }?, ## Normalise the gradient by its magnitude element normalise { @@ -2990,7 +2990,7 @@ diagnostic_gradient_vector_field = ( ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, diagnostic_output_options, diagnostic_vector_stat_options, @@ -3019,7 +3019,7 @@ diagnostic_fe_divergence_scalar_field = ( ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, diagnostic_output_options, diagnostic_scalar_stat_options, @@ -4299,7 +4299,7 @@ scalar_field_choice = )+, ## Solver - only necessary in combination with an explicit no_normal_stress free_surface bc under Velocity. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, prognostic_scalar_output_options, prognostic_scalar_stat_options, @@ -5479,7 +5479,7 @@ scalar_field_choice = velocity_mesh_choice, ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, diagnostic_scalar_field }| @@ -5585,7 +5585,7 @@ scalar_field_choice = empty }?, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }?, diagnostic_scalar_field }| @@ -5723,7 +5723,7 @@ scalar_field_choice = }?, diagnostic_scalar_field_no_adapt, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -5743,7 +5743,7 @@ scalar_field_choice = internal_algorithm, diagnostic_scalar_field_no_adapt, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -5950,7 +5950,7 @@ vector_field_choice = internal_algorithm, mesh_choice, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, diagnostic_vector_field }| @@ -6197,7 +6197,7 @@ vector_field_choice = empty }?, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }?, diagnostic_vector_field }| @@ -7499,18 +7499,60 @@ solenoidal_options = element interpolated_field { ( element continuous { - ## Lump the mass matrix for the assembly of the projection matrix - ## (not for the initial galerkin projection) - ## - ## Required when using interpolating continuous fields - element lump_mass_matrix { - ## Lump on the submesh. - ## This only works for simplex meshes and is only - ## strictly valid on 2d meshes. - element use_submesh { - empty - }? - } + ( + ## Solver options for the (inner) solve. + element full_schur_complement { + ## Specify the inner matrix (IM) to form the projection schur complement (C^T*IM^{-1}*C). + ## Use the full mass matrix (the only option). + element inner_matrix { + attribute name { "FullMassMatrix" }, + ## Solver options to solve the inner (mass) matrix in the schur complement solve, also used + ## to solve the correction solve afterwards. + element solver { + linear_solver_options_sym_scalar + } + }, + ( + ## Specify the preconditioner matrix to use on the schur complement. + ## + ## Assemble a Schur Complement using a lumped mass inner matrix. + element preconditioner_matrix { + attribute name { "LumpedSchurComplement" }, + element lump_on_submesh { + empty + }? + }| + ## Specify the preconditioner matrix to use on the schur complement. + ## + ## DiagonalSchurComplement = C_P^T * [(M)_diagonal]^-1 * C + element preconditioner_matrix { + attribute name { "DiagonalSchurComplement" }, + empty + }| + ## Specify the preconditioner matrix to use on the schur complement. + ## + ## Lagrange Mass Matrix. + element preconditioner_matrix { + attribute name { "LagrangeMassMatrix" }, + empty + }| + ## Specify the preconditioner matrix to use on the schur complement. + element preconditioner_matrix { + attribute name { "NoPreconditionerMatrix" }, + empty + } + ) + }| + ## Lump the mass matrix for the assembly of the projection matrix + element lump_mass_matrix { + ## Lump on the submesh. + ## This only works for simplex meshes and is only + ## strictly valid on 2d meshes. + element use_submesh { + empty + }? + } + ) }| element discontinuous { ## Lump the mass matrix for the assembly of the projection matrix @@ -7554,9 +7596,20 @@ solenoidal_options = } ) }, - element reference_node { - integer - }?, + ( + ## Reference node (node at which lagrange multiplier = 0). Note that the node number must be less than the total number of nodes. + ## If running in parallel, the node number must be less than the total number of nodes of the first processor. + ## ** Note - it is also an option to remove the null-space of the residual vector. This + ## option is available under solvers + element reference_node { + integer + }| + ## Input coordinates of desired reference node. If a node does not exist at these + ## coordinates, the nearest vertex will be selected. + element reference_coordinates { + real_dim_vector + } + )?, ## **UNDER DEVELOPMENT** ## This searches the CMC matrix diagonal looking for nodes that are less than the maximum value time epsilon(0.0) (i.e. nodes that are effectively zero). ## It then zeros that row and column and places a one on the diagonal and a zero on the rhs. @@ -7591,10 +7644,10 @@ solenoidal_options = empty } )?, - ## Solver options for the linear solve. + ## Solver options for the (outer) linear solve. ## This method requires the inversion of a projection matrix. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } } diff --git a/schemas/fluidity_options.rng b/schemas/fluidity_options.rng index fe9317a5ee..bfb988280b 100644 --- a/schemas/fluidity_options.rng +++ b/schemas/fluidity_options.rng @@ -3933,7 +3933,7 @@ formulations Solver options are necessary if you're not lumping your mass or if you're field isn't dg - + @@ -3956,7 +3956,7 @@ formulations Solver - + @@ -3985,7 +3985,7 @@ formulations Solver - + @@ -5458,7 +5458,7 @@ requires a distinct name for the options dictionary. Solver - only necessary in combination with an explicit no_normal_stress free_surface bc under Velocity. - + @@ -6890,7 +6890,7 @@ strong Dirichlet boundary condition of 0 on all boundaries. Solver - + @@ -7030,7 +7030,7 @@ less accurate but faster and might give smoother result. - + @@ -7200,7 +7200,7 @@ or with a CG but with the continuity tested with the CV dual. - + @@ -7221,7 +7221,7 @@ where _c and _i denote the compressible and incompressible phases respectively.< - + @@ -7463,7 +7463,7 @@ and the Foam Velocity is obtained by taking its gradient. - + @@ -7771,7 +7771,7 @@ less accurate but faster and might give smoother result. - + @@ -9337,20 +9337,74 @@ vector fields, such as incompressible velocity! Options for the mass matrix of the field being interpolated - - Lump the mass matrix for the assembly of the projection matrix -(not for the initial galerkin projection) + + + Solver options for the (inner) solve. + + Specify the inner matrix (IM) to form the projection schur complement (C^T*IM^{-1}*C). +Use the full mass matrix (the only option). + + FullMassMatrix + + + Solver options to solve the inner (mass) matrix in the schur complement solve, also used +to solve the correction solve afterwards. + + + + + + Specify the preconditioner matrix to use on the schur complement. -Required when using interpolating continuous fields - - - Lump on the submesh. +Assemble a Schur Complement using a lumped mass inner matrix. + + LumpedSchurComplement + + + + + + + + + Specify the preconditioner matrix to use on the schur complement. + +DiagonalSchurComplement = C_P^T * [(M)_diagonal]^-1 * C + + DiagonalSchurComplement + + + + + Specify the preconditioner matrix to use on the schur complement. + +Lagrange Mass Matrix. + + LagrangeMassMatrix + + + + + Specify the preconditioner matrix to use on the schur complement. + + NoPreconditionerMatrix + + + + + + + Lump the mass matrix for the assembly of the projection matrix + + + Lump on the submesh. This only works for simplex meshes and is only strictly valid on 2d meshes. - - - - + + + + + @@ -9403,9 +9457,20 @@ considered when selecting the lagrange multiplier solver options. - - - + + + Reference node (node at which lagrange multiplier = 0). Note that the node number must be less than the total number of nodes. +If running in parallel, the node number must be less than the total number of nodes of the first processor. +** Note - it is also an option to remove the null-space of the residual vector. This +option is available under solvers + + + + Input coordinates of desired reference node. If a node does not exist at these +coordinates, the nearest vertex will be selected. + + + @@ -9450,9 +9515,9 @@ vector fields, such as pressure to incompressible velocity! - Solver options for the linear solve. + Solver options for the (outer) linear solve. This method requires the inversion of a projection matrix. - + diff --git a/schemas/mesh_options.rnc b/schemas/mesh_options.rnc index e162a0c508..08173f6b27 100644 --- a/schemas/mesh_options.rnc +++ b/schemas/mesh_options.rnc @@ -64,7 +64,7 @@ mesh_info = element solver { ## Solver is used to apply projections to div-conforming ## space, and in elliptic solvers - linear_solver_options_sym + linear_solver_options_sym_scalar } }?, ## Make mesh periodic diff --git a/schemas/mesh_options.rng b/schemas/mesh_options.rng index 7dfc2033db..b4e531ac86 100644 --- a/schemas/mesh_options.rng +++ b/schemas/mesh_options.rng @@ -94,7 +94,7 @@ solver. - + Solver is used to apply projections to div-conforming space, and in elliptic solvers diff --git a/schemas/prescribed_field_options.rnc b/schemas/prescribed_field_options.rnc index a09560f635..9d82a7eeb9 100644 --- a/schemas/prescribed_field_options.rnc +++ b/schemas/prescribed_field_options.rnc @@ -1,4 +1,3 @@ -# Default child of prescribed scalar field # This is a choice of ways of inputing the prescribed field prescribed_scalar_field = ( @@ -170,4 +169,4 @@ prescribed_values_dim_minus_one_tensor_field = input_choice_dim_minus_one_tensor_field } )+ - ) \ No newline at end of file + ) diff --git a/schemas/prescribed_field_options.rng b/schemas/prescribed_field_options.rng index 1b4f2c82c7..9920dbc715 100644 --- a/schemas/prescribed_field_options.rng +++ b/schemas/prescribed_field_options.rng @@ -1,9 +1,6 @@ - + diff --git a/schemas/prognostic_field_options.rnc b/schemas/prognostic_field_options.rnc index d5a2a83c94..dc8f7e1ed7 100644 --- a/schemas/prognostic_field_options.rnc +++ b/schemas/prognostic_field_options.rnc @@ -644,7 +644,7 @@ prognostic_velocity_field = element length_scale_type { 'scalar'|'tensor' }, ## Solver options are required for calculation of Helmholtz-smoothed velocity. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, ## Spatially varying Smagorinsky coefficient. element scalar_field { @@ -1332,7 +1332,7 @@ prognostic_velocity_field = }?, ## Solver options are necessary if you're not lumping your mass or if you're field isn't dg element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }?, ## Choose whether the surface tension term in the momentum equation is integrated by parts or not element integrate_by_parts { @@ -1426,7 +1426,7 @@ prognostic_pressure_field = ) }, ( - ## Reference node (Node at which pressure = 0.) Note that the node number must be less than the total number of nodes. + ## Reference node (node at which pressure = 0). Note that the node number must be less than the total number of nodes. ## If running in parallel, the node number must be less than the total number of nodes of the first processor. ## ** Note - it is also an option to remove the null-space of the residual vector. This ## option is available under solvers @@ -1488,7 +1488,7 @@ prognostic_pressure_field = ## Solver options to solve the inner (mass) matrix in the schur complement solve, also used ## to solve the correction solve afterwards. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } }| ## Specify the inner matrix (IM) to form the projection schur complement (C^T*IM^{-1}*C). @@ -1547,7 +1547,7 @@ prognostic_pressure_field = }, ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, ( ## Initial condition for WholeMesh @@ -1648,7 +1648,7 @@ prognostic_geostrophic_pressure_field = ( ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } ), ( @@ -1704,7 +1704,7 @@ prognostic_vertical_balance_pressure_field = ( ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, ## Apply a strong dirichlet boundary condition to VerticalBalancePressure. ## This is normally be a homogeneous bc on the top surface. @@ -1815,7 +1815,7 @@ prognostic_foam_velocity_potential_field = ( ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } ), ( diff --git a/schemas/prognostic_field_options.rng b/schemas/prognostic_field_options.rng index 4c4d2fe5eb..16c823245e 100644 --- a/schemas/prognostic_field_options.rng +++ b/schemas/prognostic_field_options.rng @@ -780,7 +780,7 @@ The 'tensor' length scale uses the metric from the adaptivity process. Using thi Solver options are required for calculation of Helmholtz-smoothed velocity. - + @@ -1644,7 +1644,7 @@ Note that the default value is zero. Solver options are necessary if you're not lumping your mass or if you're field isn't dg - + @@ -1765,7 +1765,7 @@ the continuity equation is never integrated by parts. - Reference node (Node at which pressure = 0.) Note that the node number must be less than the total number of nodes. + Reference node (node at which pressure = 0). Note that the node number must be less than the total number of nodes. If running in parallel, the node number must be less than the total number of nodes of the first processor. ** Note - it is also an option to remove the null-space of the residual vector. This option is available under solvers @@ -1839,7 +1839,7 @@ Make sure you've not lumped your mass in the velocity spatial_discretisation if Solver options to solve the inner (mass) matrix in the schur complement solve, also used to solve the correction solve afterwards. - + @@ -1915,7 +1915,7 @@ viscosity tensors. Solver - + @@ -2047,7 +2047,7 @@ be used with the solver/remove_null_space option. Solver - + @@ -2118,7 +2118,7 @@ this only makes sense when excluding coriolis. Solver - + @@ -2241,7 +2241,7 @@ Only required for continuous_galerkin spatial_discretisations! Solver - + diff --git a/schemas/solvers.rnc b/schemas/solvers.rnc index 9fce5ef6b1..830d3f50b7 100644 --- a/schemas/solvers.rnc +++ b/schemas/solvers.rnc @@ -36,6 +36,37 @@ linear_solver_options_asym_scalar = generic_solver_options_scalar ) +linear_solver_options_sym_scalar = + ( + ## Iterative (Krylov) method to solve the linear discretised equation + ## Given are the most frequently used methods. The solution is done + ## by the PETSc library. Many more methods are provided. + ## + ( + kspcg_options| + kspgmres_options| + ksppreonly_options| + ksprichardson_options| + kspother_options + ), + ## Preconditioner to be used in combination with the iterative method. + ( + pcsor_options| + pceisenstat_options| + pcilu_options| + pcicc_options| + pclu_options| + pcmg_options| + pcprometheus_options| + pchypre_options| + pcbjacobi_options| + pcasm_options| + pcksp_options| + pcother_options + ), + generic_solver_options_scalar + ) + linear_solver_options_asym_vector = ( ## Iterative (Krylov) method to solve the linear discretised equation @@ -67,7 +98,7 @@ linear_solver_options_asym_vector = generic_solver_options_vector ) -linear_solver_options_sym = +linear_solver_options_sym_vector = ( ## Iterative (Krylov) method to solve the linear discretised equation ## Given are the most frequently used methods. The solution is done @@ -78,24 +109,24 @@ linear_solver_options_sym = kspgmres_options| ksppreonly_options| ksprichardson_options| - kspother_options + kspother_options ), ## Preconditioner to be used in combination with the iterative method. ( pcsor_options| pceisenstat_options| pcilu_options| - pcicc_options| pclu_options| pcmg_options| pcprometheus_options| pchypre_options| pcbjacobi_options| pcasm_options| - pcksp_options| + pcgamg_options| + pcfieldsplit_options| pcother_options ), - generic_solver_options_scalar + generic_solver_options_vector ) # #################################################################### @@ -380,7 +411,7 @@ pcksp_options = } ) -# this is a copy linear_solver_options_sym, but with preconditioner "ksp" +# this is a copy linear_solver_options_sym_scalar, but with preconditioner "ksp" # removed to avoid infinite recursion pc_ksp_solver_options = ( @@ -679,7 +710,7 @@ galerkin_projection_mass_options = ## Solver options. Required for a continuous consistent mass ## projection. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } ) @@ -698,6 +729,6 @@ galerkin_projection_mass_options_submesh = ## Solver options. Required for a continuous consistent mass ## projection. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } ) diff --git a/schemas/solvers.rng b/schemas/solvers.rng index 7ba862b981..1a4fdf1ab2 100644 --- a/schemas/solvers.rng +++ b/schemas/solvers.rng @@ -37,6 +37,35 @@ by the PETSc library. Many more methods are provided. + + + Iterative (Krylov) method to solve the linear discretised equation +Given are the most frequently used methods. The solution is done +by the PETSc library. Many more methods are provided. + + + + + + + + + Preconditioner to be used in combination with the iterative method. + + + + + + + + + + + + + + + Iterative (Krylov) method to solve the linear discretised equation @@ -66,7 +95,7 @@ by the PETSc library. Many more methods are provided. - + Iterative (Krylov) method to solve the linear discretised equation Given are the most frequently used methods. The solution is done @@ -83,17 +112,17 @@ by the PETSc library. Many more methods are provided. - - + + - + @@ -746,7 +775,7 @@ run for multiple solves. Solver options. Required for a continuous consistent mass projection. - + @@ -767,7 +796,7 @@ strictly valid on 2d meshes. Solver options. Required for a continuous consistent mass projection. - + diff --git a/schemas/test_advection_diffusion_options.rnc b/schemas/test_advection_diffusion_options.rnc index 7c333c4690..c34ef05d62 100644 --- a/schemas/test_advection_diffusion_options.rnc +++ b/schemas/test_advection_diffusion_options.rnc @@ -1086,7 +1086,7 @@ prognostic_velocity_field = }?, ## Solver options are necessary if you're not lumping your mass or if you're field isn't dg element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }?, ## Choose whether the surface tension term in the momentum equation is integrated by parts or not element integrate_by_parts { @@ -1110,7 +1110,7 @@ prognostic_scalar_field = ( ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, ( ## Initial condition for WholeMesh @@ -1425,7 +1425,7 @@ prognostic_pressure_field = element inner_matrix { attribute name { "FullMassMatrix" }, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } }| ## Specify the inner matrix (IM) to form the projection schur complement (C^T*IM^{-1}*C). @@ -1473,7 +1473,7 @@ prognostic_pressure_field = }, ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, ( ## Initial condition for WholeMesh @@ -1549,7 +1549,7 @@ prognostic_geostrophic_pressure_field = ( ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } ), ( @@ -1607,7 +1607,7 @@ prognostic_balance_pressure_field = # Solver block is the same as prognostic_scalar_field ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, # Alas, no initial_condition either, so we'd better not checkpointing it... ## Disables checkpointing of this field @@ -1743,7 +1743,7 @@ diagnostic_cv_gradient_vector_field = }?, ## Solver options are necessary if you're not lumping your mass or if you're field isn't dg element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }?, ## Normalise the gradient by its magnitude element normalise { @@ -1762,7 +1762,7 @@ diagnostic_gradient_vector_field = ( ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, diagnostic_output_options, diagnostic_vector_stat_options, @@ -1789,7 +1789,7 @@ diagnostic_fe_divergence_scalar_field = ( ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, diagnostic_output_options, diagnostic_scalar_stat_options, @@ -2476,7 +2476,7 @@ adaptivity_preprocessing = real_dim_symmetric_tensor }, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } ) @@ -3939,7 +3939,7 @@ scalar_field_choice = velocity_mesh_choice, ## Solver element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }, diagnostic_scalar_field }| @@ -4046,7 +4046,7 @@ scalar_field_choice = empty }?, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }?, diagnostic_scalar_field }| @@ -4866,7 +4866,7 @@ vector_field_choice = empty }?, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }?, diagnostic_vector_field }| @@ -4890,7 +4890,7 @@ vector_field_choice = internal_algorithm, velocity_mesh_choice, element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }?, diagnostic_vector_field }| @@ -5312,7 +5312,7 @@ linear_solver_options_asym = generic_solver_options ) -linear_solver_options_sym = +linear_solver_options_sym_scalar = ( ## Iterative (Krylov) method to solve the linear discretised equation ## Given are the most frequently used methods. The solution is done @@ -5461,7 +5461,7 @@ linear_solver_options_sym = generic_solver_options ) -# this is a copy linear_solver_options_sym, but with preconditioner "ksp" +# this is a copy linear_solver_options_sym_scalar, but with preconditioner "ksp" # removed to avoid infinite recursion pc_ksp_solver_options = ( @@ -6603,7 +6603,7 @@ scalar_equation_choice = attribute name { string } } ) - }| + } ) ) @@ -7060,7 +7060,7 @@ continuous_discontinuous_projection = ## conservation properties are affected by the tolerance of the ## linear solve. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar }| ## Lump the mass matrix on the left hand side of the galerkin projection. ## Hence solver options aren't necessary. @@ -7171,7 +7171,7 @@ solenoidal_options = ## Solver options for the linear solve. ## This method requires the inversion of a projection matrix. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } } @@ -7190,7 +7190,7 @@ balanced_interpolation = element balanced_interpolation { ## The solver options for computing the balanced velocity to subtract. element solver { - linear_solver_options_sym + linear_solver_options_sym_scalar } } diff --git a/schemas/test_advection_diffusion_options.rng b/schemas/test_advection_diffusion_options.rng index 1e0c90cf3c..a558d76578 100644 --- a/schemas/test_advection_diffusion_options.rng +++ b/schemas/test_advection_diffusion_options.rng @@ -1260,35 +1260,6 @@ only if the mass matrix is lumped (lump_mass_matrix). - - - Elastic parameters for elastic and visco-elastic materials -For a linearly elastic solid -In legacy elastic solids were run with SOLIDS = 2 from solidity_options.inp. -In gem ONEMU and CONMU = .TRUE. -and RMUPZZ taken as the isotropic Young`s modulus -UNDER DEVELOPMENT - - currently only works for lagrangian meshes - - only single materials tested so far - - momentum equations assembled in solid3d so not all - discretisation options above are valid - - Elasticity - - - 2 - - - - - - - - - - - - SurfaceTension @@ -1312,7 +1283,7 @@ UNDER DEVELOPMENT Solver options are necessary if you're not lumping your mass or if you're field isn't dg - + @@ -1325,56 +1296,6 @@ UNDER DEVELOPMENT - - - Cohesion for plastic materials -UNDER DEVELOPMENT - - currently only works for lagrangian meshes - - only single materials tested so far - - momentum equations assembled in solid3d so not all - discretisation options above are valid - - Cohesion - - - 0 - - - - - - - - - - - - - - - Friction Angle for plastic materials -UNDER DEVELOPMENT - - currently only works for lagrangian meshes - - only single materials tested so far - - momentum equations assembled in solid3d so not all - discretisation options above are valid - - FrictionAngle - - - 0 - - - - - - - - - - - - @@ -1389,7 +1310,7 @@ UNDER DEVELOPMENT Solver - + @@ -1744,7 +1665,7 @@ Make sure you've not lumped your mass in the velocity spatial_discretisation if FullMassMatrix - + @@ -1805,7 +1726,7 @@ or if density varies a lot or if not using a Boussinesque approximation) Solver - + @@ -1899,7 +1820,7 @@ want this. Solver - + @@ -1973,7 +1894,7 @@ requires a distinct name for the options dictionary. --> Solver - + @@ -2135,7 +2056,7 @@ formulations Solver options are necessary if you're not lumping your mass or if you're field isn't dg - + @@ -2157,7 +2078,7 @@ formulations Solver - + @@ -2184,7 +2105,7 @@ formulations Solver - + @@ -2946,7 +2867,7 @@ fields. - + @@ -3947,8 +3868,7 @@ Requires a continuous mesh. - PhaseVolumeFraction: -Required in porous_media problem type + PhaseVolumeFraction: 0 @@ -3966,22 +3886,6 @@ Required in porous_media problem type - - Electrical Potential: -Required in electrokinetic, electrothermal -and electrochemical problems -(sub-option of porous_media problem type) - - 0 - - - ElectricalPotential - - - - - - @@ -7549,111 +7244,6 @@ velocity. - - - - - Obtain values from point and radius file. - -First line of file is free to use (for comments) -Second line must contain the number of particles -Third and fourth line are again for comments. -Following lines include 10 columns, corresponding to -the particle's x, y, and z positions, followed by the radius, then -velocities in x, y, and z directions, followed by angular velocities -in the x, y, and z directions. - - from_input_file - - - - - - - Two python scripts must be provided. The script is cycled over each particle. -One script for particle position (output is tuple of position coords) -Second script is for particle radius (output is tuple of position coords) -Third script is for particle translational velocity. -Fourth script is for particle angular velocity. (Note: particles -have a no slip boundary condition at the surface, so this angular velocity -WILL matter to the flow.) -Python functions should be of the form: - def val(X, t): - Function code - return # Return value -where X is a tuple of length geometry dimension. - X[0] contains the number of the particle (in real format) - - python_script - - - - - use_simple_dynamics - - - - - - - - - - - - - - - - - - - - - - - - - Using y3D to model dynamics. Filename of input file for y3D must -be specified. - - use_y3D - - - - - - - Using femdem 2D to model dynamics. Filename of input file must -be specified. - - use_2Dfemdem - - - - - - - Using femdem 3D to model dynamics. Filename of input file must -be specified. - - use_3Dfemdem - - - - - - - - - - - - Cap the min and max values of this field when using @@ -8266,13 +7856,6 @@ Whatever field is selected must be present. - - Option to solve for electrical potential from -electrokinetic, electrochemical or electrothermal sources - - ElectricalPotential - - @@ -8787,7 +8370,7 @@ deviations to nodes that have less than their bounds. This method requires the inversion of a mass matrix. Note that conservation properties are affected by the tolerance of the linear solve. - + Lump the mass matrix on the left hand side of the galerkin projection. @@ -8914,7 +8497,7 @@ vector fields, such as pressure to incompressible velocity! Solver options for the linear solve. This method requires the inversion of a projection matrix. - + @@ -8932,7 +8515,7 @@ interpolation, it will be balanced afterwards also. Note well: this only makes sense for velocity. The solver options for computing the balanced velocity to subtract. - + From 83da6e80df1b514e246433bf2c647d08e5fbd580 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Wed, 7 Jul 2021 01:03:23 -0400 Subject: [PATCH 02/12] Try turning on github actions. --- .github/workflows/ubuntu.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ubuntu.yml b/.github/workflows/ubuntu.yml index f77d783cd8..26fcdf165f 100644 --- a/.github/workflows/ubuntu.yml +++ b/.github/workflows/ubuntu.yml @@ -4,6 +4,7 @@ on: push: branches: - main + - solenoidal_full_schur pull_request: branches: - main From e193effea9b55b8c34dff5ac70393ab572d8b044 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Wed, 7 Jul 2021 01:24:33 -0400 Subject: [PATCH 03/12] TO REVERT: also turn on github actions for the prefixed branch. --- .github/workflows/ubuntu.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ubuntu.yml b/.github/workflows/ubuntu.yml index 26fcdf165f..7b63de6495 100644 --- a/.github/workflows/ubuntu.yml +++ b/.github/workflows/ubuntu.yml @@ -5,6 +5,7 @@ on: branches: - main - solenoidal_full_schur + - cianwilson/solenoidal_full_schur pull_request: branches: - main From a99269686d4f72cbae0b1d5f741d678de21f824a Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Wed, 7 Jul 2021 01:46:19 -0400 Subject: [PATCH 04/12] Needed to run make makefiles one more time apparently. --- assemble/Makefile.dependencies | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/assemble/Makefile.dependencies b/assemble/Makefile.dependencies index 5a12ff33f3..00a1be0e87 100644 --- a/assemble/Makefile.dependencies +++ b/assemble/Makefile.dependencies @@ -350,7 +350,7 @@ Full_Projection.o ../include/full_projection.mod: Full_Projection.F90 \ ../include/boundary_conditions.mod \ ../include/boundary_conditions_from_options.mod \ ../include/data_structures.mod ../include/elements.mod ../include/fdebug.h \ - ../include/fields.mod ../include/fldebug.mod \ + ../include/field_options.mod ../include/fields.mod ../include/fldebug.mod \ ../include/global_parameters.mod ../include/halos.mod \ ../include/multigrid.mod ../include/parallel_tools.mod \ ../include/petsc_legacy.h ../include/petsc_solve_state_module.mod \ From 4b6d391d23ada2b4ea3cead925fa872e1d5b2710 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Wed, 7 Jul 2021 01:59:56 -0400 Subject: [PATCH 05/12] Beginnings of a test (doesn't work yet due to pschur problems) --- tests/mmat-interpolation-p2p1cv/Makefile | 7 + .../change_options.py | 34 + .../mmat-interpolation-p2p1cv.xml | 350 ++++++++++ .../mmat-interpolation.flml | 614 ++++++++++++++++++ .../src/2d_square.geo | 15 + 5 files changed, 1020 insertions(+) create mode 100644 tests/mmat-interpolation-p2p1cv/Makefile create mode 100755 tests/mmat-interpolation-p2p1cv/change_options.py create mode 100644 tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml create mode 100644 tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml create mode 100644 tests/mmat-interpolation-p2p1cv/src/2d_square.geo diff --git a/tests/mmat-interpolation-p2p1cv/Makefile b/tests/mmat-interpolation-p2p1cv/Makefile new file mode 100644 index 0000000000..b2e7c898e8 --- /dev/null +++ b/tests/mmat-interpolation-p2p1cv/Makefile @@ -0,0 +1,7 @@ +input: clean + gmsh -2 src/2d_square.geo + cp src/2d_square.msh . + +clean: + rm -rf *.err *.log *.stat *.vtu *.node *.ele *.edge *checkpoint* *convergence* *.log-0 *.err-0 src/2d_square.msh *.msh \ + matrixdump matrixdump.info diff --git a/tests/mmat-interpolation-p2p1cv/change_options.py b/tests/mmat-interpolation-p2p1cv/change_options.py new file mode 100755 index 0000000000..0224050f6b --- /dev/null +++ b/tests/mmat-interpolation-p2p1cv/change_options.py @@ -0,0 +1,34 @@ +#!/usr/bin/env python3 + +import sys + +file_name = ["mmat-interpolation_2_checkpoint.flml"] + +for f in range(len(file_name)): + try: + flml_file = open(file_name[f], "r") + flml_options = flml_file.readlines() + flml_file.close() + except: + sys.stderr.write("Error: change_options failed to read options file\n") + sys.exit(1) + + for i in range(len(flml_options)): + if "" in flml_options[i]: + while not "" in flml_options[i] and i < len(flml_options): + flml_options[i] = "" + i += 1 + if i < len(flml_options): + flml_options[i] = "" + break + + try: + flml_file = open(file_name[f], "w") + flml_file.writelines(flml_options) + flml_file.close() + except: + sys.stderr.write("Error: change_options failed to write options file\n") + sys.exit(1) + + + diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml new file mode 100644 index 0000000000..23be75ccb3 --- /dev/null +++ b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml @@ -0,0 +1,350 @@ + + + Multimaterial adaptivity and interpolation test + + flml 2dadapt + + fluidity -v2 -l mmat-interpolation.flml && mv fluidity.log-0 first.log && mv fluidity.err-0 first.err && ./change_options.py && fluidity -v2 -l mmat-interpolation_2_checkpoint.flml && mv fluidity.log-0 checkpoint.log && mv fluidity.err-0 checkpoint.err + + + +import os +files = os.listdir("./") +solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files + + +from fluidity_tools import stat_parser as stat +material1integralstart = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +material1integralend = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +material1maxstart = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] + + +from fluidity_tools import stat_parser as stat +material1maxend = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +material1minstart = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] + + +from fluidity_tools import stat_parser as stat +material1minend = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +material2integralstart = stat("mmat-interpolation.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +material2integralend = stat("mmat-interpolation.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +material3integralstart = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +material3integralend = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +material3maxstart = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] + + +from fluidity_tools import stat_parser as stat +material3maxend = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +material3minstart = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] + + +from fluidity_tools import stat_parser as stat +material3minend = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +divergenceminstart = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] + + +from fluidity_tools import stat_parser as stat +divergenceminend = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +divergencemaxstart = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] + + +from fluidity_tools import stat_parser as stat +divergencemaxend = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +cflmaxstart = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] + + +from fluidity_tools import stat_parser as stat +cflmaxend = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +lambdamaxstart = stat("mmat-interpolation.stat")["Material1"]["Lambda"]["max"][0] + + +from fluidity_tools import stat_parser as stat +lambdamaxend = stat("mmat-interpolation.stat")["Material1"]["Lambda"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +lambdaminstart = stat("mmat-interpolation.stat")["Material1"]["Lambda"]["min"][0] + + +from fluidity_tools import stat_parser as stat +lambdaminend = stat("mmat-interpolation.stat")["Material1"]["Lambda"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +material1integralstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +material1integralend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +material1maxstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] + + +from fluidity_tools import stat_parser as stat +material1maxend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +material1minstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] + + +from fluidity_tools import stat_parser as stat +material1minend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +material2integralstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +material2integralend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +material3integralstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +material3integralend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +material3maxstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] + + +from fluidity_tools import stat_parser as stat +material3maxend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +material3minstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] + + +from fluidity_tools import stat_parser as stat +material3minend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +divergenceminstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] + + +from fluidity_tools import stat_parser as stat +divergenceminend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +divergencemaxstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] + + +from fluidity_tools import stat_parser as stat +divergencemaxend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +cflmaxstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] + + +from fluidity_tools import stat_parser as stat +cflmaxend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +lambdamaxstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["Lambda"]["max"][0] + + +from fluidity_tools import stat_parser as stat +lambdamaxend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["Lambda"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +lambdaminstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["Lambda"]["min"][0] + + +from fluidity_tools import stat_parser as stat +lambdaminend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["Lambda"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +absolutedifference = stat("mmat-interpolation.stat")["Material3"]["ScalarAbsoluteDifference"]["max"][1] + + + + + assert(solvers_converged) + + +print('mass loss = ', abs(material1integralstart-material1integralend)) +assert abs(material1integralstart-material1integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + + +print('mass loss = ', abs(material2integralstart-material2integralend)) +assert abs(material2integralstart-material2integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + + +print('mass loss = ', abs(material3integralstart-material3integralend)) +assert abs(material3integralstart-material3integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + + +print('divergence tolerance = ', max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend))) +assert max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend)) < 1.E-8 + + +assert abs(material1maxstart-material1maxend) < 1.E-10 + + +assert abs(material3maxstart-material3maxend) < 1.E-10 + + +assert abs(material1minstart-material1minend) < 1.E-10 + + +assert abs(material3minstart-material3minend) < 1.E-10 + + +assert abs(divergencemaxstart) < 1.E-8 + + +assert abs(divergenceminstart) < 1.E-8 + + +assert abs(divergencemaxend) < 1.E-8 + + +assert abs(divergenceminend) < 1.E-8 + + +assert abs(cflmaxstart-0.1) < 1.E-10 + + +assert abs(cflmaxend-0.1) < 1.E-10 + + +print('mass loss = ', abs(material1integralstart_checkpoint-material1integralend_checkpoint)) +assert abs(material1integralstart_checkpoint-material1integralend_checkpoint) < max(abs(divergencemaxstart_checkpoint), abs(divergencemaxend_checkpoint), abs(divergenceminstart_checkpoint), abs(divergenceminend_checkpoint), 1.E-8) + + +print('mass loss = ', abs(material2integralstart_checkpoint-material2integralend_checkpoint)) +assert abs(material2integralstart_checkpoint-material2integralend_checkpoint) < max(abs(divergencemaxstart_checkpoint), abs(divergencemaxend_checkpoint), abs(divergenceminstart_checkpoint), abs(divergenceminend_checkpoint), 1.E-8) + + +print('mass loss = ', abs(material3integralstart_checkpoint-material3integralend_checkpoint)) +assert abs(material3integralstart_checkpoint-material3integralend_checkpoint) < max(abs(divergencemaxstart_checkpoint), abs(divergencemaxend_checkpoint), abs(divergenceminstart_checkpoint), abs(divergenceminend_checkpoint), 1.E-8) + + +print('divergence tolerance = ', max(abs(divergencemaxstart_checkpoint), abs(divergencemaxend_checkpoint), abs(divergenceminstart_checkpoint), abs(divergenceminend_checkpoint))) +assert max(abs(divergencemaxstart_checkpoint), abs(divergencemaxend_checkpoint), abs(divergenceminstart_checkpoint), abs(divergenceminend_checkpoint)) < 1.E-8 + + +assert abs(material1maxstart_checkpoint-material1maxend_checkpoint) < 1.E-10 + + +assert abs(material3maxstart_checkpoint-material3maxend_checkpoint) < 1.E-10 + + +assert abs(material1minstart_checkpoint-material1minend_checkpoint) < 1.E-10 + + +assert abs(material3minstart_checkpoint-material3minend_checkpoint) < 1.E-10 + + +assert abs(divergencemaxstart_checkpoint) < 1.E-8 + + +assert abs(divergenceminstart_checkpoint) < 1.E-8 + + +assert abs(divergencemaxend_checkpoint) < 1.E-8 + + +assert abs(divergenceminend_checkpoint) < 1.E-8 + + +assert abs(cflmaxstart_checkpoint-0.1) < 1.E-10 + + +assert abs(cflmaxend_checkpoint-0.1) < 1.E-10 + + +print('mass loss = ', abs(material1integralstart_checkpoint-material1integralend)) +assert abs(material1integralstart_checkpoint-material1integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + + +print('mass loss = ', abs(material2integralstart_checkpoint-material2integralend)) +assert abs(material2integralstart_checkpoint-material2integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + + +print('mass loss = ', abs(material3integralstart_checkpoint-material3integralend)) +assert abs(material3integralstart_checkpoint-material3integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + + + + +assert abs(lambdamaxstart) > 1.E-7 + + +assert abs(lambdaminstart) > 1.E-7 + + +assert abs(lambdamaxend) > 1.E-7 + + +assert abs(lambdaminend) > 1.E-7 + + +assert abs(lambdamaxstart_checkpoint) < 1.E-7 + + +assert abs(lambdaminstart_checkpoint) < 1.E-7 + + +assert abs(lambdamaxend_checkpoint) > 1.E-7 + + +assert abs(lambdaminend_checkpoint) > 1.E-7 + + +assert abs(absolutedifference) > 1.E-7 + + + diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml b/tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml new file mode 100644 index 0000000000..e4e09da101 --- /dev/null +++ b/tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml @@ -0,0 +1,614 @@ + + + + mmat-interpolation + + + multimaterial + + + + 2 + + + + + + + + + + + + + + + 2 + + + + + + + + + + + + + + + + + + 5 + + + 1 + + + + + + vtk + + + + 10 + 0.36440821846106625 + + + + + + 1 + + + + + + + + 0.0 + + + 0.00427561316453 + co0p1-doubleres = 0.000893158065237 +co0p1-halfres = 0.00427561316453 + + + 14.576328738442649 + 14.576328738442649 +will get cut off by final_timestep + + + 25 + + + + 0.1 + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + def val(X,t): + from math import sin, cos + # Shear rotation about origin. + return (sin(X[0])*cos(X[1]), -1.0*cos(X[0])*sin(X[1])) + + + + + + + + + + + + + + + + + + + + + + + + + + 1.e-6 + + + 100 + + + + + + + + + + + + + + + + + + + + + + + + 1.E-10 + + + 10000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + 0.0 + + + + 0.0 + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.2*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + 8 9 10 11 + + + + + + + + + + + + + + + + + + + + + + + + 0.11 + + + + + + + + + + + + + + + + + + + 10000 + + + + + + + 1.E-10 + + + 10000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + + + + + + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + 0.0 + + + + 0.0 + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.55*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + 8 9 10 11 + + + + + + + + + + + + + + + + + + + + + + + + 0.11 + + + + + + + + + + + + + + + + + + + 10000 + + + + + + + 1.E-10 + + + 10000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.55*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + + + + + + + + + + + + + 10 + + + 5000 + + + + + + 1.0 0.0 0.0 1.0 + + + + + + + + 0.01 0.0 0.0 0.01 + + + + + + + 1.5 0.0 0.0 1.5 + + + + + + diff --git a/tests/mmat-interpolation-p2p1cv/src/2d_square.geo b/tests/mmat-interpolation-p2p1cv/src/2d_square.geo new file mode 100644 index 0000000000..70cb64b6af --- /dev/null +++ b/tests/mmat-interpolation-p2p1cv/src/2d_square.geo @@ -0,0 +1,15 @@ +Point(1) = {0.0,0.0,0,0.062831853071795868}; +Point(2) = {3.1415926535897931,0.0,0,0.062831853071795868}; +Point(3) = {3.1415926535897931,3.1415926535897931,0,0.062831853071795868}; +Point(4) = {0.0,3.1415926535897931,0,0.062831853071795868}; +Line(1) = {4,3}; +Line(2) = {3,2}; +Line(3) = {2,1}; +Line(4) = {1,4}; +Line Loop(5) = {1,2,3,4}; +Plane Surface(6) = {5}; +Physical Surface(7) = {6}; +Physical Line(8) = {3}; +Physical Line(9) = {2}; +Physical Line(10) = {1}; +Physical Line(11) = {4}; From 778b98f54161b0d8a7be3d1015593af2ca74c9e0 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Wed, 7 Jul 2021 19:03:04 -0400 Subject: [PATCH 06/12] Some fixes to Solenoidal interpolation where the dg logic wasn't quite right. Also some fixes to two tests that failed the updated logic because they were themselves inconsistent. --- assemble/Solenoidal_interpolation.F90 | 8 ++++---- .../hyperc-shear-input.flml | 4 ++-- .../hyperc-shear-input.flml | 8 ++++++++ 3 files changed, 14 insertions(+), 6 deletions(-) diff --git a/assemble/Solenoidal_interpolation.F90 b/assemble/Solenoidal_interpolation.F90 index 0c11003535..5d33b57e06 100644 --- a/assemble/Solenoidal_interpolation.F90 +++ b/assemble/Solenoidal_interpolation.F90 @@ -152,8 +152,8 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & discontinuous = have_option(trim(l_option_path)//& "/interpolated_field/discontinuous") - if (dg) then - assert(discontinuous) + if (discontinuous.and..not.dg) then + FLExit("Using discontinuous solenoidal projection options but field is not discontinuous.") end if assemble_lumped_mass = have_option(trim(l_option_path)//& @@ -313,7 +313,7 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & if(assemble_lumped_mass.and.lump_mass_on_submesh) then ! get the mass matrix lumped on the submesh call compute_lumped_mass_on_submesh(local_state, field_lumped_mass) - else if (discontinuous) then + else if (discontinuous.and..not.assemble_lumped_mass) then ! get the inverse discontinuous mass matrix call construct_inverse_mass_matrix_dg(inverse_field_mass, v_field, coordinate) end if @@ -350,7 +350,7 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & call apply_dirichlet_conditions_inverse_mass(inverse_field_lumped_mass_vector, v_field) call assemble_masslumped_cmc(cmc_m, ctp_m, inverse_field_lumped_mass_vector, ct_m) - else if(dg) then + else if(discontinuous) then call assemble_cmc_dg(cmc_m, ctp_m, ct_m, inverse_field_mass) else FLExit("Not possible to not lump the mass if not dg or full_schur.") diff --git a/tests/hyperc-shear-moving-mesh-divfree-pseudolagrangian/hyperc-shear-input.flml b/tests/hyperc-shear-moving-mesh-divfree-pseudolagrangian/hyperc-shear-input.flml index 615405ff65..33ec76d810 100644 --- a/tests/hyperc-shear-moving-mesh-divfree-pseudolagrangian/hyperc-shear-input.flml +++ b/tests/hyperc-shear-moving-mesh-divfree-pseudolagrangian/hyperc-shear-input.flml @@ -110,9 +110,9 @@ - + - + diff --git a/tests/hyperc-shear-moving-mesh-divfree/hyperc-shear-input.flml b/tests/hyperc-shear-moving-mesh-divfree/hyperc-shear-input.flml index 45936b5e2d..5afc7bc3d9 100644 --- a/tests/hyperc-shear-moving-mesh-divfree/hyperc-shear-input.flml +++ b/tests/hyperc-shear-moving-mesh-divfree/hyperc-shear-input.flml @@ -21,6 +21,14 @@ + + + 0 + + + + discontinuous + From 67b43fc0a600551c2a75d13f39899527ffe538d8 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Wed, 7 Jul 2021 19:31:22 -0400 Subject: [PATCH 07/12] Try using the row halo from ct_p instead of the (not necessarily in existence) preconditioner_matrix. Also allow correct velocity to take in a option path as it also has assumed locations for velocity solver options. Add a test that appears to pass that uses Schur complement solenoidal projection. Need to look at this more. --- assemble/Full_Projection.F90 | 18 +++++---------- assemble/Momentum_CG.F90 | 21 +++++++++++------- assemble/Solenoidal_interpolation.F90 | 8 ++++--- .../mmat-interpolation.flml | 22 ++++++++++--------- 4 files changed, 36 insertions(+), 33 deletions(-) diff --git a/assemble/Full_Projection.F90 b/assemble/Full_Projection.F90 index fd91fab365..1c39936bf4 100644 --- a/assemble/Full_Projection.F90 +++ b/assemble/Full_Projection.F90 @@ -77,7 +77,7 @@ subroutine petsc_solve_full_projection(x,ctp_m,inner_m,ct_m,rhs,pmat, velocity, ! as DiagonalSchurComplement, this comes in as C(Big_m_id)C, where Big_M_id is the inverse diagonal of the full ! momentum matrix. If preconditioner is set to ScaledPressureMassMatrix, this comes in as the pressure mass matrix, ! scaled by the inverse of viscosity. - type(csr_matrix), intent(inout) :: pmat + type(csr_matrix), pointer, intent(inout) :: pmat type(vector_field), intent(in) :: velocity ! used to retrieve strong diricihlet bcs ! state, and inner_mesh are used to setup mg preconditioner of inner solve type(state_type), intent(in):: state @@ -189,7 +189,7 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so ! Fluidity RHS: type(scalar_field), intent(inout) :: rhs ! Preconditioner matrix: - type(csr_matrix), intent(inout) :: preconditioner_matrix + type(csr_matrix), pointer, intent(inout) :: preconditioner_matrix ! Stabilization matrix: type(csr_matrix), optional, intent(in) :: auxiliary_matrix ! Option paths: @@ -224,7 +224,7 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so integer, dimension(:,:), pointer :: save_gnn2unn type(integer_set), dimension(velocity%dim):: boundary_row_set integer reference_node, i, rotation_stat, ref_stat - logical parallel, have_auxiliary_matrix, have_preconditioner_matrix + logical parallel, have_auxiliary_matrix logical :: apply_reference_node, apply_reference_node_from_coordinates, reference_node_owned @@ -295,7 +295,7 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so halo=div_matrix_comp%sparsity%column_halo) call allocate(petsc_numbering_p, & nnodes=block_size(div_matrix_comp,1), nfields=1, & - halo=preconditioner_matrix%sparsity%row_halo, ghost_nodes=ghost_nodes) + halo=div_matrix_comp%sparsity%row_halo, ghost_nodes=ghost_nodes) ! - why is this using the row halo of the preconditioner matrix when there might be rows missing? ! - same question about the nnodes use of the rows of the block of the divergence matrix? @@ -439,13 +439,7 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so inner_solver_option_path, petsc_numbering=petsc_numbering_u, startfromzero_in=.true.) ! leaving out petsc_numbering and mesh, so "iteration_vtus" monitor won't work! - ! FIXME: broken logic here, shouldn't depend on option path - ! Assemble preconditioner matrix in petsc format (if required): - have_preconditioner_matrix=.not.(have_option(trim(option_path)//& - "/prognostic/scheme/use_projection_method& - &/full_schur_complement/preconditioner_matrix::NoPreconditionerMatrix")) - - if(have_preconditioner_matrix) then + if(associated(preconditioner_matrix)) then pmat=csr2petsc(preconditioner_matrix, petsc_numbering_p, petsc_numbering_p) else pmat=A @@ -472,7 +466,7 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so call MatDestroy(G_t_comp,ierr) ! Destroy Compressible Divergence Operator. call MatDestroy(G_t_incomp, ierr) ! Destroy Incompressible Divergence Operator. call MatDestroy(G, ierr) ! Destroy Gradient Operator (i.e. transpose of incompressible div). - if(have_preconditioner_matrix) call MatDestroy(pmat,ierr) ! Destroy preconditioning matrix. + if(associated(preconditioner_matrix)) call MatDestroy(pmat,ierr) ! Destroy preconditioning matrix. if(have_auxiliary_matrix) call MatDestroy(S,ierr) ! Destroy stabilization matrix call deallocate( petsc_numbering_u ) diff --git a/assemble/Momentum_CG.F90 b/assemble/Momentum_CG.F90 index d08ff6c040..f1c340652d 100644 --- a/assemble/Momentum_CG.F90 +++ b/assemble/Momentum_CG.F90 @@ -2574,7 +2574,7 @@ subroutine correct_masslumped_velocity(u, inverse_masslump, ct_m, delta_p) end subroutine correct_masslumped_velocity - subroutine correct_velocity_cg(u, mass, ct_m, delta_p, state) + subroutine correct_velocity_cg(u, mass, ct_m, delta_p, state, option_path) !!< Given the pressure correction delta_p, correct the velocity. !!< !!< U_new = U_old + M_l^{-1} * C * delta_P @@ -2583,6 +2583,7 @@ subroutine correct_velocity_cg(u, mass, ct_m, delta_p, state) type(block_csr_matrix), intent(in) :: ct_m type(scalar_field), intent(in) :: delta_p type(state_type), intent(in) :: state + character(len=*), optional, intent(in) :: option_path ! Correction to u one dimension at a time. type(vector_field) :: delta_u1, delta_u2 @@ -2591,13 +2592,17 @@ subroutine correct_velocity_cg(u, mass, ct_m, delta_p, state) call allocate(delta_u1, u%dim, u%mesh, "Delta_U1") call allocate(delta_u2, u%dim, u%mesh, "Delta_U2") - delta_u2%option_path = trim(delta_p%option_path)//& - &"/prognostic/scheme/use_projection_method"//& - &"/full_schur_complement/inner_matrix[0]" - if (.not. have_option(trim(delta_u2%option_path)//"/solver")) then - ! inner solver options are optional (for FullMomemtumMatrix), if not - ! present use the same as those for the initial velocity solve - delta_u2%option_path = u%option_path + if(present(option_path)) then + delta_u2%option_path = option_path + else + delta_u2%option_path = trim(delta_p%option_path)//& + &"/prognostic/scheme/use_projection_method"//& + &"/full_schur_complement/inner_matrix[0]" + if (.not. have_option(trim(delta_u2%option_path)//"/solver")) then + ! inner solver options are optional (for FullMomemtumMatrix), if not + ! present use the same as those for the initial velocity solve + delta_u2%option_path = u%option_path + end if end if ! compute delta_u1=grad delta_p diff --git a/assemble/Solenoidal_interpolation.F90 b/assemble/Solenoidal_interpolation.F90 index 5d33b57e06..6cc37cae7f 100644 --- a/assemble/Solenoidal_interpolation.F90 +++ b/assemble/Solenoidal_interpolation.F90 @@ -416,7 +416,7 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & call petsc_solve_full_projection(lagrange, ctp_m, field_mass, ct_m, & projec_rhs, pschur, v_field, & local_state, v_field%mesh, & - option_path=trim(l_option_path)//"/lagrange_multiplier",& + option_path=trim(l_option_path)//"/lagrange_multiplier", & inner_option_path=trim(l_option_path)//& "/interpolated_field/continuous/full_schur_complement/inner_matrix[0]",& auxiliary_matrix=schur_aux) @@ -424,7 +424,7 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & call petsc_solve_full_projection(lagrange, ctp_m, field_mass, ct_m, & projec_rhs, pschur, v_field, & local_state, v_field%mesh, & - option_path=trim(l_option_path)//"/lagrange_multiplier",& + option_path=trim(l_option_path)//"/lagrange_multiplier", & inner_option_path=trim(l_option_path)//& "/interpolated_field/continuous/full_schur_complement/inner_matrix[0]") end if @@ -437,7 +437,9 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & end if if(full_schur) then - call correct_velocity_cg(v_field, field_mass, ct_m, lagrange, local_state) + call correct_velocity_cg(v_field, field_mass, ct_m, lagrange, local_state, & + option_path=trim(l_option_path)//& + "/interpolated_field/continuous/full_schur_complement/inner_matrix[0]") else if(assemble_lumped_mass) then call correct_masslumped_velocity(v_field, inverse_field_lumped_mass_vector, ct_m, lagrange) else if(discontinuous) then diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml b/tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml index e4e09da101..9b091e53f0 100644 --- a/tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml +++ b/tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml @@ -154,16 +154,15 @@ will get cut off by final_timestep - - - - + + - 1.e-6 + 1E-10 - 100 + 10000 + @@ -183,13 +182,16 @@ will get cut off by final_timestep - - + + - 1.E-10 + 1.E-8 + + 1e-16 + - 10000 + 1000 From 2ff1b3043320c2184814936b69b022d27af168be Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 8 Jul 2021 01:33:20 -0400 Subject: [PATCH 08/12] Switching from mg to gamg for failing tests. This change appears to be necessary but unrelated to other changes in this branch. --- ...s-square-convection-1e4-p1p1-parallel.flml | 24 +++++-------------- ...s-square-convection-1e4-p2p1-parallel.flml | 24 +++++-------------- 2 files changed, 12 insertions(+), 36 deletions(-) diff --git a/tests/Stokes_square_convection_1e4_p1p1_parallel/Stokes-square-convection-1e4-p1p1-parallel.flml b/tests/Stokes_square_convection_1e4_p1p1_parallel/Stokes-square-convection-1e4-p1p1-parallel.flml index 1f9c047a77..e07e02c695 100644 --- a/tests/Stokes_square_convection_1e4_p1p1_parallel/Stokes-square-convection-1e4-p1p1-parallel.flml +++ b/tests/Stokes_square_convection_1e4_p1p1_parallel/Stokes-square-convection-1e4-p1p1-parallel.flml @@ -138,23 +138,7 @@ - - - - - - 1.0e-6 - - - 10000 - - - - - - - - + @@ -225,13 +209,17 @@ - + 1.0e-6 1000 + + + + diff --git a/tests/Stokes_square_convection_1e4_p2p1_parallel/Stokes-square-convection-1e4-p2p1-parallel.flml b/tests/Stokes_square_convection_1e4_p2p1_parallel/Stokes-square-convection-1e4-p2p1-parallel.flml index 76788bc7fc..1b5bdcee3e 100644 --- a/tests/Stokes_square_convection_1e4_p2p1_parallel/Stokes-square-convection-1e4-p2p1-parallel.flml +++ b/tests/Stokes_square_convection_1e4_p2p1_parallel/Stokes-square-convection-1e4-p2p1-parallel.flml @@ -142,23 +142,7 @@ - - - - - - 1.0E-12 - - - 1000 - - - - - - - - + @@ -228,13 +212,17 @@ - + 1.0e-12 10000 + + + + From 3a48663e2bc791053dde67f0b505961542045176 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 8 Jul 2021 01:34:12 -0400 Subject: [PATCH 09/12] Add a more comprehensive test of full schur solenoidal for all pc options. diag and lump appear to be best with tolerances increased for no and mass due to poor convergence --- .../mmat-interpolation-diagpc.flml | 620 ++++++++++++++++ .../mmat-interpolation-lumppc.flml | 622 ++++++++++++++++ .../mmat-interpolation-masspc.flml | 620 ++++++++++++++++ ...tion.flml => mmat-interpolation-nopc.flml} | 20 +- .../mmat-interpolation-p2p1cv.xml | 701 +++++++++++++----- 5 files changed, 2385 insertions(+), 198 deletions(-) create mode 100644 tests/mmat-interpolation-p2p1cv/mmat-interpolation-diagpc.flml create mode 100644 tests/mmat-interpolation-p2p1cv/mmat-interpolation-lumppc.flml create mode 100644 tests/mmat-interpolation-p2p1cv/mmat-interpolation-masspc.flml rename tests/mmat-interpolation-p2p1cv/{mmat-interpolation.flml => mmat-interpolation-nopc.flml} (96%) diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-diagpc.flml b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-diagpc.flml new file mode 100644 index 0000000000..7403317f95 --- /dev/null +++ b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-diagpc.flml @@ -0,0 +1,620 @@ + + + + mmat-interpolation-diagpc + + + multimaterial + + + + 2 + + + + + + + + + + + + + + + 2 + + + + + + + + + + + + + + + + + + 5 + + + 1 + + + + + + vtk + + + + 10 + 0.36440821846106625 + + + + + + 1 + + + + + + + + 0.0 + + + 0.00427561316453 + co0p1-doubleres = 0.000893158065237 +co0p1-halfres = 0.00427561316453 + + + 14.576328738442649 + 14.576328738442649 +will get cut off by final_timestep + + + 25 + + + + 0.1 + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + def val(X,t): + from math import sin, cos + # Shear rotation about origin. + return (sin(X[0])*cos(X[1]), -1.0*cos(X[0])*sin(X[1])) + + + + + + + + + + + + + + + + + + + + + + + + + + 1E-12 + + + 100 + + + + + + + + + + + + + + + + + + + + + + + + + 1.E-6 + + + 1e-16 + + + 1000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + 0.0 + + + + 0.0 + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.2*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + 8 9 10 11 + + + + + + + + + + + + + + + + + + + + + + + + 0.11 + + + + + + + + + + + + + + + + + + + 10000 + + + + + + + 1.E-10 + + + 10000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + + + + + + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + 0.0 + + + + 0.0 + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.55*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + 8 9 10 11 + + + + + + + + + + + + + + + + + + + + + + + + 0.11 + + + + + + + + + + + + + + + + + + + 10000 + + + + + + + 1.E-10 + + + 10000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.55*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + + + + + + + + + + + + + 10 + + + 5000 + + + + + + 1.0 0.0 0.0 1.0 + + + + + + + + 0.01 0.0 0.0 0.01 + + + + + + + 1.5 0.0 0.0 1.5 + + + + + + diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-lumppc.flml b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-lumppc.flml new file mode 100644 index 0000000000..3053041be8 --- /dev/null +++ b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-lumppc.flml @@ -0,0 +1,622 @@ + + + + mmat-interpolation-lumppc + + + multimaterial + + + + 2 + + + + + + + + + + + + + + + 2 + + + + + + + + + + + + + + + + + + 5 + + + 1 + + + + + + vtk + + + + 10 + 0.36440821846106625 + + + + + + 1 + + + + + + + + 0.0 + + + 0.00427561316453 + co0p1-doubleres = 0.000893158065237 +co0p1-halfres = 0.00427561316453 + + + 14.576328738442649 + 14.576328738442649 +will get cut off by final_timestep + + + 25 + + + + 0.1 + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + def val(X,t): + from math import sin, cos + # Shear rotation about origin. + return (sin(X[0])*cos(X[1]), -1.0*cos(X[0])*sin(X[1])) + + + + + + + + + + + + + + + + + + + + + + + + + + 1E-12 + + + 100 + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.E-6 + + + 1e-16 + + + 1000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + 0.0 + + + + 0.0 + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.2*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + 8 9 10 11 + + + + + + + + + + + + + + + + + + + + + + + + 0.11 + + + + + + + + + + + + + + + + + + + 10000 + + + + + + + 1.E-10 + + + 10000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + + + + + + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + 0.0 + + + + 0.0 + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.55*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + 8 9 10 11 + + + + + + + + + + + + + + + + + + + + + + + + 0.11 + + + + + + + + + + + + + + + + + + + 10000 + + + + + + + 1.E-10 + + + 10000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.55*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + + + + + + + + + + + + + 10 + + + 5000 + + + + + + 1.0 0.0 0.0 1.0 + + + + + + + + 0.01 0.0 0.0 0.01 + + + + + + + 1.5 0.0 0.0 1.5 + + + + + + diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-masspc.flml b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-masspc.flml new file mode 100644 index 0000000000..f31f15ce41 --- /dev/null +++ b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-masspc.flml @@ -0,0 +1,620 @@ + + + + mmat-interpolation-masspc + + + multimaterial + + + + 2 + + + + + + + + + + + + + + + 2 + + + + + + + + + + + + + + + + + + 5 + + + 1 + + + + + + vtk + + + + 10 + 0.36440821846106625 + + + + + + 1 + + + + + + + + 0.0 + + + 0.00427561316453 + co0p1-doubleres = 0.000893158065237 +co0p1-halfres = 0.00427561316453 + + + 14.576328738442649 + 14.576328738442649 +will get cut off by final_timestep + + + 25 + + + + 0.1 + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + def val(X,t): + from math import sin, cos + # Shear rotation about origin. + return (sin(X[0])*cos(X[1]), -1.0*cos(X[0])*sin(X[1])) + + + + + + + + + + + + + + + + + + + + + + + + + + 1E-12 + + + 100 + + + + + + + + + + + + + + + + + + + + + + + + + 1.E-4 + + + 1e-16 + + + 1000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + 0.0 + + + + 0.0 + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.2*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + 8 9 10 11 + + + + + + + + + + + + + + + + + + + + + + + + 0.11 + + + + + + + + + + + + + + + + + + + 10000 + + + + + + + 1.E-10 + + + 10000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + + + + + + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + 0.0 + + + + + 0.0 + + + + 0.0 + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.55*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + 8 9 10 11 + + + + + + + + + + + + + + + + + + + + + + + + 0.11 + + + + + + + + + + + + + + + + + + + 10000 + + + + + + + 1.E-10 + + + 10000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + def val(X,t): + from math import sqrt, pi + dx1 = X[0]-pi/2 + dx2 = X[1]-0.55*(1.0+pi) + r=sqrt(dx1*dx1+dx2*dx2) + if (r<=(pi/5)): + return 1.0 + else: + return 0.0 + + + + + + + + + + + + + + + + + 10 + + + 5000 + + + + + + 1.0 0.0 0.0 1.0 + + + + + + + + 0.01 0.0 0.0 0.01 + + + + + + + 1.5 0.0 0.0 1.5 + + + + + + diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-nopc.flml similarity index 96% rename from tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml rename to tests/mmat-interpolation-p2p1cv/mmat-interpolation-nopc.flml index 9b091e53f0..d4361c4e96 100644 --- a/tests/mmat-interpolation-p2p1cv/mmat-interpolation.flml +++ b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-nopc.flml @@ -1,7 +1,7 @@ - mmat-interpolation + mmat-interpolation-nopc multimaterial @@ -154,13 +154,15 @@ will get cut off by final_timestep - - + + + + - 1E-10 + 1E-12 - 10000 + 100 @@ -182,10 +184,10 @@ will get cut off by final_timestep - + - 1.E-8 + 1.E-4 1e-16 @@ -195,7 +197,9 @@ will get cut off by final_timestep - + + + diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml index 23be75ccb3..a6dbfe8436 100644 --- a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml +++ b/tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml @@ -4,7 +4,11 @@ flml 2dadapt - fluidity -v2 -l mmat-interpolation.flml && mv fluidity.log-0 first.log && mv fluidity.err-0 first.err && ./change_options.py && fluidity -v2 -l mmat-interpolation_2_checkpoint.flml && mv fluidity.log-0 checkpoint.log && mv fluidity.err-0 checkpoint.err + fluidity -v2 -l mmat-interpolation-nopc.flml && mv fluidity.log-0 nopc.log && mv fluidity.err-0 nopc.err && + fluidity -v2 -l mmat-interpolation-masspc.flml && mv fluidity.log-0 masspc.log && mv fluidity.err-0 masspc.err && + fluidity -v2 -l mmat-interpolation-lumppc.flml && mv fluidity.log-0 lumppc.log && mv fluidity.err-0 lumppc.err && + fluidity -v2 -l mmat-interpolation-diagpc.flml && mv fluidity.log-0 diagpc.log && mv fluidity.err-0 diagpc.err + @@ -12,339 +16,656 @@ import os files = os.listdir("./") solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files - + from fluidity_tools import stat_parser as stat -material1integralstart = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] +nopc_material1integralstart = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat -material1integralend = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] +nopc_material1integralend = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat -material1maxstart = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] +nopc_material1maxstart = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat -material1maxend = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] +nopc_material1maxend = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat -material1minstart = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] +nopc_material1minstart = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat -material1minend = stat("mmat-interpolation.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] +nopc_material1minend = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat -material2integralstart = stat("mmat-interpolation.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] +nopc_material2integralstart = stat("mmat-interpolation-nopc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat -material2integralend = stat("mmat-interpolation.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] +nopc_material2integralend = stat("mmat-interpolation-nopc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat -material3integralstart = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] +nopc_material3integralstart = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat -material3integralend = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] +nopc_material3integralend = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat -material3maxstart = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] +nopc_material3maxstart = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat -material3maxend = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] +nopc_material3maxend = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat -material3minstart = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] +nopc_material3minstart = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat -material3minend = stat("mmat-interpolation.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] +nopc_material3minend = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat -divergenceminstart = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] +nopc_divergenceminstart = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] - + from fluidity_tools import stat_parser as stat -divergenceminend = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] +nopc_divergenceminend = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] - + from fluidity_tools import stat_parser as stat -divergencemaxstart = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] +nopc_divergencemaxstart = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] - + from fluidity_tools import stat_parser as stat -divergencemaxend = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] +nopc_divergencemaxend = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] - + from fluidity_tools import stat_parser as stat -cflmaxstart = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] +nopc_cflmaxstart = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] - + from fluidity_tools import stat_parser as stat -cflmaxend = stat("mmat-interpolation.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] +nopc_cflmaxend = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] - + from fluidity_tools import stat_parser as stat -lambdamaxstart = stat("mmat-interpolation.stat")["Material1"]["Lambda"]["max"][0] +nopc_lambdamaxstart = stat("mmat-interpolation-nopc.stat")["Material1"]["Lambda"]["max"][0] - + from fluidity_tools import stat_parser as stat -lambdamaxend = stat("mmat-interpolation.stat")["Material1"]["Lambda"]["max"][-1] +nopc_lambdamaxend = stat("mmat-interpolation-nopc.stat")["Material1"]["Lambda"]["max"][-1] - + from fluidity_tools import stat_parser as stat -lambdaminstart = stat("mmat-interpolation.stat")["Material1"]["Lambda"]["min"][0] +nopc_lambdaminstart = stat("mmat-interpolation-nopc.stat")["Material1"]["Lambda"]["min"][0] - + from fluidity_tools import stat_parser as stat -lambdaminend = stat("mmat-interpolation.stat")["Material1"]["Lambda"]["min"][-1] +nopc_lambdaminend = stat("mmat-interpolation-nopc.stat")["Material1"]["Lambda"]["min"][-1] - + from fluidity_tools import stat_parser as stat -material1integralstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] +nopc_absolutedifference = stat("mmat-interpolation-nopc.stat")["Material3"]["ScalarAbsoluteDifference"]["max"][1] - + from fluidity_tools import stat_parser as stat -material1integralend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] +masspc_material1integralstart = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat -material1maxstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] +masspc_material1integralend = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat -material1maxend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] +masspc_material1maxstart = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat -material1minstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] +masspc_material1maxend = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat -material1minend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] +masspc_material1minstart = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat -material2integralstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] +masspc_material1minend = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat -material2integralend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] +masspc_material2integralstart = stat("mmat-interpolation-masspc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat -material3integralstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] +masspc_material2integralend = stat("mmat-interpolation-masspc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat -material3integralend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] +masspc_material3integralstart = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat -material3maxstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] +masspc_material3integralend = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat -material3maxend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] +masspc_material3maxstart = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat -material3minstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] +masspc_material3maxend = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat -material3minend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] +masspc_material3minstart = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat -divergenceminstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] +masspc_material3minend = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat -divergenceminend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] +masspc_divergenceminstart = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] - + from fluidity_tools import stat_parser as stat -divergencemaxstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] +masspc_divergenceminend = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] - + from fluidity_tools import stat_parser as stat -divergencemaxend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] +masspc_divergencemaxstart = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] - + from fluidity_tools import stat_parser as stat -cflmaxstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] +masspc_divergencemaxend = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] - + from fluidity_tools import stat_parser as stat -cflmaxend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] +masspc_cflmaxstart = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] - + from fluidity_tools import stat_parser as stat -lambdamaxstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["Lambda"]["max"][0] +masspc_cflmaxend = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] - + from fluidity_tools import stat_parser as stat -lambdamaxend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["Lambda"]["max"][-1] +masspc_lambdamaxstart = stat("mmat-interpolation-masspc.stat")["Material1"]["Lambda"]["max"][0] - + from fluidity_tools import stat_parser as stat -lambdaminstart_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["Lambda"]["min"][0] +masspc_lambdamaxend = stat("mmat-interpolation-masspc.stat")["Material1"]["Lambda"]["max"][-1] - + from fluidity_tools import stat_parser as stat -lambdaminend_checkpoint = stat("mmat-interpolation_checkpoint.stat")["Material1"]["Lambda"]["min"][-1] +masspc_lambdaminstart = stat("mmat-interpolation-masspc.stat")["Material1"]["Lambda"]["min"][0] - + from fluidity_tools import stat_parser as stat -absolutedifference = stat("mmat-interpolation.stat")["Material3"]["ScalarAbsoluteDifference"]["max"][1] +masspc_lambdaminend = stat("mmat-interpolation-masspc.stat")["Material1"]["Lambda"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +masspc_absolutedifference = stat("mmat-interpolation-masspc.stat")["Material3"]["ScalarAbsoluteDifference"]["max"][1] + + +from fluidity_tools import stat_parser as stat +lumppc_material1integralstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_material1integralend = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_material1maxstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_material1maxend = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_material1minstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_material1minend = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_material2integralstart = stat("mmat-interpolation-lumppc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_material2integralend = stat("mmat-interpolation-lumppc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_material3integralstart = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_material3integralend = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_material3maxstart = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_material3maxend = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_material3minstart = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_material3minend = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_divergenceminstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_divergenceminend = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_divergencemaxstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_divergencemaxend = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_cflmaxstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_cflmaxend = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_lambdamaxstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["Lambda"]["max"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_lambdamaxend = stat("mmat-interpolation-lumppc.stat")["Material1"]["Lambda"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_lambdaminstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["Lambda"]["min"][0] + + +from fluidity_tools import stat_parser as stat +lumppc_lambdaminend = stat("mmat-interpolation-lumppc.stat")["Material1"]["Lambda"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +lumppc_absolutedifference = stat("mmat-interpolation-lumppc.stat")["Material3"]["ScalarAbsoluteDifference"]["max"][1] + + +from fluidity_tools import stat_parser as stat +diagpc_material1integralstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_material1integralend = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_material1maxstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_material1maxend = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_material1minstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_material1minend = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_material2integralstart = stat("mmat-interpolation-diagpc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_material2integralend = stat("mmat-interpolation-diagpc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_material3integralstart = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_material3integralend = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_material3maxstart = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_material3maxend = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_material3minstart = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_material3minend = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_divergenceminstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_divergenceminend = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_divergencemaxstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_divergencemaxend = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_cflmaxstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_cflmaxend = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_lambdamaxstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["Lambda"]["max"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_lambdamaxend = stat("mmat-interpolation-diagpc.stat")["Material1"]["Lambda"]["max"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_lambdaminstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["Lambda"]["min"][0] + + +from fluidity_tools import stat_parser as stat +diagpc_lambdaminend = stat("mmat-interpolation-diagpc.stat")["Material1"]["Lambda"]["min"][-1] + + +from fluidity_tools import stat_parser as stat +diagpc_absolutedifference = stat("mmat-interpolation-diagpc.stat")["Material3"]["ScalarAbsoluteDifference"]["max"][1] assert(solvers_converged) - -print('mass loss = ', abs(material1integralstart-material1integralend)) -assert abs(material1integralstart-material1integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + +print('mass loss = ', abs(nopc_material1integralstart-nopc_material1integralend)) +assert abs(nopc_material1integralstart-nopc_material1integralend) < max(abs(nopc_divergencemaxstart), abs(nopc_divergencemaxend), abs(nopc_divergenceminstart), abs(nopc_divergenceminend), 1.E-8) - -print('mass loss = ', abs(material2integralstart-material2integralend)) -assert abs(material2integralstart-material2integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + +print('mass loss = ', abs(nopc_material2integralstart-nopc_material2integralend)) +assert abs(nopc_material2integralstart-nopc_material2integralend) < max(abs(nopc_divergencemaxstart), abs(nopc_divergencemaxend), abs(nopc_divergenceminstart), abs(nopc_divergenceminend), 1.E-8) - -print('mass loss = ', abs(material3integralstart-material3integralend)) -assert abs(material3integralstart-material3integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + +print('mass loss = ', abs(nopc_material3integralstart-nopc_material3integralend)) +assert abs(nopc_material3integralstart-nopc_material3integralend) < max(abs(nopc_divergencemaxstart), abs(nopc_divergencemaxend), abs(nopc_divergenceminstart), abs(nopc_divergenceminend), 1.E-8) - -print('divergence tolerance = ', max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend))) -assert max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend)) < 1.E-8 + +print('divergence tolerance = ', max(abs(nopc_divergencemaxstart), abs(nopc_divergencemaxend), abs(nopc_divergenceminstart), abs(nopc_divergenceminend))) +assert max(abs(nopc_divergencemaxstart), abs(nopc_divergencemaxend), abs(nopc_divergenceminstart), abs(nopc_divergenceminend)) < 1.E-4 - -assert abs(material1maxstart-material1maxend) < 1.E-10 + +assert abs(nopc_material1maxstart-nopc_material1maxend) < 1.E-10 - -assert abs(material3maxstart-material3maxend) < 1.E-10 + +assert abs(nopc_material3maxstart-nopc_material3maxend) < 1.E-10 - -assert abs(material1minstart-material1minend) < 1.E-10 + +assert abs(nopc_material1minstart-nopc_material1minend) < 1.E-10 - -assert abs(material3minstart-material3minend) < 1.E-10 + +assert abs(nopc_material3minstart-nopc_material3minend) < 1.E-10 - -assert abs(divergencemaxstart) < 1.E-8 + +assert abs(nopc_divergencemaxstart) < 1.E-4 - -assert abs(divergenceminstart) < 1.E-8 + +assert abs(nopc_divergenceminstart) < 1.E-4 - -assert abs(divergencemaxend) < 1.E-8 + +assert abs(nopc_divergencemaxend) < 1.E-4 - -assert abs(divergenceminend) < 1.E-8 + +assert abs(nopc_divergenceminend) < 1.E-4 - -assert abs(cflmaxstart-0.1) < 1.E-10 + +assert abs(nopc_cflmaxstart-0.1) < 1.E-10 - -assert abs(cflmaxend-0.1) < 1.E-10 + +assert abs(nopc_cflmaxend-0.1) < 1.E-10 - -print('mass loss = ', abs(material1integralstart_checkpoint-material1integralend_checkpoint)) -assert abs(material1integralstart_checkpoint-material1integralend_checkpoint) < max(abs(divergencemaxstart_checkpoint), abs(divergencemaxend_checkpoint), abs(divergenceminstart_checkpoint), abs(divergenceminend_checkpoint), 1.E-8) + +print('mass loss = ', abs(masspc_material1integralstart-masspc_material1integralend)) +assert abs(masspc_material1integralstart-masspc_material1integralend) < max(abs(masspc_divergencemaxstart), abs(masspc_divergencemaxend), abs(masspc_divergenceminstart), abs(masspc_divergenceminend), 1.E-8) - -print('mass loss = ', abs(material2integralstart_checkpoint-material2integralend_checkpoint)) -assert abs(material2integralstart_checkpoint-material2integralend_checkpoint) < max(abs(divergencemaxstart_checkpoint), abs(divergencemaxend_checkpoint), abs(divergenceminstart_checkpoint), abs(divergenceminend_checkpoint), 1.E-8) + +print('mass loss = ', abs(masspc_material2integralstart-masspc_material2integralend)) +assert abs(masspc_material2integralstart-masspc_material2integralend) < max(abs(masspc_divergencemaxstart), abs(masspc_divergencemaxend), abs(masspc_divergenceminstart), abs(masspc_divergenceminend), 1.E-8) - -print('mass loss = ', abs(material3integralstart_checkpoint-material3integralend_checkpoint)) -assert abs(material3integralstart_checkpoint-material3integralend_checkpoint) < max(abs(divergencemaxstart_checkpoint), abs(divergencemaxend_checkpoint), abs(divergenceminstart_checkpoint), abs(divergenceminend_checkpoint), 1.E-8) + +print('mass loss = ', abs(masspc_material3integralstart-masspc_material3integralend)) +assert abs(masspc_material3integralstart-masspc_material3integralend) < max(abs(masspc_divergencemaxstart), abs(masspc_divergencemaxend), abs(masspc_divergenceminstart), abs(masspc_divergenceminend), 1.E-8) - -print('divergence tolerance = ', max(abs(divergencemaxstart_checkpoint), abs(divergencemaxend_checkpoint), abs(divergenceminstart_checkpoint), abs(divergenceminend_checkpoint))) -assert max(abs(divergencemaxstart_checkpoint), abs(divergencemaxend_checkpoint), abs(divergenceminstart_checkpoint), abs(divergenceminend_checkpoint)) < 1.E-8 + +print('divergence tolerance = ', max(abs(masspc_divergencemaxstart), abs(masspc_divergencemaxend), abs(masspc_divergenceminstart), abs(masspc_divergenceminend))) +assert max(abs(masspc_divergencemaxstart), abs(masspc_divergencemaxend), abs(masspc_divergenceminstart), abs(masspc_divergenceminend)) < 1.E-4 - -assert abs(material1maxstart_checkpoint-material1maxend_checkpoint) < 1.E-10 + +assert abs(masspc_material1maxstart-masspc_material1maxend) < 1.E-10 - -assert abs(material3maxstart_checkpoint-material3maxend_checkpoint) < 1.E-10 + +assert abs(masspc_material3maxstart-masspc_material3maxend) < 1.E-10 - -assert abs(material1minstart_checkpoint-material1minend_checkpoint) < 1.E-10 + +assert abs(masspc_material1minstart-masspc_material1minend) < 1.E-10 - -assert abs(material3minstart_checkpoint-material3minend_checkpoint) < 1.E-10 + +assert abs(masspc_material3minstart-masspc_material3minend) < 1.E-10 - -assert abs(divergencemaxstart_checkpoint) < 1.E-8 + +assert abs(masspc_divergencemaxstart) < 1.E-4 - -assert abs(divergenceminstart_checkpoint) < 1.E-8 + +assert abs(masspc_divergenceminstart) < 1.E-4 - -assert abs(divergencemaxend_checkpoint) < 1.E-8 + +assert abs(masspc_divergencemaxend) < 1.E-4 - -assert abs(divergenceminend_checkpoint) < 1.E-8 + +assert abs(masspc_divergenceminend) < 1.E-4 - -assert abs(cflmaxstart_checkpoint-0.1) < 1.E-10 + +assert abs(masspc_cflmaxstart-0.1) < 1.E-10 - -assert abs(cflmaxend_checkpoint-0.1) < 1.E-10 + +assert abs(masspc_cflmaxend-0.1) < 1.E-10 - -print('mass loss = ', abs(material1integralstart_checkpoint-material1integralend)) -assert abs(material1integralstart_checkpoint-material1integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + +print('mass loss = ', abs(lumppc_material1integralstart-lumppc_material1integralend)) +assert abs(lumppc_material1integralstart-lumppc_material1integralend) < max(abs(lumppc_divergencemaxstart), abs(lumppc_divergencemaxend), abs(lumppc_divergenceminstart), abs(lumppc_divergenceminend), 1.E-8) - -print('mass loss = ', abs(material2integralstart_checkpoint-material2integralend)) -assert abs(material2integralstart_checkpoint-material2integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + +print('mass loss = ', abs(lumppc_material2integralstart-lumppc_material2integralend)) +assert abs(lumppc_material2integralstart-lumppc_material2integralend) < max(abs(lumppc_divergencemaxstart), abs(lumppc_divergencemaxend), abs(lumppc_divergenceminstart), abs(lumppc_divergenceminend), 1.E-8) - -print('mass loss = ', abs(material3integralstart_checkpoint-material3integralend)) -assert abs(material3integralstart_checkpoint-material3integralend) < max(abs(divergencemaxstart), abs(divergencemaxend), abs(divergenceminstart), abs(divergenceminend), 1.E-8) + +print('mass loss = ', abs(lumppc_material3integralstart-lumppc_material3integralend)) +assert abs(lumppc_material3integralstart-lumppc_material3integralend) < max(abs(lumppc_divergencemaxstart), abs(lumppc_divergencemaxend), abs(lumppc_divergenceminstart), abs(lumppc_divergenceminend), 1.E-8) + +print('divergence tolerance = ', max(abs(lumppc_divergencemaxstart), abs(lumppc_divergencemaxend), abs(lumppc_divergenceminstart), abs(lumppc_divergenceminend))) +assert max(abs(lumppc_divergencemaxstart), abs(lumppc_divergencemaxend), abs(lumppc_divergenceminstart), abs(lumppc_divergenceminend)) < 1.E-6 + + +assert abs(lumppc_material1maxstart-lumppc_material1maxend) < 1.E-10 + + +assert abs(lumppc_material3maxstart-lumppc_material3maxend) < 1.E-10 + + +assert abs(lumppc_material1minstart-lumppc_material1minend) < 1.E-10 + + +assert abs(lumppc_material3minstart-lumppc_material3minend) < 1.E-10 + + +assert abs(lumppc_divergencemaxstart) < 1.E-6 + + +assert abs(lumppc_divergenceminstart) < 1.E-6 + + +assert abs(lumppc_divergencemaxend) < 1.E-6 + + +assert abs(lumppc_divergenceminend) < 1.E-6 + + +assert abs(lumppc_cflmaxstart-0.1) < 1.E-10 + + +assert abs(lumppc_cflmaxend-0.1) < 1.E-10 + + +print('mass loss = ', abs(diagpc_material1integralstart-diagpc_material1integralend)) +assert abs(diagpc_material1integralstart-diagpc_material1integralend) < max(abs(diagpc_divergencemaxstart), abs(diagpc_divergencemaxend), abs(diagpc_divergenceminstart), abs(diagpc_divergenceminend), 1.E-8) + + +print('mass loss = ', abs(diagpc_material2integralstart-diagpc_material2integralend)) +assert abs(diagpc_material2integralstart-diagpc_material2integralend) < max(abs(diagpc_divergencemaxstart), abs(diagpc_divergencemaxend), abs(diagpc_divergenceminstart), abs(diagpc_divergenceminend), 1.E-8) + + +print('mass loss = ', abs(diagpc_material3integralstart-diagpc_material3integralend)) +assert abs(diagpc_material3integralstart-diagpc_material3integralend) < max(abs(diagpc_divergencemaxstart), abs(diagpc_divergencemaxend), abs(diagpc_divergenceminstart), abs(diagpc_divergenceminend), 1.E-8) + + +print('divergence tolerance = ', max(abs(diagpc_divergencemaxstart), abs(diagpc_divergencemaxend), abs(diagpc_divergenceminstart), abs(diagpc_divergenceminend))) +assert max(abs(diagpc_divergencemaxstart), abs(diagpc_divergencemaxend), abs(diagpc_divergenceminstart), abs(diagpc_divergenceminend)) < 1.E-6 + + +assert abs(diagpc_material1maxstart-diagpc_material1maxend) < 1.E-10 + + +assert abs(diagpc_material3maxstart-diagpc_material3maxend) < 1.E-10 + + +assert abs(diagpc_material1minstart-diagpc_material1minend) < 1.E-10 + + +assert abs(diagpc_material3minstart-diagpc_material3minend) < 1.E-10 + + +assert abs(diagpc_divergencemaxstart) < 1.E-6 + + +assert abs(diagpc_divergenceminstart) < 1.E-6 + + +assert abs(diagpc_divergencemaxend) < 1.E-6 + + +assert abs(diagpc_divergenceminend) < 1.E-6 + + +assert abs(diagpc_cflmaxstart-0.1) < 1.E-10 + + +assert abs(diagpc_cflmaxend-0.1) < 1.E-10 + - -assert abs(lambdamaxstart) > 1.E-7 + +assert abs(nopc_lambdamaxstart) > 1.E-7 + + +assert abs(nopc_lambdaminstart) > 1.E-7 + + +assert abs(nopc_lambdamaxend) > 1.E-7 + + +assert abs(nopc_lambdaminend) > 1.E-7 + + +assert abs(nopc_absolutedifference) > 1.E-7 + + +assert abs(masspc_lambdamaxstart) > 1.E-7 + + +assert abs(masspc_lambdaminstart) > 1.E-7 + + +assert abs(masspc_lambdamaxend) > 1.E-7 + + +assert abs(masspc_lambdaminend) > 1.E-7 + + +assert abs(masspc_absolutedifference) > 1.E-7 + + +assert abs(lumppc_lambdamaxstart) > 1.E-7 + + +assert abs(lumppc_lambdaminstart) > 1.E-7 - -assert abs(lambdaminstart) > 1.E-7 + +assert abs(lumppc_lambdamaxend) > 1.E-7 - -assert abs(lambdamaxend) > 1.E-7 + +assert abs(lumppc_lambdaminend) > 1.E-7 - -assert abs(lambdaminend) > 1.E-7 + +assert abs(lumppc_absolutedifference) > 1.E-7 - -assert abs(lambdamaxstart_checkpoint) < 1.E-7 + +assert abs(diagpc_lambdamaxstart) > 1.E-7 - -assert abs(lambdaminstart_checkpoint) < 1.E-7 + +assert abs(diagpc_lambdaminstart) > 1.E-7 - -assert abs(lambdamaxend_checkpoint) > 1.E-7 + +assert abs(diagpc_lambdamaxend) > 1.E-7 - -assert abs(lambdaminend_checkpoint) > 1.E-7 + +assert abs(diagpc_lambdaminend) > 1.E-7 - -assert abs(absolutedifference) > 1.E-7 + +assert abs(diagpc_absolutedifference) > 1.E-7 From fb15c9efc250e530df6de6a9bdea651615998d5d Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 8 Jul 2021 12:29:32 -0400 Subject: [PATCH 10/12] Use the preconditioner matrix if it's available for the pressure petsc numbering so that we get a second order halo when it's required. --- assemble/Full_Projection.F90 | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/assemble/Full_Projection.F90 b/assemble/Full_Projection.F90 index 1c39936bf4..d6d63bd17c 100644 --- a/assemble/Full_Projection.F90 +++ b/assemble/Full_Projection.F90 @@ -293,17 +293,18 @@ subroutine petsc_solve_setup_full_projection(y,A,b,ksp,petsc_numbering_p,name,so nnodes=block_size(div_matrix_comp,2), nfields=blocks(div_matrix_comp,2), & group_size=inner_m%row_numbering%group_size, & halo=div_matrix_comp%sparsity%column_halo) - call allocate(petsc_numbering_p, & - nnodes=block_size(div_matrix_comp,1), nfields=1, & - halo=div_matrix_comp%sparsity%row_halo, ghost_nodes=ghost_nodes) - - ! - why is this using the row halo of the preconditioner matrix when there might be rows missing? - ! - same question about the nnodes use of the rows of the block of the divergence matrix? - ! - and how can ghost_nodes be appropriate for both this and the auxiliary_matrix? - ! this definitely appears to be inappropriate for the auxiliary matrix (hence there's a new one added - ! below) so two questions: - ! 1. is it definitely appropriate for all its other used (the divergence matrix and the pressure vectors)? - ! 2. can it be made appropriate for the auxiliary matrix at the same time as being appropriate for the current uses? + ! the preconditioner matrix might have a second order halo, which it will need if it's present + ! so let's use that if it's available to set up the petsc numbering + if(associated(preconditioner_matrix)) then + call allocate(petsc_numbering_p, & + nnodes=block_size(div_matrix_comp,1), nfields=1, & + halo=preconditioner_matrix%sparsity%row_halo, ghost_nodes=ghost_nodes) + else + ! if it's not available though we should default to the first order halo from the div matrix + call allocate(petsc_numbering_p, & + nnodes=block_size(div_matrix_comp,1), nfields=1, & + halo=div_matrix_comp%sparsity%row_halo, ghost_nodes=ghost_nodes) + end if ! the rows of the gradient matrix (ct_m^T) and columns of ctp_m ! corresponding to dirichlet bcs have not been zeroed From d0f1b07a508508e930b69001edde563c9910010e Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 8 Jul 2021 12:30:53 -0400 Subject: [PATCH 11/12] Make this test parallel to catch solenoidal halo problems and rename it appropriately. --- tests/mmat-interpolation-p2p1cv/Makefile | 7 ---- .../change_options.py | 34 ------------------- .../Makefile | 7 ++++ .../mmat-interpolation-diagpc.flml | 0 .../mmat-interpolation-lumppc.flml | 0 .../mmat-interpolation-masspc.flml | 0 .../mmat-interpolation-nopc.flml | 0 .../mmat-interpolation-parallel-p2p1cv.xml} | 14 +++++--- .../src/2d_square.geo | 0 9 files changed, 16 insertions(+), 46 deletions(-) delete mode 100644 tests/mmat-interpolation-p2p1cv/Makefile delete mode 100755 tests/mmat-interpolation-p2p1cv/change_options.py create mode 100644 tests/mmat-interpolation-parallel-p2p1cv/Makefile rename tests/{mmat-interpolation-p2p1cv => mmat-interpolation-parallel-p2p1cv}/mmat-interpolation-diagpc.flml (100%) rename tests/{mmat-interpolation-p2p1cv => mmat-interpolation-parallel-p2p1cv}/mmat-interpolation-lumppc.flml (100%) rename tests/{mmat-interpolation-p2p1cv => mmat-interpolation-parallel-p2p1cv}/mmat-interpolation-masspc.flml (100%) rename tests/{mmat-interpolation-p2p1cv => mmat-interpolation-parallel-p2p1cv}/mmat-interpolation-nopc.flml (100%) rename tests/{mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml => mmat-interpolation-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml} (96%) rename tests/{mmat-interpolation-p2p1cv => mmat-interpolation-parallel-p2p1cv}/src/2d_square.geo (100%) diff --git a/tests/mmat-interpolation-p2p1cv/Makefile b/tests/mmat-interpolation-p2p1cv/Makefile deleted file mode 100644 index b2e7c898e8..0000000000 --- a/tests/mmat-interpolation-p2p1cv/Makefile +++ /dev/null @@ -1,7 +0,0 @@ -input: clean - gmsh -2 src/2d_square.geo - cp src/2d_square.msh . - -clean: - rm -rf *.err *.log *.stat *.vtu *.node *.ele *.edge *checkpoint* *convergence* *.log-0 *.err-0 src/2d_square.msh *.msh \ - matrixdump matrixdump.info diff --git a/tests/mmat-interpolation-p2p1cv/change_options.py b/tests/mmat-interpolation-p2p1cv/change_options.py deleted file mode 100755 index 0224050f6b..0000000000 --- a/tests/mmat-interpolation-p2p1cv/change_options.py +++ /dev/null @@ -1,34 +0,0 @@ -#!/usr/bin/env python3 - -import sys - -file_name = ["mmat-interpolation_2_checkpoint.flml"] - -for f in range(len(file_name)): - try: - flml_file = open(file_name[f], "r") - flml_options = flml_file.readlines() - flml_file.close() - except: - sys.stderr.write("Error: change_options failed to read options file\n") - sys.exit(1) - - for i in range(len(flml_options)): - if "" in flml_options[i]: - while not "" in flml_options[i] and i < len(flml_options): - flml_options[i] = "" - i += 1 - if i < len(flml_options): - flml_options[i] = "" - break - - try: - flml_file = open(file_name[f], "w") - flml_file.writelines(flml_options) - flml_file.close() - except: - sys.stderr.write("Error: change_options failed to write options file\n") - sys.exit(1) - - - diff --git a/tests/mmat-interpolation-parallel-p2p1cv/Makefile b/tests/mmat-interpolation-parallel-p2p1cv/Makefile new file mode 100644 index 0000000000..b2c5f4305c --- /dev/null +++ b/tests/mmat-interpolation-parallel-p2p1cv/Makefile @@ -0,0 +1,7 @@ +input: clean + gmsh -2 src/2d_square.geo + cp src/2d_square.msh . + +clean: + rm -rf mmat-interpolation-*_* *-flredecomp.flml *.err-* *.log-* *.stat *.halo *.[p]vtu *.node *.ele *.edge *checkpoint* *convergence* src/2d_square.msh *.msh \ + matrixdump matrixdump.info diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-diagpc.flml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-diagpc.flml similarity index 100% rename from tests/mmat-interpolation-p2p1cv/mmat-interpolation-diagpc.flml rename to tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-diagpc.flml diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-lumppc.flml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-lumppc.flml similarity index 100% rename from tests/mmat-interpolation-p2p1cv/mmat-interpolation-lumppc.flml rename to tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-lumppc.flml diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-masspc.flml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-masspc.flml similarity index 100% rename from tests/mmat-interpolation-p2p1cv/mmat-interpolation-masspc.flml rename to tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-masspc.flml diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-nopc.flml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-nopc.flml similarity index 100% rename from tests/mmat-interpolation-p2p1cv/mmat-interpolation-nopc.flml rename to tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-nopc.flml diff --git a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml similarity index 96% rename from tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml rename to tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml index a6dbfe8436..cc13b67f3f 100644 --- a/tests/mmat-interpolation-p2p1cv/mmat-interpolation-p2p1cv.xml +++ b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml @@ -3,11 +3,15 @@ Multimaterial adaptivity and interpolation test flml 2dadapt - - fluidity -v2 -l mmat-interpolation-nopc.flml && mv fluidity.log-0 nopc.log && mv fluidity.err-0 nopc.err && - fluidity -v2 -l mmat-interpolation-masspc.flml && mv fluidity.log-0 masspc.log && mv fluidity.err-0 masspc.err && - fluidity -v2 -l mmat-interpolation-lumppc.flml && mv fluidity.log-0 lumppc.log && mv fluidity.err-0 lumppc.err && - fluidity -v2 -l mmat-interpolation-diagpc.flml && mv fluidity.log-0 diagpc.log && mv fluidity.err-0 diagpc.err + + mpiexec ../../bin/flredecomp -i 1 -o 2 -v -l mmat-interpolation-nopc mmat-interpolation-nopc-flredecomp && + mpiexec ../../bin/flredecomp -i 1 -o 2 -v -l mmat-interpolation-masspc mmat-interpolation-masspc-flredecomp && + mpiexec ../../bin/flredecomp -i 1 -o 2 -v -l mmat-interpolation-lumppc mmat-interpolation-lumppc-flredecomp && + mpiexec ../../bin/flredecomp -i 1 -o 2 -v -l mmat-interpolation-diagpc mmat-interpolation-diagpc-flredecomp && + mpiexec ../../bin/fluidity -v2 -l mmat-interpolation-nopc-flredecomp.flml && mv fluidity.log-0 nopc.log-0 && mv fluidity.err-0 nopc.err-0 && + mpiexec ../../bin/fluidity -v2 -l mmat-interpolation-masspc-flredecomp.flml && mv fluidity.log-0 masspc.log-0 && mv fluidity.err-0 masspc.err-0 && + mpiexec ../../bin/fluidity -v2 -l mmat-interpolation-lumppc-flredecomp.flml && mv fluidity.log-0 lumppc.log-0 && mv fluidity.err-0 lumppc.err-0 && + mpiexec ../../bin/fluidity -v2 -l mmat-interpolation-diagpc-flredecomp.flml && mv fluidity.log-0 diagpc.log-0 && mv fluidity.err-0 diagpc.err-0 diff --git a/tests/mmat-interpolation-p2p1cv/src/2d_square.geo b/tests/mmat-interpolation-parallel-p2p1cv/src/2d_square.geo similarity index 100% rename from tests/mmat-interpolation-p2p1cv/src/2d_square.geo rename to tests/mmat-interpolation-parallel-p2p1cv/src/2d_square.geo From 279e8c6a20615125cac9b050a2740ad3f7d91883 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 4 Nov 2022 01:37:51 +0000 Subject: [PATCH 12/12] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- schemas/fluidity_options.rnc | 2 +- schemas/solvers.rnc | 2 +- .../mmat-interpolation-parallel-p2p1cv.xml | 204 +++++++++--------- 3 files changed, 104 insertions(+), 104 deletions(-) diff --git a/schemas/fluidity_options.rnc b/schemas/fluidity_options.rnc index c0864b9887..c0f1d131e1 100644 --- a/schemas/fluidity_options.rnc +++ b/schemas/fluidity_options.rnc @@ -7502,7 +7502,7 @@ solenoidal_options = ( ## Solver options for the (inner) solve. element full_schur_complement { - ## Specify the inner matrix (IM) to form the projection schur complement (C^T*IM^{-1}*C). + ## Specify the inner matrix (IM) to form the projection schur complement (C^T*IM^{-1}*C). ## Use the full mass matrix (the only option). element inner_matrix { attribute name { "FullMassMatrix" }, diff --git a/schemas/solvers.rnc b/schemas/solvers.rnc index 0db8c56869..becc4978d4 100644 --- a/schemas/solvers.rnc +++ b/schemas/solvers.rnc @@ -47,7 +47,7 @@ linear_solver_options_sym_scalar = kspgmres_options| ksppreonly_options| ksprichardson_options| - kspother_options + kspother_options ), ## Preconditioner to be used in combination with the iterative method. ( diff --git a/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml index cc13b67f3f..fc39e2a224 100644 --- a/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml +++ b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml @@ -8,10 +8,10 @@ mpiexec ../../bin/flredecomp -i 1 -o 2 -v -l mmat-interpolation-masspc mmat-interpolation-masspc-flredecomp && mpiexec ../../bin/flredecomp -i 1 -o 2 -v -l mmat-interpolation-lumppc mmat-interpolation-lumppc-flredecomp && mpiexec ../../bin/flredecomp -i 1 -o 2 -v -l mmat-interpolation-diagpc mmat-interpolation-diagpc-flredecomp && - mpiexec ../../bin/fluidity -v2 -l mmat-interpolation-nopc-flredecomp.flml && mv fluidity.log-0 nopc.log-0 && mv fluidity.err-0 nopc.err-0 && + mpiexec ../../bin/fluidity -v2 -l mmat-interpolation-nopc-flredecomp.flml && mv fluidity.log-0 nopc.log-0 && mv fluidity.err-0 nopc.err-0 && mpiexec ../../bin/fluidity -v2 -l mmat-interpolation-masspc-flredecomp.flml && mv fluidity.log-0 masspc.log-0 && mv fluidity.err-0 masspc.err-0 && mpiexec ../../bin/fluidity -v2 -l mmat-interpolation-lumppc-flredecomp.flml && mv fluidity.log-0 lumppc.log-0 && mv fluidity.err-0 lumppc.err-0 && - mpiexec ../../bin/fluidity -v2 -l mmat-interpolation-diagpc-flredecomp.flml && mv fluidity.log-0 diagpc.log-0 && mv fluidity.err-0 diagpc.err-0 + mpiexec ../../bin/fluidity -v2 -l mmat-interpolation-diagpc-flredecomp.flml && mv fluidity.log-0 diagpc.log-0 && mv fluidity.err-0 diagpc.err-0 @@ -20,403 +20,403 @@ import os files = os.listdir("./") solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files - + from fluidity_tools import stat_parser as stat nopc_material1integralstart = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat nopc_material1integralend = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat nopc_material1maxstart = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat nopc_material1maxend = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat nopc_material1minstart = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat nopc_material1minend = stat("mmat-interpolation-nopc.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat nopc_material2integralstart = stat("mmat-interpolation-nopc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat nopc_material2integralend = stat("mmat-interpolation-nopc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat nopc_material3integralstart = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat nopc_material3integralend = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat nopc_material3maxstart = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat nopc_material3maxend = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat nopc_material3minstart = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat nopc_material3minend = stat("mmat-interpolation-nopc.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat nopc_divergenceminstart = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] - + from fluidity_tools import stat_parser as stat nopc_divergenceminend = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] - + from fluidity_tools import stat_parser as stat nopc_divergencemaxstart = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] - + from fluidity_tools import stat_parser as stat nopc_divergencemaxend = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] - + from fluidity_tools import stat_parser as stat nopc_cflmaxstart = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] - + from fluidity_tools import stat_parser as stat nopc_cflmaxend = stat("mmat-interpolation-nopc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] - + from fluidity_tools import stat_parser as stat nopc_lambdamaxstart = stat("mmat-interpolation-nopc.stat")["Material1"]["Lambda"]["max"][0] - + from fluidity_tools import stat_parser as stat nopc_lambdamaxend = stat("mmat-interpolation-nopc.stat")["Material1"]["Lambda"]["max"][-1] - + from fluidity_tools import stat_parser as stat nopc_lambdaminstart = stat("mmat-interpolation-nopc.stat")["Material1"]["Lambda"]["min"][0] - + from fluidity_tools import stat_parser as stat nopc_lambdaminend = stat("mmat-interpolation-nopc.stat")["Material1"]["Lambda"]["min"][-1] - + from fluidity_tools import stat_parser as stat nopc_absolutedifference = stat("mmat-interpolation-nopc.stat")["Material3"]["ScalarAbsoluteDifference"]["max"][1] - + from fluidity_tools import stat_parser as stat masspc_material1integralstart = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat masspc_material1integralend = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat masspc_material1maxstart = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat masspc_material1maxend = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat masspc_material1minstart = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat masspc_material1minend = stat("mmat-interpolation-masspc.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat masspc_material2integralstart = stat("mmat-interpolation-masspc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat masspc_material2integralend = stat("mmat-interpolation-masspc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat masspc_material3integralstart = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat masspc_material3integralend = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat masspc_material3maxstart = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat masspc_material3maxend = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat masspc_material3minstart = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat masspc_material3minend = stat("mmat-interpolation-masspc.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat masspc_divergenceminstart = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] - + from fluidity_tools import stat_parser as stat masspc_divergenceminend = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] - + from fluidity_tools import stat_parser as stat masspc_divergencemaxstart = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] - + from fluidity_tools import stat_parser as stat masspc_divergencemaxend = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] - + from fluidity_tools import stat_parser as stat masspc_cflmaxstart = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] - + from fluidity_tools import stat_parser as stat masspc_cflmaxend = stat("mmat-interpolation-masspc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] - + from fluidity_tools import stat_parser as stat masspc_lambdamaxstart = stat("mmat-interpolation-masspc.stat")["Material1"]["Lambda"]["max"][0] - + from fluidity_tools import stat_parser as stat masspc_lambdamaxend = stat("mmat-interpolation-masspc.stat")["Material1"]["Lambda"]["max"][-1] - + from fluidity_tools import stat_parser as stat masspc_lambdaminstart = stat("mmat-interpolation-masspc.stat")["Material1"]["Lambda"]["min"][0] - + from fluidity_tools import stat_parser as stat masspc_lambdaminend = stat("mmat-interpolation-masspc.stat")["Material1"]["Lambda"]["min"][-1] - + from fluidity_tools import stat_parser as stat masspc_absolutedifference = stat("mmat-interpolation-masspc.stat")["Material3"]["ScalarAbsoluteDifference"]["max"][1] - + from fluidity_tools import stat_parser as stat lumppc_material1integralstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat lumppc_material1integralend = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat lumppc_material1maxstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat lumppc_material1maxend = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat lumppc_material1minstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat lumppc_material1minend = stat("mmat-interpolation-lumppc.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat lumppc_material2integralstart = stat("mmat-interpolation-lumppc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat lumppc_material2integralend = stat("mmat-interpolation-lumppc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat lumppc_material3integralstart = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat lumppc_material3integralend = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat lumppc_material3maxstart = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat lumppc_material3maxend = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat lumppc_material3minstart = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat lumppc_material3minend = stat("mmat-interpolation-lumppc.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat lumppc_divergenceminstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] - + from fluidity_tools import stat_parser as stat lumppc_divergenceminend = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] - + from fluidity_tools import stat_parser as stat lumppc_divergencemaxstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] - + from fluidity_tools import stat_parser as stat lumppc_divergencemaxend = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] - + from fluidity_tools import stat_parser as stat lumppc_cflmaxstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] - + from fluidity_tools import stat_parser as stat lumppc_cflmaxend = stat("mmat-interpolation-lumppc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] - + from fluidity_tools import stat_parser as stat lumppc_lambdamaxstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["Lambda"]["max"][0] - + from fluidity_tools import stat_parser as stat lumppc_lambdamaxend = stat("mmat-interpolation-lumppc.stat")["Material1"]["Lambda"]["max"][-1] - + from fluidity_tools import stat_parser as stat lumppc_lambdaminstart = stat("mmat-interpolation-lumppc.stat")["Material1"]["Lambda"]["min"][0] - + from fluidity_tools import stat_parser as stat lumppc_lambdaminend = stat("mmat-interpolation-lumppc.stat")["Material1"]["Lambda"]["min"][-1] - + from fluidity_tools import stat_parser as stat lumppc_absolutedifference = stat("mmat-interpolation-lumppc.stat")["Material3"]["ScalarAbsoluteDifference"]["max"][1] - + from fluidity_tools import stat_parser as stat diagpc_material1integralstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat diagpc_material1integralend = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat diagpc_material1maxstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat diagpc_material1maxend = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat diagpc_material1minstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat diagpc_material1minend = stat("mmat-interpolation-diagpc.stat")["Material1"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat diagpc_material2integralstart = stat("mmat-interpolation-diagpc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat diagpc_material2integralend = stat("mmat-interpolation-diagpc.stat")["Material2"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat diagpc_material3integralstart = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][0] - + from fluidity_tools import stat_parser as stat diagpc_material3integralend = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["integral"][-1] - + from fluidity_tools import stat_parser as stat diagpc_material3maxstart = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["max"][0] - + from fluidity_tools import stat_parser as stat diagpc_material3maxend = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["max"][-1] - + from fluidity_tools import stat_parser as stat diagpc_material3minstart = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["min"][0] - + from fluidity_tools import stat_parser as stat diagpc_material3minend = stat("mmat-interpolation-diagpc.stat")["Material3"]["MaterialVolumeFraction"]["min"][-1] - + from fluidity_tools import stat_parser as stat diagpc_divergenceminstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeDivergence"]["min"][0] - + from fluidity_tools import stat_parser as stat diagpc_divergenceminend = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeDivergence"]["min"][-1] - + from fluidity_tools import stat_parser as stat diagpc_divergencemaxstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeDivergence"]["max"][0] - + from fluidity_tools import stat_parser as stat diagpc_divergencemaxend = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeDivergence"]["max"][-1] - + from fluidity_tools import stat_parser as stat diagpc_cflmaxstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][0] - + from fluidity_tools import stat_parser as stat diagpc_cflmaxend = stat("mmat-interpolation-diagpc.stat")["Material1"]["ControlVolumeCFLNumber"]["max"][-1] - + from fluidity_tools import stat_parser as stat diagpc_lambdamaxstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["Lambda"]["max"][0] - + from fluidity_tools import stat_parser as stat diagpc_lambdamaxend = stat("mmat-interpolation-diagpc.stat")["Material1"]["Lambda"]["max"][-1] - + from fluidity_tools import stat_parser as stat diagpc_lambdaminstart = stat("mmat-interpolation-diagpc.stat")["Material1"]["Lambda"]["min"][0] - + from fluidity_tools import stat_parser as stat diagpc_lambdaminend = stat("mmat-interpolation-diagpc.stat")["Material1"]["Lambda"]["min"][-1] - + from fluidity_tools import stat_parser as stat diagpc_absolutedifference = stat("mmat-interpolation-diagpc.stat")["Material3"]["ScalarAbsoluteDifference"]["max"][1]