diff --git a/.github/workflows/ubuntu.yml b/.github/workflows/ubuntu.yml index 8aab030e36..41c500aa4d 100644 --- a/.github/workflows/ubuntu.yml +++ b/.github/workflows/ubuntu.yml @@ -4,6 +4,8 @@ on: push: branches: - main + - solenoidal_full_schur + - cianwilson/solenoidal_full_schur pull_request: branches: - main diff --git a/assemble/Full_Projection.F90 b/assemble/Full_Projection.F90 index 972386592b..02f767ea1d 100644 --- a/assemble/Full_Projection.F90 +++ b/assemble/Full_Projection.F90 @@ -36,6 +36,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 @@ -60,7 +61,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) !! @@ -74,11 +75,12 @@ 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 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 @@ -87,6 +89,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 @@ -99,11 +102,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' @@ -118,7 +135,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.) @@ -136,7 +153,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) !-------------------------------------------------------------------------------------------------------- @@ -170,11 +187,11 @@ 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 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 @@ -201,39 +218,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 - logical parallel, have_auxiliary_matrix, have_preconditioner_matrix + integer reference_node, i, rotation_stat, ref_stat + logical parallel, have_auxiliary_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)) @@ -277,17 +291,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=preconditioner_matrix%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 @@ -419,12 +434,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! - ! 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 @@ -451,7 +461,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/Makefile.dependencies b/assemble/Makefile.dependencies index 8982285d40..66c3c32391 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 \ @@ -730,11 +730,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_CG.F90 b/assemble/Momentum_CG.F90 index 4aeceb5863..e8ff294189 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/Momentum_Equation.F90 b/assemble/Momentum_Equation.F90 index 16f6fbc3a1..ad1662e378 100644 --- a/assemble/Momentum_Equation.F90 +++ b/assemble/Momentum_Equation.F90 @@ -491,7 +491,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 859e6607a8..c9026335ae 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") - lump_mass = have_option(trim(l_option_path)//& - "/interpolated_field/discontinuous/lump_mass_matrix") & - .or. .not.dg + if (discontinuous.and..not.dg) then + FLExit("Using discontinuous solenoidal projection options but field is not discontinuous.") + end if - lump_on_submesh = have_option(trim(l_option_path)//& + assemble_lumped_mass = have_option(trim(l_option_path)//& + "/interpolated_field/discontinuous/lump_mass_matrix") & + .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_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") @@ -182,9 +251,30 @@ subroutine solenoidal_interpolation_fields(v_field, coordinate, & call allocate(ct_rhs, lagrange_mesh, "DivergenceRHS") call zero(ct_rhs) - 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.or.assemble_schur_aux.or.assemble_cdiagmc) then + cmc_m_sparsity = make_sparsity_transpose(lagrange_mesh, v_field%mesh, "LagrangeProjectionSparsity") + end if + + 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.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 + + 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) + 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 + 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 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(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.") + 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,86 @@ 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, & + 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(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 b3640699e8..fbd20af72b 100644 --- a/femtools/FEFields.F90 +++ b/femtools/FEFields.F90 @@ -9,6 +9,7 @@ module fefields use elements, only: element_type use parallel_tools use sparse_tools + use sparse_tools_petsc use parallel_fields use transform_elements, only: transform_to_physical, element_volume use fetools, only: shape_shape, shape_rhs, shape_vector_rhs @@ -27,6 +28,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 @@ -213,7 +218,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 @@ -228,7 +233,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 @@ -260,7 +265,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 6c75b632a4..de5f28b6e8 100644 --- a/femtools/Makefile.dependencies +++ b/femtools/Makefile.dependencies @@ -367,7 +367,8 @@ FEFields.o ../include/fefields.mod: FEFields.F90 \ ../include/halos.mod ../include/mpi_interfaces.mod \ ../include/parallel_fields.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 cb3fd56872..6d13c4a757 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 583d05585f..181a78012a 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 c7eeefccd7..6cf669199f 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 9cda3550ee..4b220e3e1d 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 9d20da7812..c0f1d131e1 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 eed62ce573..435a1b18f3 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 30d0e24c9a..004c427b22 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 81af73b7eb..32931059d0 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 a6ce6a05ac..413b63cefa 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 = ( 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 3bbccf35c1..425b168fe9 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 a5d514b03e..9553c4cacb 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 e81caead39..becc4978d4 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 @@ -85,17 +116,17 @@ linear_solver_options_sym = 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 05048a6b12..aa2c516203 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 5f79541f92..d23fb5a9d9 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 fa751b4370..e90169fd28 100644 --- a/schemas/test_advection_diffusion_options.rng +++ b/schemas/test_advection_diffusion_options.rng @@ -828,7 +828,7 @@ alternating directions. Devised by C.Pain. Classical interior penalty scheme see, e.g., SIAM Journal on Numerical Analysis -Vol. 39, No. 5 (2002), pp. 1749-1779 +Vol. 39, No. 5 (2002), pp. 1749-1779 Penalty_parameter The penalty term Int [u][v] dS on element boundaries @@ -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 - + @@ -2277,7 +2198,7 @@ for ABL simulations use a log profile Reynolds stresses profile usually a function of height, -assumes that the remaining stresses are negligible +assumes that the remaining stresses are negligible @@ -2797,7 +2718,7 @@ enabled) and file (if /io/convergence_file is enabled) - Exclude this field from convergence testing and file + Exclude this field from convergence testing and file @@ -2832,7 +2753,7 @@ i.e. excluding the components - Exclude this field entirely from convergence testing and file + Exclude this field entirely from convergence testing and file @@ -2882,7 +2803,7 @@ i.e. excluding the components - Exclude this field entirely from convergence testing and file + Exclude this field entirely from convergence testing and file @@ -2946,7 +2867,7 @@ fields. - + @@ -2961,7 +2882,7 @@ fields. one specifies the absolute interpolation error in the units of the field that is being adapted, e.g. you can specify -the error to be 1.3 units +the error to be 1.3 units 0 @@ -3034,7 +2955,7 @@ Typically, tau will be ~= 6-8 * |e|_H1. one specifies the absolute interpolation error in the units of the field that is being adapted, e.g. you can specify -the error to be 1.3 units +the error to be 1.3 units 0 @@ -3115,7 +3036,7 @@ Typically, tau will be ~= 6-8 * |e|_H1. one specifies the absolute interpolation error in the units of the field that is being adapted, e.g. you can specify -the error to be 1.3 units +the error to be 1.3 units 1 @@ -3173,7 +3094,7 @@ Source: Fluidity/ICOM manual draft version 1.2 one specifies the absolute interpolation error in the units of the field that is being adapted, e.g. you can specify -the error to be 1.3 units +the error to be 1.3 units 1 @@ -3231,7 +3152,7 @@ Source: Fluidity/ICOM manual draft version 1.2 one specifies the absolute interpolation error in the units of the field that is being adapted, e.g. you can specify -the error to be 1.3 units +the error to be 1.3 units 2 @@ -3289,7 +3210,7 @@ Source: Fluidity/ICOM manual draft version 1.2 one specifies the absolute interpolation error in the units of the field that is being adapted, e.g. you can specify -the error to be 1.3 units +the error to be 1.3 units 2 @@ -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 - - - - - - Diagnostic scalar fields below this @@ -4322,24 +4226,6 @@ Adapting to this field is not recommended - - - 0 - - - SolidConcentration - - - - - - - - - - - - This scalar field is meant to replace DENTRAF. Basically, if you use new options, DENTRAF is no longer needed @@ -4356,82 +4242,6 @@ No repointing is done from this field to DENTRAF. - - Add field to be used by Solid_configuration to -Visualize the solids and MaterialVolumeFraction together - - 0 - - - VisualizeSolidFluid - - - - - - - - - Add field to be used by Solid_configuration to -Visualize the solid_Concentration - - 0 - - - VisualizeSolid - - - - - - - - - Add field to be used by Solid_configuration to -map the solid_Concentration from particle mesh to -the fluid mesh. - - 0 - - - ParticleScalar - - - - - - - - - Add field to be used by Explicit_ALE to -visualize functional values before iterations start. - - 0 - - - FunctionalBegin - - - - - - - - - Add field to be used by Explicit_ALE to -visualize functional values at each iteration. - - 0 - - - FunctionalIter - - - - - - - add a MaterialVolume scalar_field to calculate the spatially varying volume of a material (requires a MaterialVolumeFraction) @@ -4830,7 +4640,7 @@ strong Dirichlet boundary condition of 0 on all boundaries. Solver - + @@ -4863,30 +4673,6 @@ Limitations: - - Volume of the vehicles - -used in Traffic Modelling - - 0 - - - SolidPhase - - - IDENT = -42 - - - - - - - - - - - - Absolute Difference between two scalar fields. @@ -4990,13 +4776,13 @@ of this field is continuous. Lump the mass matrix of the galerkin projection -less accurate but faster and might give smoother result. +less accurate but faster and might give smoother result. - + @@ -5556,22 +5342,22 @@ ocean biology is being simulated. Prescribed vector fields below this Prescribed vector fields below this @@ -5879,97 +5665,6 @@ Limitations: - - Solid Velocity field. Used to generate the momentum source - - 1 - - - SolidVelocity - - - - - - - - - Same as Solid Velocity field but it is on the Particle mesh. -It is used to map the velocities coming from an external program like -FEMDEM or DEM to the fluid mesh. - - 1 - - - ParticleVector - - - - - - - - - Same as Solid Velocity field but it is on the Particle mesh. -It is used to map the velocities coming from an external program like -FEMDEM or DEM to the fluid mesh. - - 1 - - - ParticleForce - - - - - - - - - Same as Solid Velocity field but it is on the Particle mesh. -It is used to map the velocities coming from an external program like -FEMDEM or DEM to the fluid mesh. - - 1 - - - SolidForce - - - - - - - - - - 1 - - - VelocityPlotForSolids - - - - - - - - - Same as Solid Velocity field but it is on the Particle mesh. -It is used to map the velocities coming from an external program like -FEMDEM or DEM to the fluid mesh. - - 1 - - - FunctionalGradient - - - - - - - LinearMomentum field. p = \rho*u @@ -6121,13 +5816,13 @@ of this field is continuous. Lump the mass matrix of the galerkin projection -less accurate but faster and might give smoother result. +less accurate but faster and might give smoother result. - + @@ -6157,7 +5852,7 @@ is necessary. - + @@ -6244,11 +5939,11 @@ This diagnostic depends on ImbalancedVelocity. Prescribed scalar fields below this @@ -6294,7 +5989,7 @@ in a multimaterial simulation. Diagnostic tensor fields below this @@ -6651,7 +6346,7 @@ http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/P - + Iterative (Krylov) method to solve the linear discretised equation Given are the most frequently used methods. The solution is done @@ -6835,7 +6530,7 @@ http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/P @@ -7255,7 +6950,7 @@ higher order solution is pivotted. This enables DG-specific timestepping options, such as -explicit advection subcycling. +explicit advection subcycling. @@ -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 @@ -8032,7 +7622,7 @@ divided by the total volume of the domain (continuous galerkin) or within the control volume (control_volumes). - if select normalise the volume fractions will be divided by the total volume of the domain + if select normalise the volume fractions will be divided by the total volume of the domain @@ -8043,7 +7633,7 @@ divided by the total volume of the domain e.g. the values 0 1 2 3 will return 4 bins and the fraction of the field in each bin with, 0<=field<1, 1<=field<2, 2<=field<3, 3<=field, -will be calculated. +will be calculated. @@ -8061,7 +7651,7 @@ Defaults to zero at machine tolerance (epsilon(0.0)) if not selected. - + @@ -8266,13 +7856,6 @@ Whatever field is selected must be present. - - Option to solve for electrical potential from -electrokinetic, electrochemical or electrothermal sources - - ElectricalPotential - - @@ -8392,17 +7975,17 @@ If you enable this you MUST enable the /geometry/ocean_boundaries option too Model of biological processes in the ocean. - A simple model of phytoplankton, zooplankton, general nutrient and detritus. + A simple model of phytoplankton, zooplankton, general nutrient and detritus. Python code specifying the source and sink relationships between the biological tracers. This is usually achieved by -importing fluidity.ocean_biology and calling a scheme from there. +importing fluidity.ocean_biology and calling a scheme from there. Do not calculate sources and sinks. -This option is generally only useful for testing. +This option is generally only useful for testing. @@ -8502,7 +8085,7 @@ This option is generally only useful for testing. Discontinuous galerkin formulation. You can also use this formulation with a continuous field in which case a simple -galerkin formulation will result. +galerkin formulation will result. @@ -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. @@ -8889,7 +8472,7 @@ options too and it must be on the same mesh as used for the solenoidal projection above. Note well: this only really makes sense for scalar fields linked to nondivergent -vector fields, such as pressure to incompressible velocity! +vector fields, such as pressure to incompressible velocity! Pressure @@ -8903,7 +8486,7 @@ options too and it must be on the same mesh as used for the solenoidal projection above. Note well: this only really makes sense for scalar fields linked to nondivergent -vector fields, such as pressure to incompressible velocity! +vector fields, such as pressure to incompressible velocity! @@ -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. - + 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 + + + + 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 + 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-parallel-p2p1cv/mmat-interpolation-diagpc.flml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-diagpc.flml new file mode 100644 index 0000000000..7403317f95 --- /dev/null +++ b/tests/mmat-interpolation-parallel-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-parallel-p2p1cv/mmat-interpolation-lumppc.flml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-lumppc.flml new file mode 100644 index 0000000000..3053041be8 --- /dev/null +++ b/tests/mmat-interpolation-parallel-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-parallel-p2p1cv/mmat-interpolation-masspc.flml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-masspc.flml new file mode 100644 index 0000000000..f31f15ce41 --- /dev/null +++ b/tests/mmat-interpolation-parallel-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-parallel-p2p1cv/mmat-interpolation-nopc.flml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-nopc.flml new file mode 100644 index 0000000000..d4361c4e96 --- /dev/null +++ b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-nopc.flml @@ -0,0 +1,620 @@ + + + + mmat-interpolation-nopc + + + 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-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml new file mode 100644 index 0000000000..fc39e2a224 --- /dev/null +++ b/tests/mmat-interpolation-parallel-p2p1cv/mmat-interpolation-parallel-p2p1cv.xml @@ -0,0 +1,675 @@ + + + Multimaterial adaptivity and interpolation test + + flml 2dadapt + + 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 + + + + +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] + + + + + assert(solvers_converged) + + +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(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(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(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(nopc_material1maxstart-nopc_material1maxend) < 1.E-10 + + +assert abs(nopc_material3maxstart-nopc_material3maxend) < 1.E-10 + + +assert abs(nopc_material1minstart-nopc_material1minend) < 1.E-10 + + +assert abs(nopc_material3minstart-nopc_material3minend) < 1.E-10 + + +assert abs(nopc_divergencemaxstart) < 1.E-4 + + +assert abs(nopc_divergenceminstart) < 1.E-4 + + +assert abs(nopc_divergencemaxend) < 1.E-4 + + +assert abs(nopc_divergenceminend) < 1.E-4 + + +assert abs(nopc_cflmaxstart-0.1) < 1.E-10 + + +assert abs(nopc_cflmaxend-0.1) < 1.E-10 + + +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(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(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(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(masspc_material1maxstart-masspc_material1maxend) < 1.E-10 + + +assert abs(masspc_material3maxstart-masspc_material3maxend) < 1.E-10 + + +assert abs(masspc_material1minstart-masspc_material1minend) < 1.E-10 + + +assert abs(masspc_material3minstart-masspc_material3minend) < 1.E-10 + + +assert abs(masspc_divergencemaxstart) < 1.E-4 + + +assert abs(masspc_divergenceminstart) < 1.E-4 + + +assert abs(masspc_divergencemaxend) < 1.E-4 + + +assert abs(masspc_divergenceminend) < 1.E-4 + + +assert abs(masspc_cflmaxstart-0.1) < 1.E-10 + + +assert abs(masspc_cflmaxend-0.1) < 1.E-10 + + +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(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(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(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(lumppc_lambdamaxend) > 1.E-7 + + +assert abs(lumppc_lambdaminend) > 1.E-7 + + +assert abs(lumppc_absolutedifference) > 1.E-7 + + +assert abs(diagpc_lambdamaxstart) > 1.E-7 + + +assert abs(diagpc_lambdaminstart) > 1.E-7 + + +assert abs(diagpc_lambdamaxend) > 1.E-7 + + +assert abs(diagpc_lambdaminend) > 1.E-7 + + +assert abs(diagpc_absolutedifference) > 1.E-7 + + + diff --git a/tests/mmat-interpolation-parallel-p2p1cv/src/2d_square.geo b/tests/mmat-interpolation-parallel-p2p1cv/src/2d_square.geo new file mode 100644 index 0000000000..70cb64b6af --- /dev/null +++ b/tests/mmat-interpolation-parallel-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};