diff --git a/assemble/Diagnostic_fields_wrapper.F90 b/assemble/Diagnostic_fields_wrapper.F90 index a3f2a6c5dc..a66137512a 100644 --- a/assemble/Diagnostic_fields_wrapper.F90 +++ b/assemble/Diagnostic_fields_wrapper.F90 @@ -355,6 +355,13 @@ subroutine calculate_diagnostic_variables(state, exclude_nonrecalculated) end if end if + v_field => extract_vector_field(state(i), "InitialCoordinate", stat) + if(stat == 0) then + if(recalculate(trim(v_field%option_path))) then + call calculate_diagnostic_variable(state(i), "InitialCoordinate", v_field) + end if + end if + s_field => extract_scalar_field(state(i), "UniversalNumber", stat) if(stat == 0) then call calculate_diagnostic_variable(state(i), "UniversalNumber", s_field) diff --git a/assemble/Makefile.dependencies b/assemble/Makefile.dependencies index 5809574f9e..477b54d21c 100644 --- a/assemble/Makefile.dependencies +++ b/assemble/Makefile.dependencies @@ -23,9 +23,10 @@ Adapt_State.o ../include/adapt_state_module.mod: Adapt_State.F90 \ ../include/detector_parallel.mod ../include/diagnostic_fields_wrapper.mod \ ../include/diagnostic_fields_wrapper_new.mod \ ../include/diagnostic_variables.mod \ - ../include/discrete_properties_module.mod ../include/edge_length_module.mod \ - ../include/elements.mod ../include/eventcounter.mod ../include/fdebug.h \ - ../include/fefields.mod ../include/field_options.mod ../include/fields.mod \ + ../include/discrete_properties_module.mod ../include/dqmom.mod \ + ../include/edge_length_module.mod ../include/elements.mod \ + ../include/eventcounter.mod ../include/fdebug.h ../include/fefields.mod \ + ../include/field_options.mod ../include/fields.mod \ ../include/fields_halos.mod ../include/fldebug.mod ../include/futils.mod \ ../include/global_parameters.mod ../include/hadapt_extrude.mod \ ../include/hadapt_metric_based_extrude.mod ../include/halos.mod \ @@ -214,9 +215,10 @@ Diagnostic_Fields_Matrices.o ../include/diagnostic_fields_matrices.mod: \ Diagnostic_fields_wrapper.o ../include/diagnostic_fields_wrapper.mod: \ Diagnostic_fields_wrapper.F90 ../include/diagnostic_fields.mod \ - ../include/diagnostic_fields_matrices.mod ../include/equation_of_state.mod \ - ../include/fdebug.h ../include/fetools.mod ../include/field_derivatives.mod \ - ../include/field_options.mod ../include/fields.mod ../include/fldebug.mod \ + ../include/diagnostic_fields_matrices.mod ../include/dqmom.mod \ + ../include/equation_of_state.mod ../include/fdebug.h ../include/fetools.mod \ + ../include/field_derivatives.mod ../include/field_options.mod \ + ../include/fields.mod ../include/fldebug.mod \ ../include/free_surface_module.mod ../include/futils.mod \ ../include/geostrophic_pressure.mod ../include/global_parameters.mod \ ../include/momentum_diagnostic_fields.mod \ @@ -526,13 +528,16 @@ Mba3d_Integration.o ../include/mba3d_integration.mod: Mba3d_Integration.F90 \ @true MeshMovement.o ../include/meshmovement.mod: MeshMovement.F90 \ - ../include/element_numbering.mod ../include/elements.mod \ - ../include/eventcounter.mod ../include/fdebug.h ../include/fefields.mod \ - ../include/fetools.mod ../include/field_derivatives.mod \ - ../include/fields.mod ../include/fields_base.mod ../include/fldebug.mod \ + ../include/boundary_conditions.mod ../include/element_numbering.mod \ + ../include/elements.mod ../include/eventcounter.mod ../include/fdebug.h \ + ../include/fefields.mod ../include/fetools.mod \ + ../include/field_derivatives.mod ../include/fields.mod \ + ../include/fields_base.mod ../include/fldebug.mod \ ../include/global_numbering.mod ../include/global_parameters.mod \ - ../include/shape_functions.mod ../include/solvers.mod \ - ../include/sparse_matrices_fields.mod ../include/sparse_tools.mod \ + ../include/parallel_fields.mod ../include/petsc_legacy.h \ + ../include/profiler.mod ../include/shape_functions.mod \ + ../include/solvers.mod ../include/sparse_matrices_fields.mod \ + ../include/sparse_tools.mod ../include/sparse_tools_petsc.mod \ ../include/sparsity_patterns.mod ../include/sparsity_patterns_meshes.mod \ ../include/state_module.mod ../include/unittest_tools.mod \ ../include/vector_tools.mod ../include/vtk_interfaces.mod diff --git a/assemble/Makefile.in b/assemble/Makefile.in index 6046b232eb..23a3d364f4 100644 --- a/assemble/Makefile.in +++ b/assemble/Makefile.in @@ -75,8 +75,8 @@ OBJS = Momentum_Diagnostic_Fields.o \ Petsc_Solve_State.o Adaptivity_1D.o \ Zoltan_global_variables.o Zoltan_integration.o Zoltan_callbacks.o Zoltan_detectors.o \ Diagnostic_Children.o Reduced_Model_Runtime.o Turbine.o Implicit_Solids.o \ - Manifold_Projections.o Hybridized_Helmholtz.o Burgers_Assembly.o Pressure_Dirichlet_BCS_CV.o \ - Porous_Media.o + Manifold_Projections.o Hybridized_Helmholtz.o Burgers_Assembly.o Pressure_Dirichlet_BCS_CV.o \ + Porous_Media.o laplacian_smoother.o lineal_smoother.o lineal_torsional_smoother.o centroid_relaxer.o .SUFFIXES: .F90 .c .cpp .o .a diff --git a/assemble/MeshMovement.F90 b/assemble/MeshMovement.F90 index c6d4237631..ef8d70ff90 100644 --- a/assemble/MeshMovement.F90 +++ b/assemble/MeshMovement.F90 @@ -2,6 +2,7 @@ module meshmovement + use iso_c_binding use global_parameters use fldebug use vector_tools, only: solve, mat_diag_mat, eigendecomposition_symmetric @@ -15,14 +16,19 @@ module meshmovement use eventcounter use fetools use unittest_tools + use parallel_fields use fields + use profiler use state_module + use boundary_conditions, only : get_boundary_condition, get_boundary_condition_count,& + extract_surface_field use vtk_interfaces use sparse_matrices_fields use solvers use fefields use field_derivatives use sparsity_patterns + use sparse_tools_petsc use sparsity_patterns_meshes implicit none @@ -38,10 +44,116 @@ end subroutine set_debug_level subroutine reset_debug_level end subroutine reset_debug_level end interface - + + interface + subroutine lap_smoother(dim, num_nodes,& + num_elements, num_surf_elements, num_owned_nodes,& + mapping, connectivity, phys_mesh, smooth_mesh,& + comp_mesh, surf_connectivity, matrix,& + options, debug_level) bind(c) + use iso_c_binding +#ifdef HAVE_PETSC_MODULES + use petsc +#endif + use solvers, only: solver_options_type + implicit none +#include "petsc_legacy.h" + + integer (c_int), value :: dim, num_nodes, num_elements, & + num_surf_elements, num_owned_nodes, debug_level + integer (c_int), dimension(num_nodes) :: mapping + integer (c_int), dimension((dim+1)*num_elements) :: connectivity + real (c_double), dimension(num_nodes) :: phys_mesh, smooth_mesh, comp_mesh + integer (c_int), dimension(dim*num_elements) :: surf_connectivity + Mat :: matrix + type(solver_options_type) :: options + + end subroutine lap_smoother + end interface + + interface + subroutine lin_smoother(dim, num_nodes,& + num_elements, num_surf_elements, num_owned_nodes,& + mapping, connectivity, phys_mesh, smooth_mesh,& + comp_mesh, surf_connectivity,findrm, colm, matrix,& + options, debug_level) bind(c) + use iso_c_binding +#ifdef HAVE_PETSC_MODULES + use petsc +#endif + use solvers, only: solver_options_type + implicit none +#include "petsc_legacy.h" + + integer (c_int), value :: dim, num_nodes, num_elements, & + num_surf_elements, num_owned_nodes, debug_level + integer (c_int), dimension(num_nodes,dim) :: mapping + integer (c_int), dimension((dim+1)*num_elements) :: connectivity + real (c_double), dimension(num_nodes) :: phys_mesh, smooth_mesh, comp_mesh + integer (c_int), dimension(dim*num_elements) :: surf_connectivity + integer (c_int), dimension(num_nodes+1) :: findrm + integer (c_int), dimension(findrm(num_nodes+1)-1) :: colm + Mat :: matrix + type(solver_options_type) :: options + + end subroutine lin_smoother + end interface + + interface + subroutine lin_tor_smoother(dim, num_nodes,& + num_elements, num_surf_elements, num_owned_nodes,& + mapping, connectivity, phys_mesh, smooth_mesh,& + comp_mesh, surf_connectivity,findrm, colm, matrix,& + options, debug_level) bind(c) + use iso_c_binding +#ifdef HAVE_PETSC_MODULES + use petsc +#endif + use solvers, only: solver_options_type + implicit none +#include "petsc_legacy.h" + + integer (c_int), value :: dim, num_nodes, num_elements, & + num_surf_elements, num_owned_nodes, debug_level + integer (c_int), dimension(num_nodes,dim) :: mapping + integer (c_int), dimension((dim+1)*num_elements) :: connectivity + real (c_double), dimension(num_nodes) :: phys_mesh, smooth_mesh, comp_mesh + integer (c_int), dimension(dim*num_elements) :: surf_connectivity + integer (c_int), dimension(num_nodes+1) :: findrm + integer (c_int), dimension(findrm(num_nodes+1)-1) :: colm + Mat :: matrix + type(solver_options_type) :: options + end subroutine lin_tor_smoother + end interface + + interface + subroutine cent_relaxer(dim, num_nodes,& + num_elements, num_surf_elements, num_owned_nodes,& + mapping, connectivity, phys_mesh, smooth_mesh,& + comp_mesh, surf_connectivity,findrm, colm) bind(c) + use iso_c_binding + + implicit none + + integer (c_int), value :: dim, num_nodes, num_elements, & + num_surf_elements, num_owned_nodes + integer (c_int), dimension(num_nodes) :: mapping + integer (c_int), dimension((dim+1)*num_elements) :: connectivity + real (c_double), dimension(num_nodes) :: phys_mesh, smooth_mesh, comp_mesh + integer (c_int), dimension(dim*num_elements) :: surf_connectivity + integer (c_int), dimension(num_nodes+1) :: findrm + integer (c_int), dimension(findrm(num_nodes+1)-1) :: colm + + + end subroutine cent_relaxer + end interface private - public :: move_mesh_imposed_velocity, move_mesh_pseudo_lagrangian, movemeshy + public :: move_mesh_imposed_velocity, move_mesh_pseudo_lagrangian, movemeshy,& + move_mesh_laplacian_smoothing, move_mesh_initialise_laplacian_smoothing,& + move_mesh_lineal_smoothing, move_mesh_initialise_lineal_smoothing,& + move_mesh_lineal_torsional_smoothing, move_mesh_initialise_lineal_torsional_smoothing,& + move_mesh_centroid_relaxer,move_mesh_initialise_centroid_relaxer contains subroutine movemeshy(state,move_option,TimeStep) @@ -1174,7 +1286,7 @@ subroutine apply_lagrangian_multiplier(matrix,cdim,blocking) end subroutine apply_lagrangian_multiplier subroutine move_mesh_imposed_velocity(states) - type(state_type), dimension(:), intent(inout) :: states + type(state_type), dimension(:), intent(inout) :: states type(vector_field), pointer :: coordinate, old_coordinate, new_coordinate type(vector_field), pointer :: velocity @@ -1196,26 +1308,27 @@ subroutine move_mesh_imposed_velocity(states) new_coordinate => extract_vector_field(states(1), "IteratedCoordinate") call get_option("/timestepping/timestep", dt) + found_velocity = .false. do i = 1, size(states) - velocity => extract_vector_field(states(i), "Velocity", stat) - if(stat==0 .and. .not. velocity%aliased) then - call get_option(trim(velocity%option_path)//"/prognostic/temporal_discretisation/relaxation", itheta, stat) - if(found_velocity.and.(stat==0)) then - FLExit("Only one prognostic velocity allowed with imposed mesh movement.") - else - found_velocity = (stat==0) - end if - end if + velocity => extract_vector_field(states(i), "Velocity", stat) + if(stat==0 .and. .not. velocity%aliased) then + call get_option(trim(velocity%option_path)//"/prognostic/temporal_discretisation/relaxation", itheta, stat) + if(found_velocity.and.(stat==0)) then + FLExit("Only one prognostic velocity allowed with imposed mesh movement.") + else + found_velocity = (stat==0) + end if + end if end do if(.not.found_velocity) then - itheta = 0.5 + itheta = 0.5 end if call set(new_coordinate, old_coordinate) call addto(new_coordinate, grid_velocity, scale=dt) - + call set(coordinate, new_coordinate, old_coordinate, itheta) end subroutine move_mesh_imposed_velocity @@ -1280,6 +1393,563 @@ subroutine move_mesh_pseudo_lagrangian(states) end subroutine move_mesh_pseudo_lagrangian + + subroutine move_mesh_laplacian_smoothing(states, diagnostic_only) + type(state_type), dimension(:), intent(inout) :: states + logical, optional, intent(in) :: diagnostic_only + + type(vector_field), pointer :: coordinate, old_coordinate, new_coordinate,& + initial_coordinate + type(vector_field), pointer :: velocity + type(vector_field), pointer :: grid_velocity + type(mesh_type), pointer :: surface_mesh + + integer :: i, stat + real :: itheta, dt + logical :: found_velocity + type(solver_options_type) :: solver_options + type(petsc_csr_matrix) :: A + + !!! This routine drives a call to C code which does the actual assembly and + !!! solution for a Laplacian smoothed grid velocity. + + !!! The boundary conditions from the grid velocity field specified in diamond + !!! should be satisfied by this solution. + + if (.not. have_option("/mesh_adaptivity/mesh_movement/laplacian_smoothing")) return + + ewrite(1,*) 'Entering move_mesh_laplacian_smoothing' + + grid_velocity => extract_vector_field(states(1), "GridVelocity") + if (have_option('/mesh_adaptivity/mesh_movement/solver')) then + call populate_solver_options_struct_from_path(solver_options,& + '/mesh_adaptivity/mesh_movement/solver') + solver_options%ksptype=trim(solver_options%ksptype)//c_null_char + solver_options%pctype=trim(solver_options%pctype)//c_null_char + solver_options%hypretype=trim(solver_options%hypretype)//c_null_char + else + solver_options%ksptype=KSPGMRES//c_null_char + solver_options%pctype=PCSOR//c_null_char + solver_options%rtol=1.0e-7 + solver_options%atol=1.0e-7 + solver_options%max_its=1000 + end if + + call set_boundary_values(grid_velocity) + + coordinate => extract_vector_field(states(1), "Coordinate") + old_coordinate => extract_vector_field(states(1), "OldCoordinate") + new_coordinate => extract_vector_field(states(1), "IteratedCoordinate") + initial_coordinate => extract_vector_field(states(1), "InitialCoordinate") + surface_mesh => extract_mesh(states(1),"FullCoordinateSurfaceMesh") + + call profiler_tic("laplacian_smoother") + + call get_option("/timestepping/timestep", dt) + + + found_velocity = .false. + do i = 1, size(states) + velocity => extract_vector_field(states(i), "Velocity", stat) + if(stat==0 .and. .not. velocity%aliased) then + call get_option(trim(velocity%option_path)//"/prognostic/temporal_discretisation/relaxation", itheta, stat) + if(found_velocity.and.(stat==0)) then + FLExit("Only one prognostic velocity allowed with imposed mesh movement.") + else + found_velocity = (stat==0) + end if + end if + end do + if(.not.found_velocity) then + itheta = 0.5 + end if + + call set(coordinate, old_coordinate) + call addto(coordinate, grid_velocity, scale = dt) + call allocate(A, & + get_csr_sparsity_firstorder(states, coordinate%mesh, coordinate%mesh),& + [1,1],"StiffnessMatrix") + call zero(A) + + !!! actual call out to the C code. + +!Figure out indexing between C and Fortran + call lap_smoother(mesh_dim(coordinate), node_count(coordinate),& + element_count(coordinate), surface_element_count(coordinate),& + nowned_nodes(coordinate), universal_numbering(coordinate)-1,& + coordinate%mesh%ndglno, coordinate%val, new_coordinate%val,& + initial_coordinate%val, surface_mesh%ndglno,& + A%M, solver_options, debug_level()) + + call deallocate(A) + + + + call halo_update(new_coordinate) + + !!! Convert the new coordinates returned into a grid velocity + !!! and calculate the coordinate field at the theta time level. + + call set(grid_velocity,new_coordinate) + call addto(grid_velocity, old_coordinate, scale = -1.0) + call scale(grid_velocity, 1.0/dt) + + if (present_and_true(diagnostic_only)) then + call set(coordinate, old_coordinate) + call set(new_coordinate, old_coordinate) + else + call set(coordinate, new_coordinate, old_coordinate, itheta) + end if + + call profiler_toc("laplacian_smoother") + + end subroutine move_mesh_laplacian_smoothing + + subroutine move_mesh_initialise_laplacian_smoothing(states) + type(state_type), dimension(:), intent(inout) :: states + + type(vector_field), pointer :: coordinate + type(vector_field) :: initial_coordinate + type(mesh_type) :: surface_mesh + + !!! Store the initial mesh for use as the computational geometry for the + !! Laplacian mesh smoothing library. + + coordinate => extract_vector_field(states(1), "Coordinate") + + call allocate(initial_coordinate, mesh=coordinate%mesh, dim=coordinate%dim,& + name="InitialCoordinate") + call set(initial_coordinate, coordinate) + call insert(states(1), initial_coordinate, "InitialCoordinate") + call deallocate(initial_coordinate) + + !!! Add the surface mesh. This needs revisiting for parallel + + call extract_surface_mesh(surface_mesh, coordinate%mesh, "FullCoordinateSurfaceMesh") + call insert(states(1), surface_mesh, "FullCoordinateSurfaceMesh") + call deallocate(surface_mesh) + + end subroutine move_mesh_initialise_laplacian_smoothing + + subroutine move_mesh_lineal_smoothing(states, diagnostic_only) + type(state_type), dimension(:), intent(inout) :: states + logical, optional, intent(in) :: diagnostic_only + + type(vector_field), pointer :: coordinate, old_coordinate, new_coordinate,& + initial_coordinate + type(vector_field), pointer :: velocity + type(vector_field), pointer :: grid_velocity + type(mesh_type), pointer :: surface_mesh + + integer :: i, stat + real :: itheta, dt + logical :: found_velocity + type(solver_options_type) :: solver_options + type(petsc_csr_matrix) :: A + + !!! This routine drives a call to C code which does the actual assembly and + !!! solution for a Lineal Spring smoothed grid velocity. + + !!! The boundary conditions from the grid velocity field specified in diamond + !!! should be satisfied by this solution. + + if (.not. have_option("/mesh_adaptivity/mesh_movement/lineal_smoothing")) return + + ewrite(1,*) 'Entering move_mesh_lineal_smoothing' + + grid_velocity => extract_vector_field(states(1), "GridVelocity") + if (have_option('/mesh_adaptivity/mesh_movement/solver')) then + call populate_solver_options_struct_from_path(solver_options,& + '/mesh_adaptivity/mesh_movement/solver') + solver_options%ksptype=trim(solver_options%ksptype)//c_null_char + solver_options%pctype=trim(solver_options%pctype)//c_null_char + solver_options%hypretype=trim(solver_options%hypretype)//c_null_char + else + solver_options%ksptype=KSPGMRES//c_null_char + solver_options%pctype=PCSOR//c_null_char + solver_options%rtol=1.0e-7 + solver_options%atol=1.0e-7 + solver_options%max_its=1000 + end if + + call set_boundary_values(grid_velocity) + + coordinate => extract_vector_field(states(1), "Coordinate") + old_coordinate => extract_vector_field(states(1), "OldCoordinate") + new_coordinate => extract_vector_field(states(1), "IteratedCoordinate") + initial_coordinate => extract_vector_field(states(1), "InitialCoordinate") + surface_mesh => extract_mesh(states(1),"FullCoordinateSurfaceMesh") + + call profiler_tic(coordinate, "lineal_smoother") + + call get_option("/timestepping/timestep", dt) + + + found_velocity = .false. + do i = 1, size(states) + velocity => extract_vector_field(states(i), "Velocity", stat) + if(stat==0 .and. .not. velocity%aliased) then + call get_option(trim(velocity%option_path)//"/prognostic/temporal_discretisation/relaxation", itheta, stat) + if(found_velocity.and.(stat==0)) then + FLExit("Only one prognostic velocity allowed with imposed mesh movement.") + else + found_velocity = (stat==0) + end if + end if + end do + if(.not.found_velocity) then + itheta = 0.5 + end if + + call set(coordinate, old_coordinate) + call addto(coordinate, grid_velocity, scale = dt) + call allocate(A, & + get_csr_sparsity_firstorder(states, coordinate%mesh, coordinate%mesh),& + [mesh_dim(coordinate),mesh_dim(coordinate)],"StiffnessMatrix") + call zero(A) + + !!! actual call out to the C code. + + + if (.not. associated(coordinate%mesh%adj_lists%nnlist)) & + call add_nnlist(coordinate%mesh) + +!Figure out indexing between C and Fortran + call lin_smoother(mesh_dim(coordinate), node_count(coordinate),& + element_count(coordinate), surface_element_count(coordinate),& + nowned_nodes(coordinate), A%column_numbering%gnn2unn,& + coordinate%mesh%ndglno, coordinate%val, new_coordinate%val,& + initial_coordinate%val, surface_mesh%ndglno,& + coordinate%mesh%adj_lists%nnlist%findrm,& + coordinate%mesh%adj_lists%nnlist%colm, A%M,& + solver_options, debug_level()) + + + call halo_update(new_coordinate) + + !!! Convert the new coordinates returned into a grid velocity + !!! and calculate the coordinate field at the theta time level. + + call set(grid_velocity,new_coordinate) + call addto(grid_velocity, old_coordinate, scale = -1.0) + call scale(grid_velocity, 1.0/dt) + + if (present_and_true(diagnostic_only)) then + call set(coordinate, old_coordinate) + call set(new_coordinate, old_coordinate) + else + call set(coordinate, new_coordinate, old_coordinate, itheta) + end if + call deallocate(A) + + call profiler_toc(coordinate, "lineal_smoother") + + end subroutine move_mesh_lineal_smoothing + + subroutine move_mesh_initialise_lineal_smoothing(states) + type(state_type), dimension(:), intent(inout) :: states + + type(vector_field), pointer :: coordinate + type(vector_field) :: initial_coordinate + type(mesh_type) :: surface_mesh + + !!! Store the initial mesh for use as the computational geometry for the + !! Lineal mesh smoothing library. + + coordinate => extract_vector_field(states(1), "Coordinate") + + call allocate(initial_coordinate, mesh=coordinate%mesh, dim=coordinate%dim,& + name="InitialCoordinate") + call set(initial_coordinate, coordinate) + call insert(states(1), initial_coordinate, "InitialCoordinate") + call deallocate(initial_coordinate) + + !!! Add the surface mesh. This needs revisiting for parallel + + call extract_surface_mesh(surface_mesh, coordinate%mesh, "FullCoordinateSurfaceMesh") + call insert(states(1), surface_mesh, "FullCoordinateSurfaceMesh") + call deallocate(surface_mesh) + + end subroutine move_mesh_initialise_lineal_smoothing + + subroutine move_mesh_lineal_torsional_smoothing(states, diagnostic_only) + type(state_type), dimension(:), intent(inout) :: states + logical, optional, intent(in) :: diagnostic_only + + type(vector_field), pointer :: coordinate, old_coordinate, new_coordinate,& + initial_coordinate + type(vector_field), pointer :: velocity + type(vector_field), pointer :: grid_velocity + type(mesh_type), pointer :: surface_mesh + + integer :: i, stat + real :: itheta, dt + logical :: found_velocity + type(solver_options_type) :: solver_options + type(petsc_csr_matrix) :: A + + !!! This routine drives a call to C code which does the actual assembly and + !!! solution for a Lineal Torsional Spring smoothed grid velocity. + + !!! The boundary conditions from the grid velocity field specified in diamond + !!! should be satisfied by this solution. + + if (.not. have_option("/mesh_adaptivity/mesh_movement/lineal_torsional_smoothing")) return + + ewrite(1,*) 'Entering move_mesh_lineal_torsional_smoothing' + + grid_velocity => extract_vector_field(states(1), "GridVelocity") + if (have_option('/mesh_adaptivity/mesh_movement/solver')) then + call populate_solver_options_struct_from_path(solver_options,& + '/mesh_adaptivity/mesh_movement/solver') + solver_options%ksptype=trim(solver_options%ksptype)//c_null_char + solver_options%pctype=trim(solver_options%pctype)//c_null_char + solver_options%hypretype=trim(solver_options%hypretype)//c_null_char + else + solver_options%ksptype=KSPGMRES//c_null_char + solver_options%pctype=PCSOR//c_null_char + solver_options%rtol=1.0e-7 + solver_options%atol=1.0e-7 + solver_options%max_its=1000 + end if + + call set_boundary_values(grid_velocity) + + coordinate => extract_vector_field(states(1), "Coordinate") + old_coordinate => extract_vector_field(states(1), "OldCoordinate") + new_coordinate => extract_vector_field(states(1), "IteratedCoordinate") + initial_coordinate => extract_vector_field(states(1), "InitialCoordinate") + surface_mesh => extract_mesh(states(1),"FullCoordinateSurfaceMesh") + + call profiler_tic(coordinate, "lineal_torsional_smoother") + + call get_option("/timestepping/timestep", dt) + + + found_velocity = .false. + do i = 1, size(states) + velocity => extract_vector_field(states(i), "Velocity", stat) + if(stat==0 .and. .not. velocity%aliased) then + call get_option(trim(velocity%option_path)//"/prognostic/temporal_discretisation/relaxation", itheta, stat) + if(found_velocity.and.(stat==0)) then + FLExit("Only one prognostic velocity allowed with imposed mesh movement.") + else + found_velocity = (stat==0) + end if + end if + end do + if(.not.found_velocity) then + itheta = 0.5 + end if + + call set(coordinate, old_coordinate) + call addto(coordinate, grid_velocity, scale = dt) + call allocate(A, & + get_csr_sparsity_firstorder(states, coordinate%mesh, coordinate%mesh),& + [mesh_dim(coordinate),mesh_dim(coordinate)],"StiffnessMatrix") + call zero(A) + + !!! actual call out to the C code. + if (.not. associated(coordinate%mesh%adj_lists%nnlist)) & + call add_nnlist(coordinate%mesh) + +!Figure out indexing between C and Fortran + call lin_tor_smoother(mesh_dim(coordinate), node_count(coordinate),& + element_count(coordinate), surface_element_count(coordinate),& + nowned_nodes(coordinate), A%column_numbering%gnn2unn,& + coordinate%mesh%ndglno, coordinate%val, new_coordinate%val,& + initial_coordinate%val, surface_mesh%ndglno,& + coordinate%mesh%adj_lists%nnlist%findrm,& + coordinate%mesh%adj_lists%nnlist%colm, A%M,& + solver_options, debug_level()) + + call halo_update(new_coordinate) + + !!! Convert the new coordinates returned into a grid velocity + !!! and calculate the coordinate field at the theta time level. + + call set(grid_velocity,new_coordinate) + call addto(grid_velocity, old_coordinate, scale = -1.0) + call scale(grid_velocity, 1.0/dt) + + if (present_and_true(diagnostic_only)) then + call set(coordinate, old_coordinate) + call set(new_coordinate, old_coordinate) + else + call set(coordinate, new_coordinate, old_coordinate, itheta) + end if + call deallocate(A) + + call profiler_toc(coordinate, "lineal_torsional_smoother") + + end subroutine move_mesh_lineal_torsional_smoothing + + + subroutine move_mesh_initialise_lineal_torsional_smoothing(states) + type(state_type), dimension(:), intent(inout) :: states + + type(vector_field), pointer :: coordinate + type(vector_field) :: initial_coordinate + type(mesh_type) :: surface_mesh + + !!! Store the initial mesh for use as the computational geometry for the + !! centroid relaxation library. + + coordinate => extract_vector_field(states(1), "Coordinate") + + call allocate(initial_coordinate, mesh=coordinate%mesh, dim=coordinate%dim,& + name="InitialCoordinate") + call set(initial_coordinate, coordinate) + call insert(states(1), initial_coordinate, "InitialCoordinate") + call deallocate(initial_coordinate) + + !!! Add the surface mesh. This needs revisiting for parallel + + call extract_surface_mesh(surface_mesh, coordinate%mesh, "FullCoordinateSurfaceMesh") + call insert(states(1), surface_mesh, "FullCoordinateSurfaceMesh") + call deallocate(surface_mesh) + + end subroutine move_mesh_initialise_lineal_torsional_smoothing + + subroutine move_mesh_centroid_relaxer(states, diagnostic_only) + type(state_type), dimension(:), intent(inout) :: states + logical, optional, intent(in) :: diagnostic_only + + type(vector_field), pointer :: coordinate, old_coordinate, new_coordinate,& + initial_coordinate + type(vector_field), pointer :: velocity + type(vector_field), pointer :: grid_velocity + type(mesh_type), pointer :: surface_mesh + + integer :: i, stat + real :: itheta, dt + logical :: found_velocity + + !!! This routine drives a call to C code which does the actual assembly and + !!! solution for a Lineal Spring smoothed grid velocity. + + !!! The boundary conditions from the grid velocity field specified in diamond + !!! should be satisfied by this solution. + + if (.not. have_option("/mesh_adaptivity/mesh_movement/centroid_relaxer")) return + + ewrite(1,*) 'Entering move_mesh_centroid_relaxer' + + grid_velocity => extract_vector_field(states(1), "GridVelocity") + + call set_boundary_values(grid_velocity) + + coordinate => extract_vector_field(states(1), "Coordinate") + old_coordinate => extract_vector_field(states(1), "OldCoordinate") + new_coordinate => extract_vector_field(states(1), "IteratedCoordinate") + initial_coordinate => extract_vector_field(states(1), "InitialCoordinate") + surface_mesh => extract_mesh(states(1),"FullCoordinateSurfaceMesh") + + call get_option("/timestepping/timestep", dt) + + + found_velocity = .false. + do i = 1, size(states) + velocity => extract_vector_field(states(i), "Velocity", stat) + if(stat==0 .and. .not. velocity%aliased) then + call get_option(trim(velocity%option_path)//"/prognostic/temporal_discretisation/relaxation", itheta, stat) + if(found_velocity.and.(stat==0)) then + FLExit("Only one prognostic velocity allowed with imposed mesh movement.") + else + found_velocity = (stat==0) + end if + end if + end do + if(.not.found_velocity) then + itheta = 0.5 + end if + + call set(coordinate, old_coordinate) + call addto(coordinate, grid_velocity, scale = dt) + + !!! actual call out to the C code. + + + if (.not. associated(coordinate%mesh%adj_lists%nnlist)) & + call add_nnlist(coordinate%mesh) + +!Figure out indexing between C and Fortran + call cent_relaxer(mesh_dim(coordinate), node_count(coordinate),& + element_count(coordinate), surface_element_count(coordinate),& + nowned_nodes(coordinate), universal_numbering(coordinate)-1,& + coordinate%mesh%ndglno, coordinate%val, new_coordinate%val,& + initial_coordinate%val, surface_mesh%ndglno,& + coordinate%mesh%adj_lists%nnlist%findrm,& + coordinate%mesh%adj_lists%nnlist%colm) + + + call halo_update(new_coordinate) + + !!! Convert the new coordinates returned into a grid velocity + !!! and calculate the coordinate field at the theta time level. + + call set(grid_velocity,new_coordinate) + call addto(grid_velocity, old_coordinate, scale = -1.0) + call scale(grid_velocity, 1.0/dt) + + if (present_and_true(diagnostic_only)) then + call set(coordinate, old_coordinate) + call set(new_coordinate, old_coordinate) + else + call set(coordinate, new_coordinate, old_coordinate, itheta) + end if + + + end subroutine move_mesh_centroid_relaxer + + subroutine move_mesh_initialise_centroid_relaxer(states) + type(state_type), dimension(:), intent(inout) :: states + + type(vector_field), pointer :: coordinate + type(vector_field) :: initial_coordinate + type(mesh_type) :: surface_mesh + + !!! Store the initial mesh for use as the computational geometry for the + !! Lineal mesh smoothing library. + + coordinate => extract_vector_field(states(1), "Coordinate") + + call allocate(initial_coordinate, mesh=coordinate%mesh, dim=coordinate%dim,& + name="InitialCoordinate") + call set(initial_coordinate, coordinate) + call insert(states(1), initial_coordinate, "InitialCoordinate") + call deallocate(initial_coordinate) + + !!! Add the surface mesh. This needs revisiting for parallel + + call extract_surface_mesh(surface_mesh, coordinate%mesh, "FullCoordinateSurfaceMesh") + call insert(states(1), surface_mesh, "FullCoordinateSurfaceMesh") + call deallocate(surface_mesh) + + end subroutine move_mesh_initialise_centroid_relaxer + + subroutine set_boundary_values(vfield) + type(vector_field), intent(inout) :: vfield + + type(vector_field), pointer:: surface_field + integer, dimension(:), pointer :: surface_node_list + + integer :: i, node + + do i=1,get_boundary_condition_count(vfield) + + call get_boundary_condition(vfield, i, surface_node_list= surface_node_list) + surface_field => extract_surface_field(vfield,i,"value") + + do node=1, node_count(surface_field) + call set(vfield, surface_node_list(node), node_val(surface_field,node)) + end do + + end do + + end subroutine set_boundary_values + + end module meshmovement diff --git a/assemble/Momentum_CG.F90 b/assemble/Momentum_CG.F90 index 42f6a8ec5c..e62a6cffcb 100644 --- a/assemble/Momentum_CG.F90 +++ b/assemble/Momentum_CG.F90 @@ -78,7 +78,8 @@ module momentum_cg public :: construct_momentum_cg, correct_masslumped_velocity, & correct_velocity_cg, assemble_masslumped_poisson_rhs, & add_kmk_matrix, add_kmk_rhs, assemble_kmk_matrix, & - deallocate_cg_mass, assemble_poisson_rhs + deallocate_cg_mass, assemble_poisson_rhs, & + get_stress_free_field ! are we lumping the mass, absorption or source logical :: lump_mass, lump_absorption, lump_source @@ -2775,4 +2776,125 @@ subroutine add_kmk_rhs(state, rhs, pressure, dt) call scale(rhs, dt) end subroutine add_kmk_rhs + + subroutine get_stress_free_field(state, positions, vfield,& + chi, C1, C2, lame1, lame2, option_path,& + ufield) + + type(state_type) :: state + type(vector_field), intent(in) :: positions + type(vector_field), intent(inout) :: vfield + real, intent(in) :: chi, C1, C2, lame1, lame2 + character(len=*), optional, intent(in) :: option_path + type(vector_field), optional, intent(in) :: ufield + + type(petsc_csr_matrix) :: M + + integer :: i, ele + integer, dimension(:), pointer :: nodes + type(element_type), pointer :: u_shape + real, dimension(ele_loc(vfield, 1), ele_ngi(vfield, 1), vfield%dim) :: du_t + real, dimension(vfield%dim, vfield%dim, ele_ngi(vfield, 1)) :: lame1_tensor, lame2_tensor + real, dimension(ele_ngi(vfield, 1)) :: detwei + real, dimension(vfield%dim, ele_ngi(vfield, 1)) :: ugi + type(vector_field) :: rhs + type(csr_sparsity), pointer :: u_sparsity + + lame1_tensor = lame1 + lame2_tensor = lame2 + + u_sparsity => get_csr_sparsity_firstorder(state, vfield%mesh, vfield%mesh) + + call allocate(M, u_sparsity, [vfield%dim, vfield%dim],& + diagonal = .false., name = "StiffnessMatrix") + call allocate(rhs, dim = vfield%dim, mesh= vfield%mesh, name = "RHS") + + call zero(M) + call zero(rhs) + + do ele = 1, element_count(vfield) + + nodes => ele_nodes(vfield, ele) + u_shape=>ele_shape(vfield, ele) + + if (present(ufield)) then + ugi = ele_val_at_quad(ufield,ele) + else + ugi = 0.0 + end if + + call transform_to_physical(positions, ele, & + u_shape, dshape=du_t, detwei=detwei) + + call addto(M, nodes, nodes, stiffness_matrix(du_t, lame1_tensor,& + lame2_tensor, detwei**(1.0- chi))) + + do i=1,vfield%dim + call addto(M, i,i, nodes, nodes, shape_shape(u_shape,u_shape,(C1*detwei)**(1.0- chi))) + end do + call addto(rhs,nodes,shape_vector_rhs(u_shape,ugi,(C2*detwei)**(1.0- chi))) + + end do + + call apply_dirichlet_conditions(matrix=M, rhs=rhs, field=vfield) + + + call petsc_solve(vfield, M, rhs, option_path) + + call deallocate(M) + call deallocate(rhs) + + + contains + function stiffness_matrix(dshape, l1, l2, detwei) result (matrix) + !!< Calculates the element stiffness matrix. + + real, dimension(:,:,:), intent(in) :: dshape + real, dimension(size(dshape,3),size(dshape,3),size(dshape,2)), intent(in) :: l1, l2 + real, dimension(size(dshape,2)), intent(in) :: detwei + + real, dimension(size(dshape,3),size(dshape,3),size(dshape,1),size(dshape,1)) :: matrix + + real, dimension(size(dshape,3),size(dshape,2)) :: tensor_diag, tensor_entries + + integer :: iloc,jloc, gi, i, j + integer :: loc, ngi, dim + + loc=size(dshape,1) + ngi=size(dshape,2) + dim=size(dshape,3) + + tensor_diag = 0.0 + tensor_entries = 0.0 + + matrix =0.0 + + do i=1,dim + ! extract the relevent tensor entries into a vector + do j = 1, i-1 + tensor_entries(j,:) = l2(j,i,:) + end do + do j = i, dim + tensor_entries(j,:) = l2(i,j,:) + end do + matrix(i,i,:,:) = dshape_vector_dshape(dshape, 2.0*tensor_entries, dshape, detwei) + end do + + do gi=1,ngi + forall(iloc=1:loc,jloc=1:loc) + matrix(:,:,iloc,jloc) = matrix(:,:,iloc,jloc) & + +(spread(dshape(iloc,gi,:), 1, dim) & + *spread(dshape(jloc,gi,:), 2, dim) & + *2.0*l2(:,:,gi) & + +spread(dshape(iloc,gi,:), 2, dim) & + *spread(dshape(jloc,gi,:), 1, dim) & + *l1(:,:,gi)) & + *detwei(gi) + end forall + end do + + end function stiffness_matrix + + end subroutine get_stress_free_field + end module momentum_cg diff --git a/assemble/centroid_relaxer.c b/assemble/centroid_relaxer.c new file mode 100644 index 0000000000..4c51814e3b --- /dev/null +++ b/assemble/centroid_relaxer.c @@ -0,0 +1,264 @@ +static char help[] = "2D serial centroid relaxer.\n\n"; +#include +#include +#include +#include + + +void cent_relaxer(int dimension, int num_nodes, int num_elements, int num_surf_elements, int num_owned_nodes, int * mapping, int connectivity[][3], double * phys_mesh, double * smooth_mesh, double * comp_mesh, int * surf_connectivity, int* findrm, int* colm) { + + Vec F,U_h; + Mat K; + KSP ksp; + PC pc; + PetscErrorCode ierr; + PetscInt i,n,m,j,its,nodei,Ii,len_new=0,local_surf_connectivity[dimension*num_surf_elements], global_surf_connectivity[dimension*num_surf_elements], + num_nodes_col_x[num_nodes],num_nodes_col_y[num_nodes],num_nodes_col_z[num_nodes]; + PetscScalar length,x_disp[num_surf_elements],y_disp[num_surf_elements],z_disp[num_surf_elements],smoothed_x[num_nodes],smoothed_y[num_nodes],smoothed_z[num_nodes]; + + double lij(int nodei, int nodej){ + double xi,yi,zi,xj,yj,zj,length; + if (nodei == nodej){ + PetscPrintf(PETSC_COMM_WORLD,"Failure: nodei is the same as nodej.\n"); + length = 0.0; + } + else if (dimension == 1) { + xi = comp_mesh[dimension*nodei-1]; + xj = comp_mesh[dimension*nodej-1]; + length = labs(xi-xj); + } + else if (dimension == 2) { + xi = comp_mesh[dimension*nodei]; yi = comp_mesh[dimension*nodei+1]; + xj = comp_mesh[dimension*nodej]; yj = comp_mesh[dimension*nodej+1]; + length = sqrt(pow(xi-xj,2)+pow(yi-yj,2)); + } + else { + xi = comp_mesh[dimension*nodei];yi = comp_mesh[dimension*nodei+1];zi = comp_mesh[dimension*nodei+2]; + xj = comp_mesh[dimension*nodej];yj = comp_mesh[dimension*nodej+1];zj = comp_mesh[dimension*nodej+2]; + length = sqrt(pow(xi-xj,2)+pow(yi-yj,2)+pow(zi-zj,2)); + } + + return length; + } + PetscPrintf(PETSC_COMM_WORLD,"Inside cent_relaxer.\n"); + /* + MatCreate(PETSC_COMM_WORLD,&K); + PetscObjectSetName((PetscObject) K, "Stiffness Matrix"); + MatSetSizes(K,dimension*num_owned_nodes,dimension*num_owned_nodes,PETSC_DETERMINE,PETSC_DETERMINE); + MatSetUp(K); + + for(Ii=0;Iinum_owned_nodes){ + continue;} + + for(j=0;j +#include +#include +#include +#include "solver_options.h" + +void lap_smoother(int dimension, int num_nodes, int num_elements, int num_surf_elements, int num_owned_nodes, int* mapping, int * connectivity, double * phys_mesh, double* smooth_mesh, double * comp_mesh, int * surf_connectivity, Mat* K, struct solver_options* options, int debug_level) { + + Vec F[3],U_h[3]; + PetscInt i,n,m,convert[4],*holder_new, + Ii,len_new=0,j,local_surf_connectivity[dimension*num_surf_elements], global_surf_connectivity[dimension*num_surf_elements]; + PetscScalar p_value[4],A_holder[16],grad_holder[12],det,Fe,x_bound_points[num_surf_elements],Vol,A_inv_holder[4][4],length, + y_bound_points[num_surf_elements],z_bound_points[num_surf_elements],Area,Ke[4][4]; + + //Decomposing F into F_x, F_y and F_z owing to nature of Laplace's equation (independence of variables) + VecCreate(PETSC_COMM_WORLD,&F[0]); + PetscObjectSetName((PetscObject) F[0], "Fx"); + VecSetSizes(F[0],num_owned_nodes,PETSC_DECIDE); + VecSetFromOptions(F[0]); + if (dimension>1) { + VecCreate(PETSC_COMM_WORLD,&F[1]); + PetscObjectSetName((PetscObject) F[1], "Fy"); + VecSetSizes(F[1],num_owned_nodes,PETSC_DECIDE); + VecSetFromOptions(F[1]); + if (dimension>2) { + VecCreate(PETSC_COMM_WORLD,&F[2]); + PetscObjectSetName((PetscObject) F[2], "Fz"); + VecSetSizes(F[2],num_owned_nodes,PETSC_DECIDE); + VecSetFromOptions(F[2]); + } + } + + + //Assembly of Stiffness Matrix K + if (dimension == 1){ + for(Ii=0;Iinum_owned_nodes) continue; + + for(j=0;j< len_new;j++){ + + if(surf_connectivity[i] == local_surf_connectivity[j]+1) + break; + } + + if (j==len_new) + local_surf_connectivity[len_new++] = surf_connectivity[i]-1; + } + + for(i=0;i1) { + y_bound_points[n]=phys_mesh[dimension*(local_surf_connectivity[n])+1]; + if (dimension>2) { + z_bound_points[n]=phys_mesh[dimension*(local_surf_connectivity[n])+2]; + } + } + } + + for(Ii=0;Ii1) { + VecSetValue(F[1],global_surf_connectivity[Ii],y_bound_points[Ii],INSERT_VALUES); + if (dimension>2) { + VecSetValue(F[2],global_surf_connectivity[Ii],z_bound_points[Ii],INSERT_VALUES); + } + } + } + } + + MatAssemblyBegin(*K,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(*K,MAT_FINAL_ASSEMBLY); + + + for (j=0; j +#include +#include +#include +#include "solver_options.h" + + +void lin_smoother(int dimension, int num_nodes, int num_elements, int num_surf_elements, int num_owned_nodes, int * mapping, int connectivity[][3], double * phys_mesh, double * smooth_mesh, double * comp_mesh, int * surf_connectivity, int* findrm, int* colm, Mat* K, struct solver_options* options, int debug_level) { + + Vec F,U_h; + PetscInt i,n,j,its,Ii,len_new=0,local_surf_connectivity[dimension*num_surf_elements], global_surf_connectivity[dimension*num_surf_elements]; + PetscScalar x_disp[num_surf_elements],y_disp[num_surf_elements],z_disp[num_surf_elements]; + + MatSetOption(*K,MAT_SYMMETRIC,PETSC_TRUE); + + double lij(int nodei, int nodej){ + double xi,yi,zi,xj,yj,zj,length; + + if (nodei == nodej){ + PetscPrintf(PETSC_COMM_WORLD,"Failure: nodei is the same as nodej.\n"); + length = 0.0; + } + else if (dimension == 1) { + xi = comp_mesh[dimension*nodei-1]; + xj = comp_mesh[dimension*nodej-1]; + length = labs(xi-xj); + } + else if (dimension == 2) { + xi = comp_mesh[dimension*nodei]; yi = comp_mesh[dimension*nodei+1]; + xj = comp_mesh[dimension*nodej]; yj = comp_mesh[dimension*nodej+1]; + length = sqrt(pow(xi-xj,2)+pow(yi-yj,2)); + } + else { + xi = comp_mesh[dimension*nodei];yi = comp_mesh[dimension*nodei+1];zi = comp_mesh[dimension*nodei+2]; + xj = comp_mesh[dimension*nodej];yj = comp_mesh[dimension*nodej+1];zj = comp_mesh[dimension*nodej+2]; + length = sqrt(pow(xi-xj,2)+pow(yi-yj,2)+pow(zi-zj,2)); + } + + return length; + } + + /*Lineal Stiffness*/ + double * K_lin(int nodei, int nodej){ + + double xi,yi,xj,yj,x_edge,y_edge,alpha,inv_length; + + xi = comp_mesh[dimension*nodei];yi = comp_mesh[dimension*nodei+1]; + xj = comp_mesh[dimension*nodej];yj = comp_mesh[dimension*nodej+1]; + x_edge = xj - xi;y_edge = yj - yi; + alpha = atan2(y_edge,x_edge);inv_length = 1.0/lij(nodei,nodej); + + if(dimension == 2){ + static double K_lin_temp[3]; + K_lin_temp[0]=-1.0*inv_length*cos(alpha)*cos(alpha); + K_lin_temp[1]=-1.0*inv_length*cos(alpha)*sin(alpha); + K_lin_temp[2]=-1.0*inv_length*sin(alpha)*sin(alpha); + + return K_lin_temp; + } + else{ + double zi,zj,z_edge,phi; + zi=comp_mesh[dimension*nodei+2];zj=comp_mesh[dimension*nodei+2];z_edge=zj-zi;phi=acos(z_edge/lij(nodei,nodej)); + static double K_lin_temp_3d[8]; + + + K_lin_temp_3d[0]=-1.0*inv_length*cos(alpha)*cos(alpha)*sin(phi); + K_lin_temp_3d[1]=-1.0*inv_length*sin(alpha)*cos(alpha)*sin(phi); + K_lin_temp_3d[2]=-1.0*inv_length*cos(alpha)*cos(phi); + K_lin_temp_3d[3]=-1.0*inv_length*sin(alpha)*sin(alpha)*sin(phi); + K_lin_temp_3d[4]=-1.0*inv_length*sin(alpha)*cos(phi); + K_lin_temp_3d[5]=-1.0*inv_length*cos(alpha)*sin(phi)*cos(phi); + K_lin_temp_3d[6]=-1.0*inv_length*sin(alpha)*sin(phi)*cos(phi); + K_lin_temp_3d[7]=-1.0*inv_length*cos(phi)*cos(phi); + + return K_lin_temp_3d; + + } + } + + for(Ii=0;Iinum_owned_nodes){ + continue;} + + for(j=0;j +#include +#include +#include +#include "solver_options.h" + + +void lin_tor_smoother(int dimension, int num_nodes, int num_elements, int num_surf_elements, int num_owned_nodes, int * mapping, int * connectivity, double * phys_mesh, double * smooth_mesh, double * comp_mesh, int * surf_connectivity, int* findrm, int* colm, Mat* K, struct solver_options* options, int debug_level) { + + Vec F,U_h; + PetscInt i,n,j,Ii,len_new=0,local_surf_connectivity[dimension*num_surf_elements], global_surf_connectivity[dimension*num_surf_elements]; + PetscScalar x_disp[num_surf_elements],y_disp[num_surf_elements],z_disp[num_surf_elements]; + + MatSetOption(*K,MAT_SYMMETRIC,PETSC_TRUE); + double lij(int nodei, int nodej){ + double xi,yi,zi,xj,yj,zj,length; + if (nodei == nodej){ + PetscPrintf(PETSC_COMM_WORLD,"Failure: nodei is the same as nodej.\n"); + length = 0.0; + } + else if (dimension == 1) { + xi = comp_mesh[dimension*nodei-1]; + xj = comp_mesh[dimension*nodej-1]; + length = labs(xi-xj); + } + else if (dimension == 2) { + xi = comp_mesh[dimension*nodei]; yi = comp_mesh[dimension*nodei+1]; + xj = comp_mesh[dimension*nodej]; yj = comp_mesh[dimension*nodej+1]; + length = sqrt(pow(xi-xj,2)+pow(yi-yj,2)); + } + else { + xi = comp_mesh[dimension*nodei];yi = comp_mesh[dimension*nodei+1];zi = comp_mesh[dimension*nodei+2]; + xj = comp_mesh[dimension*nodej];yj = comp_mesh[dimension*nodej+1];zj = comp_mesh[dimension*nodej+2]; + length = sqrt(pow(xi-xj,2)+pow(yi-yj,2)+pow(zi-zj,2)); + } + + return length; + } + + double alpha_ang(int nodei, int nodej){ + double alpha; + double xi = comp_mesh[dimension*nodei]; double yi = comp_mesh[dimension*nodei+1]; + double xj = comp_mesh[dimension*nodej]; double yj = comp_mesh[dimension*nodej+1]; + double x_edge = xj - xi; double y_edge = yj - yi; + alpha = atan2(y_edge,x_edge); + return alpha; + } + + /*Lineal Stiffness*/ + double * K_lin(int nodei, int nodej){ + double xi,yi,xj,yj,x_edge,y_edge,alpha,inv_length; + + xi = comp_mesh[dimension*nodei];yi = comp_mesh[dimension*nodei+1]; + xj = comp_mesh[dimension*nodej];yj = comp_mesh[dimension*nodej+1]; + x_edge = xj - xi;y_edge = yj - yi; + alpha = atan2(y_edge,x_edge);inv_length = 1.0/lij(nodei,nodej); + + static double K_lin_temp[3]; + K_lin_temp[0]=-1.0*inv_length*cos(alpha)*cos(alpha); + K_lin_temp[1]=-1.0*inv_length*cos(alpha)*sin(alpha); + K_lin_temp[2]=-1.0*inv_length*sin(alpha)*sin(alpha); + + return K_lin_temp; + } + + double area(int nodei, int nodej, int nodek){ + double s = (lij(nodei,nodej) + lij(nodei,nodek) + lij(nodej,nodek))/2.0; + double area = sqrt(s*(s-lij(nodei,nodej))*(s-lij(nodei,nodek))*(s-lij(nodej,nodek))); + return area; + } + + double tors(int nodei, int nodej, int nodek){ + double tors = (pow(lij(nodei,nodej),2) * pow(lij(nodei,nodek),2))/(4*pow(area(nodei,nodej,nodek),2)); + return tors; + } + + double * C(int nodei, int nodej, int nodek){ + double* C_temp; + C_temp = malloc(sizeof(double)*9); + + C_temp[0]=tors(nodei,nodej,nodek); + C_temp[1]=0;C_temp[2]=0;C_temp[3]=0; + C_temp[4]=tors(nodej,nodek,nodei); + C_temp[5]=0;C_temp[6]=0;C_temp[7]=0; + C_temp[8]=tors(nodek,nodei,nodej); + return C_temp; + } + + double x_edge(int nodei,int nodej){ + double xi,xj,x_edge; + xi = comp_mesh[dimension*nodei]; + xj = comp_mesh[dimension*nodej]; + x_edge = xj - xi; + return x_edge; + } + + double y_edge(int nodei,int nodej){ + double yi,yj,y_edge; + yi = comp_mesh[dimension*nodei+1]; + yj = comp_mesh[dimension*nodej+1]; + y_edge = yj - yi; + return y_edge; + } + + double a_rot(int nodei, int nodej){ + double a_rot = x_edge(nodei,nodej)/(pow(lij(nodei,nodej),2)); + return a_rot; + } + + double b_rot(int nodei, int nodej){ + double b_rot = y_edge(nodei,nodej)/(pow(lij(nodei,nodej),2)); + return b_rot; + } + + double * rot(int nodei, int nodej, int nodek){ + double* rot_mat; + rot_mat = malloc(sizeof(double)*18); + + rot_mat[0]=b_rot(nodei,nodek)-b_rot(nodei,nodej); + rot_mat[1]=a_rot(nodei,nodej)-a_rot(nodei,nodek); + rot_mat[2]=b_rot(nodei,nodej); + rot_mat[3]=-1.0*a_rot(nodei,nodej); + rot_mat[4]=-1.0*b_rot(nodei,nodek); + rot_mat[5]=a_rot(nodei,nodek); + + rot_mat[6]=-1.0*b_rot(nodej,nodei); + rot_mat[7]=a_rot(nodej,nodei); + rot_mat[8]=b_rot(nodej,nodei)-b_rot(nodej,nodek); + rot_mat[9]=a_rot(nodej,nodek)-a_rot(nodej,nodei); + rot_mat[10]=b_rot(nodej,nodek); + rot_mat[11]=-1.0*a_rot(nodej,nodek); + + rot_mat[12]=b_rot(nodek,nodei); + rot_mat[13]=-1.0*a_rot(nodek,nodei); + rot_mat[14]=-1.0*b_rot(nodek,nodej); + rot_mat[15]=a_rot(nodek,nodej); + rot_mat[16]=b_rot(nodek,nodej)-b_rot(nodek,nodei); + rot_mat[17]=a_rot(nodek,nodei)-a_rot(nodek,nodej); + + return rot_mat; + } + + double * K_tor(int nodei, int nodej, int nodek){ + int i,j,k,r1=6,c1=3,c2=3; + double rot_mat[3][6],rot_trans_mat[6][3], C_mat[3][3]; + double *A_mat,*K_tor_mat,*rot_hold, *C_hold,*K_tor_mat_trim; + rot_hold = rot(nodei,nodej,nodek); + C_hold = C(nodei,nodej,nodek); + + K_tor_mat = malloc(sizeof(double)*36); + K_tor_mat_trim = malloc(sizeof(double)*24); + A_mat = malloc(sizeof(double)*18); + + rot_mat[0][0]=rot_hold[0]; rot_mat[0][1]=rot_hold[1]; rot_mat[0][2]=rot_hold[2]; rot_mat[0][3]=rot_hold[3]; rot_mat[0][4]=rot_hold[4]; rot_mat[0][5]=rot_hold[5]; + rot_mat[1][0]=rot_hold[6]; rot_mat[1][1]=rot_hold[7]; rot_mat[1][2]=rot_hold[8]; rot_mat[1][3]=rot_hold[9]; rot_mat[1][4]=rot_hold[10];rot_mat[1][5]=rot_hold[11]; + rot_mat[2][0]=rot_hold[12];rot_mat[2][1]=rot_hold[13];rot_mat[2][2]=rot_hold[14];rot_mat[2][3]=rot_hold[15];rot_mat[2][4]=rot_hold[16];rot_mat[2][5]=rot_hold[17]; + + rot_trans_mat[0][0]=rot_hold[0]; rot_trans_mat[0][1]=rot_hold[6]; rot_trans_mat[0][2]=rot_hold[12]; + rot_trans_mat[1][0]=rot_hold[1]; rot_trans_mat[1][1]=rot_hold[7]; rot_trans_mat[1][2]=rot_hold[13]; + rot_trans_mat[2][0]=rot_hold[2]; rot_trans_mat[2][1]=rot_hold[8]; rot_trans_mat[2][2]=rot_hold[14]; + rot_trans_mat[3][0]=rot_hold[3]; rot_trans_mat[3][1]=rot_hold[9]; rot_trans_mat[3][2]=rot_hold[15]; + rot_trans_mat[4][0]=rot_hold[4]; rot_trans_mat[4][1]=rot_hold[10];rot_trans_mat[4][2]=rot_hold[16]; + rot_trans_mat[5][0]=rot_hold[5]; rot_trans_mat[5][1]=rot_hold[11];rot_trans_mat[5][2]=rot_hold[17]; + + C_mat[0][0]=C_hold[0]; C_mat[0][1]=C_hold[1]; C_mat[0][2]=C_hold[2]; + C_mat[1][0]=C_hold[3]; C_mat[1][1]=C_hold[4]; C_mat[1][2]=C_hold[5]; + C_mat[2][0]=C_hold[6]; C_mat[2][1]=C_hold[7]; C_mat[2][2]=C_hold[8]; + + for(i=0;i=num_owned_nodes){continue;} + int neb_hold = connectivity[3*ele+(face+1)%3]-1; + int k1 = connectivity[3*ele+(face+2)%3]-1; + + if (Iinum_owned_nodes){ + continue;} + + for(j=0;j extract_vector_field(state, "Coordinate") + velocity => extract_vector_field(state, "IteratedVelocity") + + ! Allocate strain_rate and vorticity: + call allocate(grad_vel, s_field%mesh, name="grad_vel") + + ! Calculate strain rate second invariant: + call grad(velocity,positions,grad_vel) + + do i=1, node_count(s_field) + t=node_val(grad_vel,i) + trace=0 + do j=1,velocity%dim + trace=trace+t(j,j) + end do + do j=1,velocity%dim + t(j,j)=t(j,j)-trace/velocity%dim + end do + omega2=sum(((t-transpose(t))/2.)**2) + S2=sum(((t+transpose(t))/2.) **2) + if (normalise) then + call set(s_field,i, (omega2-s2)/(tol+omega2+s2)) + else + call set(s_field,i, (omega2-s2)/2.0) + end if + end do + + ! Clean-up: + call deallocate(grad_vel) + + ! Print min and max: + ewrite_minmax(s_field) + + + end subroutine calculate_q_criterion + subroutine calculate_sediment_concentration_dependent_viscosity(state, t_field) ! calculates viscosity based upon total sediment concentration type(state_type), intent(inout) :: state @@ -597,4 +656,36 @@ subroutine calculate_geostrophic_velocity(state, v_field) end subroutine calculate_geostrophic_velocity + subroutine calculate_stress_free_field(state,v_field) + type(state_type), intent(inout) :: state + type(vector_field), intent(inout) :: v_field + + character(len = OPTION_PATH_LEN) :: path + character(len = FIELD_NAME_LEN) :: name + type(vector_field), pointer :: positions, velocity + real :: chi, C1, C2, lame1, lame2 + + path = trim(complete_field_path(v_field%option_path)) // "/algorithm" + + call get_option(trim(path)//"/chi_parameter", chi, default=0.0) + call get_option(trim(path)//"/lambda_parameter", lame1, default=1.0) + call get_option(trim(path)//"/mu_parameter", lame2, default=1.0) + call get_option(trim(path)//"/C1", C1, default=0.0) + call get_option(trim(path)//"/C2", C2, default=0.0) + + ! Extract positions field from state: + if (have_option(trim(path)//"/coordinate_field")) then + call get_option(trim(path)//"/coordinate_field/name", name) + else + name = "Coordinate" + end if + positions => extract_vector_field(state, trim(name)) + velocity => vector_source_field(state, v_field) + + call get_stress_free_field(state, positions, v_field,& + chi, C1, C2, lame1, lame2, path, velocity) + + end subroutine calculate_stress_free_field + + end module momentum_diagnostics diff --git a/femtools/C_Solvers.c b/femtools/C_Solvers.c new file mode 100644 index 0000000000..1cb4a0ed55 --- /dev/null +++ b/femtools/C_Solvers.c @@ -0,0 +1,105 @@ +/* ! Copyright (C) 2006 Imperial College London and others. */ +/* ! */ +/* ! Please see the AUTHORS file in the main source directory for a full list */ +/* ! of copyright holders. */ +/* ! */ +/* ! Prof. C Pain */ +/* ! Applied Modelling and Computation Group */ +/* ! Department of Earth Science and Engineering */ +/* ! Imperial College London */ +/* ! */ +/* ! amcgsoftware@imperial.ac.uk */ +/* ! */ +/* ! This library is free software; you can redistribute it and/or */ +/* ! modify it under the terms of the GNU Lesser General Public */ +/* ! License as published by the Free Software Foundation, */ +/* ! version 2.1 of the License. */ +/* ! */ +/* ! This library is distributed in the hope that it will be useful, */ +/* ! but WITHOUT ANY WARRANTY; without even the implied warranty of */ +/* ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */ +/* ! Lesser General Public License for more details. */ +/* ! */ +/* ! You should have received a copy of the GNU Lesser General Public */ +/* ! License along with this library; if not, write to the Free Software */ +/* ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 */ +/* ! USA */ + +#include +#include +#include +#include +#include "solver_options.h" + +void setup_ksp_from_options_c(KSP ksp, struct solver_options *options) { + +// Initialise a PETSc KSP solver from a set of solver options +// +// This routine duplicates much of the functionality of the equivalent +// Fortran subroutine + + PC pc; + + KSPSetType(ksp,options->ksptype); + KSPGetPC(ksp,&pc); + + int size; + MPI_Comm_size(PETSC_COMM_WORLD,&size); + if (strcmp(options->pctype,PCHYPRE)) { + PCSetType(pc,options->pctype); + PCHYPRESetType(pc, options->hypretype);} + else if (strcmp(options->pctype,PCSOR) + || strcmp(options->pctype,PCEISENSTAT)) { + if (size == 1) { + PCSetType(pc,options->pctype); + } else { + PCSetType(pc, PCBJACOBI); + PCSetUp(pc); + KSP *subksps; + PC subpc; + PCBJacobiGetSubKSP(pc, NULL, NULL, &subksps); + KSPGetPC(subksps[0], &subpc); + PCSetType(subpc,options->pctype); + } + } + KSPSetInitialGuessNonzero(ksp,~options->start_from_zero); + KSPSetTolerances(ksp,options->rtol,options->atol,PETSC_DEFAULT,options->max_its); + KSPSetFromOptions(ksp); +} + +void petsc_solve_many_c(Vec *x, Mat M, Vec *b, int neqns, + struct solver_options *options, + int debug_level) { + KSP ksp; + int n, its; + KSPCreate(PETSC_COMM_WORLD,&ksp); + + KSPSetOperators(ksp,M,M); + + setup_ksp_from_options_c(ksp, options); + + char profiler_title[80]; + const char *matname; + PetscObjectGetName((PetscObject)(x[0]),&matname); + strncpy(profiler_title,matname,80); + strncat(profiler_title,"::solve",80); + + cprofiler_tic(profiler_title); + for (n=0;n1) PetscPrintf(PETSC_COMM_WORLD,"ksp iter: %D\n", its); + KSPConvergedReason reason; + KSPGetConvergedReason(ksp, &reason); + if (debug_level>1) PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason); + + KSPDestroy(&ksp); +} + +void petsc_solve_c(Vec x, Mat M, Vec b, + struct solver_options *options, + int debug_level) { + petsc_solve_many_c(&x,M,&b,1,options, debug_level); +} diff --git a/femtools/Diagnostic_Fields.F90 b/femtools/Diagnostic_Fields.F90 index 0a13a321df..f7264eaac1 100644 --- a/femtools/Diagnostic_Fields.F90 +++ b/femtools/Diagnostic_Fields.F90 @@ -31,7 +31,7 @@ module diagnostic_fields !!< A module to calculate diagnostic fields. use fldebug - use global_parameters, only:FIELD_NAME_LEN, current_time, OPTION_PATH_LEN + use global_parameters, only:FIELD_NAME_LEN, current_time, OPTION_PATH_LEN, timestep use futils use spud use Vector_Tools @@ -322,6 +322,9 @@ subroutine calculate_vector_diagnostic_variable_single_state(state, d_field_name case("DiagnosticCoordinate") call calculate_diagnostic_coordinate_field(state, d_field) + case("InitialCoordinate") + call calculate_initial_coordinate_field(state, d_field) + case default if(present(stat)) then stat = 1 @@ -410,9 +413,9 @@ subroutine calculate_CFL_number(State, CFL, dt) type(scalar_field), intent(inout) :: cfl real, intent(in), optional :: dt - type(vector_field), pointer :: U, X + type(vector_field), pointer :: U, X, ugrid real :: l_dt - integer :: ele, gi + integer :: ele, gi, stat ! Transformed quadrature weights. real, dimension(ele_ngi(CFL, 1)) :: detwei ! Inverse of the local coordinate change matrix. @@ -426,8 +429,13 @@ subroutine calculate_CFL_number(State, CFL, dt) ! current CFL element shape type(element_type), pointer :: CFL_shape + logical :: move_mesh + U=>extract_vector_field(state, "Velocity") X=>extract_vector_field(state, "Coordinate") + + move_mesh = have_option("/mesh_adaptivity/mesh_movement") + Ugrid=>extract_vector_field(state, "GridVelocity", stat) if(present(dt)) then l_dt = dt @@ -449,6 +457,9 @@ subroutine calculate_CFL_number(State, CFL, dt) ! be. I don't understand why it's this way round but the results ! appear correct. -dham CFL_q=ele_val_at_quad(U, ele) + + if (move_mesh) CFL_q=CFL_q-ele_val_at_quad(ugrid, ele) + do gi=1, size(detwei) CFL_q(:,gi)=l_dt*matmul(CFL_q(:,gi), invJ(:,:,gi)) end do @@ -1611,9 +1622,10 @@ subroutine calculate_courant_number_dg(state, courant, dt) type(scalar_field), intent(inout) :: courant real, intent(in), optional :: dt ! - type(vector_field), pointer :: u, x + type(vector_field), pointer :: u, x, ugrid real :: l_dt integer :: ele, stat + logical :: grid_vel u=>extract_vector_field(state, "NonlinearVelocity",stat) if(stat.ne.0) then @@ -1622,6 +1634,8 @@ subroutine calculate_courant_number_dg(state, courant, dt) FLExit('Missing velocity field!') end if end if + grid_vel = have_option("/mesh_adaptivity/mesh_movement") + ugrid=>extract_vector_field(state, "GridVelocity",stat) x=>extract_vector_field(state, "Coordinate") if(present(dt)) then @@ -1633,7 +1647,7 @@ subroutine calculate_courant_number_dg(state, courant, dt) call zero(courant) do ele = 1, element_count(courant) - call calculate_courant_number_dg_ele(courant,x,u,ele,l_dt) + call calculate_courant_number_dg_ele(courant,x,u,ugrid,ele,l_dt,grid_vel) end do ! the courant values at the edge of the halo are going to be incorrect @@ -1642,11 +1656,12 @@ subroutine calculate_courant_number_dg(state, courant, dt) end subroutine calculate_courant_number_dg - subroutine calculate_courant_number_dg_ele(courant,x,u,ele,dt) - type(vector_field), intent(in) :: x, u + subroutine calculate_courant_number_dg_ele(courant,x,u,ugrid,ele,dt,grid_vel) + type(vector_field), intent(in) :: x, u, ugrid type(scalar_field), intent(inout) :: courant real, intent(in) :: dt integer, intent(in) :: ele + logical :: grid_vel ! real :: Vol real :: Flux @@ -1679,6 +1694,10 @@ subroutine calculate_courant_number_dg_ele(courant,x,u,ele,dt) U_f_quad =0.5*(face_val_at_quad(U, face)& & +face_val_at_quad(U, face_2)) + if (grid_vel) then + u_f_quad=u_f_quad-face_val_at_quad(ugrid,face) + end if + call transform_facet_to_physical(X, face, & & detwei_f=detwei_f,& & normal=normal) @@ -3808,5 +3827,18 @@ subroutine calculate_diagnostic_coordinate_field(state, field) call deallocate(coordinate_field) end subroutine calculate_diagnostic_coordinate_field + + subroutine calculate_initial_coordinate_field(state, field) + type(state_type), intent(in) :: state + type(vector_field), intent(inout) :: field + + type(vector_field), pointer :: coordinate_field + + if (timestep == 0) then + coordinate_field => extract_vector_field(state, "Coordinate") + call set(field, coordinate_field) + end if + + end subroutine calculate_initial_coordinate_field end module diagnostic_fields diff --git a/femtools/Fields_Allocates.F90 b/femtools/Fields_Allocates.F90 index 1cb0a259ac..f37e1569e2 100644 --- a/femtools/Fields_Allocates.F90 +++ b/femtools/Fields_Allocates.F90 @@ -62,6 +62,7 @@ module fields_allocates public :: add_lists, extract_lists, add_nnlist, extract_nnlist, add_nelist, & & extract_nelist, add_eelist, extract_eelist, remove_lists, remove_nnlist, & & remove_nelist, remove_eelist, extract_elements, remove_boundary_conditions + public :: extract_surface_mesh interface allocate module procedure allocate_scalar_field, allocate_vector_field,& @@ -2111,6 +2112,52 @@ subroutine create_surface_mesh(surface_mesh, surface_nodes, & end if end subroutine create_surface_mesh + + subroutine extract_surface_mesh(surface_mesh, mesh, name) + !! Creates a surface mesh consisting of all surface elements of a mesh + !! The mesh uses the global numbering + type(mesh_type), intent(out):: surface_mesh + !! mesh to take surface from (should have %faces component) + type(mesh_type), intent(in):: mesh + !! name for the new surface_mesh + character(len=*), intent(in):: name + + integer, dimension(:), pointer:: suf_ndglno + integer, dimension(:), allocatable:: nod2sufnod + integer, dimension(mesh%faces%shape%loc):: glnodes + integer j, sele, sufnod, snloc + + snloc=mesh%faces%shape%loc + + allocate(nod2sufnod(1:node_count(mesh))) + nod2sufnod=0 + + ! mark surface nodes with nod2sufnod(nod)==1 + sufnod=0 + do sele=1, surface_element_count(mesh) + glnodes=face_global_nodes(mesh, sele) + do j=1, snloc + if (nod2sufnod(glnodes(j))==0) then + sufnod=sufnod+1 + nod2sufnod(glnodes(j))=1 + end if + end do + end do + + call allocate(surface_mesh, nodes=sufnod, & + elements=surface_element_count(mesh), shape=mesh%faces%shape, & + name=name) + + surface_mesh%periodic=mesh%periodic + surface_mesh%continuity=mesh%continuity + + ! map global node numbering to surface node numbering + suf_ndglno => surface_mesh%ndglno + do sele=1, surface_element_count(mesh) + suf_ndglno( (sele-1)*snloc+1:sele*snloc )=face_global_nodes(mesh, sele) + end do + + end subroutine extract_surface_mesh logical function SetContains(a, b) !!< Auxillary function that returns true if b contains a diff --git a/femtools/Makefile.in b/femtools/Makefile.in index 6cb7a2b82d..725d1da42e 100644 --- a/femtools/Makefile.in +++ b/femtools/Makefile.in @@ -107,7 +107,7 @@ OBJS = Dgtools.o Coordinates.o EventCounter.o \ GMSH_Common.o Read_GMSH.o Write_GMSH.o \ Exodusii_C_Interface.o Exodusii_F_Interface.o Exodusii_Common.o Read_Exodusii.o \ Mesh_Files.o Vertical_Extrapolation.o \ - Mesh_Quality.o Mesh_Quality_C.o + Mesh_Quality.o Mesh_Quality_C.o C_Solvers.o # objects to be included in libfemtools: diff --git a/femtools/Parallel_fields.F90 b/femtools/Parallel_fields.F90 index 68d50fa2fe..9568d9a90d 100644 --- a/femtools/Parallel_fields.F90 +++ b/femtools/Parallel_fields.F90 @@ -53,7 +53,7 @@ module parallel_fields & element_owner, node_owned, assemble_ele, & & surface_element_owned, nowned_nodes ! Apparently ifort has a problem with the generic name node_owned - public :: node_owned_mesh, zero_non_owned + public :: node_owned_mesh, zero_non_owned, universal_numbering interface node_owned module procedure node_owned_mesh, node_owned_scalar, node_owned_vector, & @@ -100,6 +100,12 @@ module parallel_fields & nowned_nodes_vector, nowned_nodes_tensor end interface nowned_nodes + interface universal_numbering + module procedure universal_numbering_mesh, universal_numbering_scalar,& + universal_numbering_vector, universal_numbering_tensor + end interface universal_numbering + + contains function halo_communicator_mesh(mesh) result(communicator) @@ -757,4 +763,48 @@ pure function nowned_nodes_tensor(t_field) result(nodes) end function nowned_nodes_tensor + + function universal_numbering_mesh(mesh) result(map) + type(mesh_type), intent(in) :: mesh + + integer, dimension(node_count(mesh)) :: map + + integer :: nhalos, i + + nhalos = halo_count(mesh) + if(nhalos > 0) then + call get_universal_numbering(mesh%halos(nhalos), map) + else + map = [(i,i=1,node_count(mesh))] + end if + + end function universal_numbering_mesh + + function universal_numbering_scalar(s_field) result(map) + type(scalar_field), intent(in) :: s_field + + integer, dimension(node_count(s_field)) :: map + + map = universal_numbering_mesh(s_field%mesh) + + end function universal_numbering_scalar + + function universal_numbering_vector(v_field) result(map) + type(vector_field), intent(in) :: v_field + + integer, dimension(node_count(v_field)) :: map + + map = universal_numbering_mesh(v_field%mesh) + + end function universal_numbering_vector + + function universal_numbering_tensor(t_field) result(map) + type(tensor_field), intent(in) :: t_field + + integer, dimension(node_count(t_field)) :: map + + map = universal_numbering_mesh(t_field%mesh) + + end function universal_numbering_tensor + end module parallel_fields diff --git a/femtools/Profiler.cpp b/femtools/Profiler.cpp index 5b5adb9b6f..2b0f61f11e 100644 --- a/femtools/Profiler.cpp +++ b/femtools/Profiler.cpp @@ -141,11 +141,18 @@ extern "C" { void cprofiler_tic_fc(const char *key, const int *key_len){ flprofiler.tic(string(key, *key_len)); } + void cprofiler_tic(const char *key){ + flprofiler.tic(string(key, strlen(key))); + } + #define cprofiler_toc_fc F77_FUNC(cprofiler_toc, CPROFILER_TOC) void cprofiler_toc_fc(const char *key, const int *key_len){ flprofiler.toc(string(key, *key_len)); } + void cprofiler_toc(const char *key){ + flprofiler.toc(string(key, strlen(key))); + } #define cprofiler_zero_fc F77_FUNC(cprofiler_zero, CPROFILER_ZERO) void cprofiler_zero_fc(){ diff --git a/femtools/Solvers.F90 b/femtools/Solvers.F90 index c1cbdbacfe..07921ed889 100644 --- a/femtools/Solvers.F90 +++ b/femtools/Solvers.F90 @@ -26,6 +26,7 @@ ! USA #include "fdebug.h" module solvers + use iso_c_binding use FLDebug use Global_Parameters use futils, only: present_and_true, int2str, free_unit, real_format @@ -51,6 +52,17 @@ module solvers #include "petsc_legacy.h" + type, bind(C) :: solver_options_type + KSPType :: ksptype + PetscInt :: restart + PCType :: pctype + character(len=80) :: hypretype + PetscReal :: atol + PetscReal :: rtol + PetscInt :: max_its + PetscBool :: start_from_zero + end type solver_options_type + ! stuff used in the PETSc monitor (see petsc_solve_callback_setup() below) integer :: petsc_monitor_iteration = 0 Vec :: petsc_monitor_x @@ -79,7 +91,8 @@ module solvers private public petsc_solve, set_solver_options, & - complete_solver_option_path, petsc_solve_needs_positions + complete_solver_option_path, petsc_solve_needs_positions,& + solver_options_type, populate_solver_options_struct_from_path ! meant for unit-testing solver code only: public petsc_solve_setup, petsc_solve_core, petsc_solve_destroy, & @@ -2192,6 +2205,29 @@ subroutine set_solver_options_tensor(field, & end subroutine set_solver_options_tensor +subroutine populate_solver_options_struct_from_path(solver_options, option_path) + type(solver_options_type), intent(out) :: solver_options + character(len=*), intent(in):: option_path + + call get_option(trim(option_path)//'/iterative_method/name', solver_options%ksptype) + if (trim(solver_options%ksptype)=="gmres") & + call get_option(trim(option_path)//'/iterative_method/restart', solver_options%restart) + call get_option(trim(option_path)//'/preconditioner/name', solver_options%pctype) + if (trim(solver_options%pctype)=="hypre") & + call get_option(trim(option_path)//'/preconditioner/hypre_type/name', solver_options%hypretype) + call get_option(trim(option_path)//'/absolute_error', & + solver_options%atol, default=0.0) + call get_option(trim(option_path)//'/relative_error',solver_options%rtol) + call get_option(trim(option_path)//'/max_iterations',solver_options%max_its) + + if (have_option(trim(option_path)//'/start_from_zero')) then + solver_options%start_from_zero=PETSC_TRUE + else + solver_options%start_from_zero=PETSC_FALSE + endif + +end subroutine populate_solver_options_struct_from_path + subroutine petsc_monitor_setup(petsc_numbering, max_its) ! sets up the petsc monitors "exact" or "iteration_vtus" type(petsc_numbering_type), intent(in):: petsc_numbering diff --git a/include/solver_options.h b/include/solver_options.h new file mode 100644 index 0000000000..95bebbe945 --- /dev/null +++ b/include/solver_options.h @@ -0,0 +1,51 @@ +/* ! Copyright (C) 2006 Imperial College London and others. */ +/* ! */ +/* ! Please see the AUTHORS file in the main source directory for a full list */ +/* ! of copyright holders. */ +/* ! */ +/* ! Prof. C Pain */ +/* ! Applied Modelling and Computation Group */ +/* ! Department of Earth Science and Engineering */ +/* ! Imperial College London */ +/* ! */ +/* ! amcgsoftware@imperial.ac.uk */ +/* ! */ +/* ! This library is free software; you can redistribute it and/or */ +/* ! modify it under the terms of the GNU Lesser General Public */ +/* ! License as published by the Free Software Foundation, */ +/* ! version 2.1 of the License. */ +/* ! */ +/* ! This library is distributed in the hope that it will be useful, */ +/* ! but WITHOUT ANY WARRANTY; without even the implied warranty of */ +/* ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */ +/* ! Lesser General Public License for more details. */ +/* ! */ +/* ! You should have received a copy of the GNU Lesser General Public */ +/* ! License along with this library; if not, write to the Free Software */ +/* ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 */ +/* ! USA */ + +#ifndef SOLVER_OPTIONS_H +#define SOLVER_OPTIONS_H +#include + +struct solver_options { + char ksptype[80]; + PetscInt restart; + char pctype[80]; + char hypretype[80]; + PetscReal atol; + PetscReal rtol; + PetscInt max_its; + PetscBool start_from_zero; +}; + +void setup_ksp_from_options_c(KSP ksp, struct solver_options *options); +void petsc_solve_many_c(Vec *x, Mat M, Vec *b, int neqns, + struct solver_options *options, + int debug_level); +void petsc_solve_c(Vec x, Mat M, Vec b, + struct solver_options *options, + int debug_level); + +#endif diff --git a/main/Fluids.F90 b/main/Fluids.F90 index 4d113d9b4b..66229a7e42 100644 --- a/main/Fluids.F90 +++ b/main/Fluids.F90 @@ -342,6 +342,23 @@ SUBROUTINE FLUIDS() ! they will be updated (inside the call) call move_mesh_free_surface(state, initialise=.true.) + !! initialise laplacian smoothing if necessary + if (have_option("/mesh_adaptivity/mesh_movement/laplacian_smoothing")) then + call move_mesh_initialise_laplacian_smoothing(state) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/lineal_smoothing")) then + call move_mesh_initialise_lineal_smoothing(state) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/lineal_torsional_smoothing")) then + call move_mesh_initialise_lineal_torsional_smoothing(state) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/centroid_relaxer")) then + call move_mesh_initialise_centroid_relaxer(state) + end if + call run_diagnostics(state) !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC @@ -411,6 +428,22 @@ SUBROUTINE FLUIDS() call initialise_diagnostics(filename, state) + !! initialise laplacian smoothing if necessary + if (have_option("/mesh_adaptivity/mesh_movement/laplacian_smoothing")) then + call move_mesh_initialise_laplacian_smoothing(state) + !! call move_mesh_laplacian_smoothing(state, diagnostic_only=.true.) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/lineal_smoothing")) then + call move_mesh_initialise_lineal_smoothing(state) + !! call move_mesh_lineal_smoothing(state, diagnostic_only=.true.) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/lineal_torsional_smoothing")) then + call move_mesh_initialise_lineal_torsional_smoothing(state) + !! call move_mesh_lineal_torsional_smoothing(state, diagnostic_only=.true.) + end if + ! Initialise ice_meltrate, read constatns, allocate surface, and calculate melt rate if (have_option("/ocean_forcing/iceshelf_meltrate/Holland08")) then call melt_surf_init(state(1)) @@ -530,6 +563,10 @@ SUBROUTINE FLUIDS() ! to ensure those properties are satisfied on the new mesh not the old one. call move_mesh_imposed_velocity(state) call move_mesh_pseudo_lagrangian(state) + call move_mesh_laplacian_smoothing(state) + call move_mesh_lineal_smoothing(state) + call move_mesh_lineal_torsional_smoothing(state) + call move_mesh_centroid_relaxer(state) call enforce_discrete_properties(state, only_prescribed=.true., & exclude_interpolated=.true., & @@ -918,6 +955,27 @@ SUBROUTINE FLUIDS() if(have_option("/io/stat/output_after_adapts")) call write_diagnostics(state, current_time, dt, timestep, not_to_move_det_yet=.true.) call run_diagnostics(state) + !! initialise laplacian smoothing if necessary + if (have_option("/mesh_adaptivity/mesh_movement/laplacian_smoothing")) then + call move_mesh_initialise_laplacian_smoothing(state) + call move_mesh_laplacian_smoothing(state, diagnostic_only=.true.) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/lineal_smoothing")) then + call move_mesh_initialise_lineal_smoothing(state) + call move_mesh_lineal_smoothing(state, diagnostic_only=.true.) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/lineal_torsional_smoothing")) then + call move_mesh_initialise_lineal_torsional_smoothing(state) + call move_mesh_lineal_torsional_smoothing(state, diagnostic_only=.true.) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/centroid_relaxer")) then + call move_mesh_initialise_centroid_relaxer(state) + call move_mesh_centroid_relaxer(state, diagnostic_only=.true.) + end if + end if else if(have_option("/mesh_adaptivity/prescribed_adaptivity")) then if(do_adapt_state_prescribed(current_time)) then @@ -933,6 +991,27 @@ SUBROUTINE FLUIDS() if(have_option("/io/stat/output_after_adapts")) call write_diagnostics(state, current_time, dt, timestep, not_to_move_det_yet=.true.) call run_diagnostics(state) + !! initialise laplacian smoothing if necessary + if (have_option("/mesh_adaptivity/mesh_movement/laplacian_smoothing")) then + call move_mesh_initialise_laplacian_smoothing(state) + call move_mesh_laplacian_smoothing(state) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/lineal_smoothing")) then + call move_mesh_initialise_lineal_smoothing(state) + call move_mesh_lineal_smoothing(state) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/lineal_torsional_smoothing")) then + call move_mesh_initialise_lineal_torsional_smoothing(state) + call move_mesh_lineal_torsional_smoothing(state) + end if + + if (have_option("/mesh_adaptivity/mesh_movement/centroid_relaxer")) then + call move_mesh_initialise_centroid_relaxer(state) + call move_mesh_centroid_relaxer(state) + end if + end if not_to_move_det_yet=.false. diff --git a/main/Makefile.dependencies b/main/Makefile.dependencies index 4172d0e0cd..23fc0651a2 100644 --- a/main/Makefile.dependencies +++ b/main/Makefile.dependencies @@ -15,20 +15,20 @@ Fluids.o ../include/fluids_module.mod: Fluids.F90 \ ../include/diagnostic_children.mod ../include/diagnostic_fields_new.mod \ ../include/diagnostic_fields_wrapper.mod \ ../include/diagnostic_variables.mod \ - ../include/discrete_properties_module.mod ../include/elements.mod \ - ../include/equation_of_state.mod ../include/eventcounter.mod \ - ../include/fdebug.h ../include/field_equations_cv.mod \ - ../include/field_priority_lists.mod ../include/fields.mod \ - ../include/fldebug.mod ../include/foam_flow_module.mod \ - ../include/free_surface_module.mod ../include/futils.mod \ - ../include/global_parameters.mod ../include/gls.mod ../include/goals.mod \ - ../include/halos.mod ../include/iceshelf_meltrate_surf_normal.mod \ - ../include/implicit_solids.mod ../include/k_epsilon.mod \ - ../include/memory_diagnostics.mod ../include/meshdiagnostics.mod \ - ../include/meshmovement.mod ../include/momentum_diagnostic_fields.mod \ - ../include/momentum_equation.mod ../include/multimaterial_module.mod \ - ../include/multiphase_module.mod ../include/parallel_tools.mod \ - ../include/populate_state_module.mod \ + ../include/discrete_properties_module.mod ../include/dqmom.mod \ + ../include/elements.mod ../include/equation_of_state.mod \ + ../include/eventcounter.mod ../include/fdebug.h \ + ../include/field_equations_cv.mod ../include/field_priority_lists.mod \ + ../include/fields.mod ../include/fldebug.mod \ + ../include/foam_flow_module.mod ../include/free_surface_module.mod \ + ../include/futils.mod ../include/global_parameters.mod ../include/gls.mod \ + ../include/goals.mod ../include/halos.mod \ + ../include/iceshelf_meltrate_surf_normal.mod ../include/implicit_solids.mod \ + ../include/k_epsilon.mod ../include/memory_diagnostics.mod \ + ../include/meshdiagnostics.mod ../include/meshmovement.mod \ + ../include/momentum_diagnostic_fields.mod ../include/momentum_equation.mod \ + ../include/multimaterial_module.mod ../include/multiphase_module.mod \ + ../include/parallel_tools.mod ../include/populate_state_module.mod \ ../include/populate_sub_state_module.mod ../include/qmesh_module.mod \ ../include/reduced_model_runtime.mod ../include/reference_counting.mod \ ../include/reserve_state_module.mod \ diff --git a/manual/adaptivity.tex b/manual/adaptivity.tex index 1c6d012b39..c6d639cc25 100644 --- a/manual/adaptivity.tex +++ b/manual/adaptivity.tex @@ -619,4 +619,90 @@ \section{The cost of adaptivity} \label{table:adaptivity_cost} \end{table} +\section{Pure mesh movement} + +It will be seen from the previous section that the interpolation of data between meshes, and the use of metric advection to ensure the mesh remains of high quality between adapt steps are both relatively expensive processes. Nevertheless they can be necessitated by the fact that the $hr$--adaptivity algorithm described above inserts new degrees of freedom into the mesh and should not generally be run after every time step on grounds of cost. An alternative is to continuously modify the mesh over the course of the model time integration step, by moving the mesh vertices (i.e pure $r$--adaptivity) in the manner of the Arbitrary Lagrangian--Eulerian framework described by \cite{doneahuerta2004}. This gives a further advantage, in that the boundaries of the geometry can be explicitly allowed to move or deform. + +\begin{figure}[h] +\centering +\includegraphics[width=0.475\textwidth]{./adaptivity_images/moving_mesh_schematic} + +\caption[Example of a moving mesh.]{A cartoon of mesh movement for a fluids problem. While a Lagrangian fluid parcel is advected by the fluid velocity, $\bmu$, the mesh vertices are transported by a grid velocity, $\bm{v}$. In this frame of reference the material derivative is a function of of the convective velocity $\bm{c}:=\bm{u}-\bm{v}$, i.e. $\DDt{a}\equiv\frac{\partial a}{\partial{t}}+\bm{c}\cdot\nabla a$.} +\label{fig:mesh_movement} +\end{figure} + +In the \fluidity implementation of mesh movement, this method requires the specification of a continous, piecewise linear vector field, called the grid velocity, which defines the rate of change of the locations of the elements of the mesh. In this formalism the material derivative at fixed Eulerian location, +\begin{equation}\left.\DDt{a}\right|_{\bmx}\equiv\frac{\partial{a}}{\partial{t}}+\bmu\cdot\nabla a +\end{equation} is replaced by one moving with the mesh, +\begin{equation} +\left.\DDt{a}\right|_{\bm{\chi}}\equiv\frac{\partial{a}}{\partial{t}}+\bmu\cdot\nabla a-\bm{v}\cdot\nabla a +\end{equation} +where the mesh coordinate evolves under the action of the grid velocity, +\begin{equation} +\frac{\partial \bm{\chi}}{\partial t} = \bm{v}(\bm{\chi},t). +\end{equation} + This grid velocity may be explicitly chosen by the user. Alternatively, a number of algorithms are available to determine an optimal grid velocity from other information. + +\subsection{A pseudo-Lagrangian approach} + +Here the mesh velocity is set equal to the fluid velocity at the start of the timestep. This is only pseudo-Lagrangian, rather than fully Lagrangian since the two velocities diverge over the course of the time integration. This means that the advective term doesn't fully cancel. However stability criteria such as the CFL number will be derived from this convective velocity, which for steady flow will be smaller than the same criterion on a static mesh. + +%%\subsection{free surface method} + +%%\subsection{Explicit ALE methods} + +\subsection{A linear elasticity/viscosity scheme} + +In this approach \citep{alauzet2014changing} the grid velocity, $\bm{v}$ satisfies a partial differential equation based on a pseudo-force balance, with an optional drag term coupling to the fluid velocity, $\bmu$: +\[\frac{\partial \sigma(\bm{v})_{ij}}{\partial X_j} + C_1 u_j -C_2 v_j = 0,\] +with a grid pseudo stress tensor, +\[ \sigma_{ij} := \lambda \frac{\partial}{\partial X_i}\frac{\partial v_k}{\partial X_k}\delta_{ij} + 2 \mu \left(\frac{\partial v_i}{\partial X_j} + \frac{\partial v_i}{\partial X_j}\right).\] + +Here the user is free to choose the coordinate system, $\bm{X}$, in which the differential calculations are performed, as well as choosing the constand coefficients $\lambda$, $\mu$, $C_1$ and $C_2$. If the coordinate system of the original mesh is used (i.e. $\bm{X}=\bm{\chi}$, the mesh coordinate) then the stress equation can be identified with the equilibrium deformation of a linear elastic solid with Lamé parameters $\lambda$ and $\mu$. Equally if the coordinates of the current mesh are used (i.e. $\bm{X}=\bm{x}$, the \fluidity default) then the equation models a highly viscous compressible fluid, with shear viscosity $2\mu$ and bulk viscosity $\lambda +4\mu/3$. For well-posedness it is necessary that $\mu>0$ and $\lambda\geq0$. As implemented, the mesh movement force equation requires Dirichlet boundary conditions for $\bm{v}$ over the entire boundary of the solution domain. This information will then be smoothly extended into the interior of the domain. + +\subsection{Laplacian Smoothing} +In this approach \citep{buell1973,field} a boundary deformation displaces internal nodes according to the solution of a modified Laplace's equation +\begin{equation} +\label{eq:modified_laplaces} +\frac{{\partial}^{2}\bm{U}}{\partial \xi^{2}} + \frac{\partial^{2}\bm{U}}{\partial \eta^{2}} = 0, +\end{equation} +where $\left(\xi,\eta\right) \in \Omega_{c}\,(\textrm{predeformed domain}),\,\left(x,y\right) \in \Omega_{p}\,(\textrm{postdeformed domain}),\, \bm{U}\left(\xi,\eta\right)\,:\,\Omega_{c} \to \Omega_{p}$. Here it is important to highlight the underlying vector nature of $\bm{U}=\left(U^{\xi},U^{\eta}\right)$ meaning \ref{eq:modified_laplaces} can be recast as +\begin{equation} +\label{eq:recast_modified_laplaces} +\frac{\partial^{2}U^{\xi}}{\partial \xi^{2}} + \frac{\partial^{2}U^{\xi}}{\partial \eta^{2}} = 0\quad \textrm{and}\quad \frac{\partial^{2}U^{\eta}}{\partial \xi^{2}} + \frac{\partial^{2}U^{\eta}}{\partial \eta^{2}}=0, +\end{equation} + +where solutions to $U^{\xi}$ and $U^{\eta}$ can be found analytically by separation of variables. Having found $\bm{U}$ the grid velocity, $\bm{v}$, mentioned above, immediately follows. + +There are a number of important facts to keep in mind when using this smoothing technique. First, solutions to Laplace's equation in the context of this method, are not coupled (see \ref{eq:recast_modified_laplaces}). This means that a boundary deformation solely in the $y$-direction, for example, will cause nodes on the interior to be displaced solely along the $y$-direction. In a few instances this behaviour is desired (e.g. holes in domains), but more often than not this lack of coupling is a shortcoming of Laplacian smoothing. Second, Laplacian smoothing has no inbuilt safeguards against mesh tangling, which when deformations are large compared to the largest dimension of the domain, often occur near the boundary. Lastly, it is worth mentioning that for this particular smoothing technique, along with the following two, that in the absence of a boundary deformation, no mesh smoothing will occur along the interior. Laplacian smoothing, practically speaking, is computationally inexpensive and in general works well for convex domains that expand outward from their centres. However, if an external boundary moves inward too much and creates a nonconvex domain this smoother will most likely produce a tangled mesh. + +\subsection{Lineal Spring Analogy} +In this approach \citep{farhat1998} a mesh is viewed as a discrete structural system described as +\begin{equation} +\label{eq:discrete_struc_system} +\bm{M}\ddot{\bm{q}} + \bm{D}\dot{\bm{q}} + \bm{K}\bm{q} = 0, +\end{equation} +where $\bm{M}$, $\bm{D}$, $\bm{K}$, $\bm{q}$ are the fictitious mass, damping and stiffness matrices, respectively and $\bm{q}$ is the displacement vector. The displacement vector is defined as +\begin{equation} +\bm{q}\left(t\right) = \bm{\xi}\left(t\right) - \bm{\xi}\left(0\right), +\end{equation} +where $\bm{\xi}$ denotes the grid point coordinates and $t$ time. In this particular implementation, the quasi-static version of \ref{eq:discrete_struc_system} is assumed +\begin{equation} +\label{eq:static_version} +\bm{K}\bm{q}=0, +\end{equation} +and $\bm{q}=\bar{\bm{q}}$ along the moving boundary $\Gamma_{m}$ which is normally furnished by the solution to an underlying physical PDE. In the greater context of $r$-based adaptivity as it appears in this manual, $\bm{q}$ is the grid velocity $\bm{v}$. The main idea behind this technique is that a fictitious spring is attached to each edge connecting two vertices ($i$ and $j$), whose spring constant $k_{ij}$ is inversely proportional to the distance separating vertices $i$ and $j$. This means that $\bm{K} = \bm{K}^{ij}_{\textrm{lineal}}$. + +There are a number of important facts to keep in mind when using this smoothing technique. First, $k_{ij}$ is not a scalar, but is in fact a matrix which couples displacement in directions other than the direction of the boundary displacement (unlike Laplacian smoothing). Second, this technique only insures that neighbouring vertices never collide which in and of itself does not ensure the mesh will not tangle. Lastly, this technique operates on the coordinates of the current mesh and not the original mesh (see Lineal elasticity/viscosity scheme above). This smoother, practically speaking, can be thought of as a more robust version of Laplacian smoothing (without much more computational effort) and assuming the boundary deformation is small relative to the largest dimension of the domain, it should produce a reliable mesh. + +\subsection{Lineal Torsional Spring Analogy} +This approach \citep{farhat1998} is a direct extension of the Lineal Spring Analogy. What sets the two methods apart is that $\bm{K}$ as it appears in \ref{eq:static_version} is no longer solely a function of elemental edge length, but also elemental area. Conceptually this means $\bm{K} = \bm{K}^{ij}_{\textrm{lineal}} + \bm{K}^{ijk}_{\textrm{torsional}}$. Elemental (triangle $\mathcal{T}_{ijk}$ in this case) area contributes to torsional spring stiffness at vertex $i$ via +\begin{equation} +\label{eq:torsional_stiffness} +C^{ijk}_{i} = \frac{l^{2}_{ij}l^{2}_{ik}}{4A^{2}_{ijk}}, +\end{equation} +where $l_{ij}$ is the distance between vertices $i$ and $j$, $l_{ik}$ is the distance between vertices $i$ and $k$, and $A_{ijk}$ is the area of the triangle $\mathcal{T}_{ijk}$. This addition does not entirely safeguard against mesh tangling, and a further step must be taken to limit the rotation, $\bm{R}^{ijk}$, of each triangle within the mesh as well (see \citep{farhat1998} for further details). + +There are two important facts to keep in mind when using this smoothing technique. First, though this smoother ensures that no vertex will ever collide with its neighbour and that the mesh will not tangle (assuming it is given a valid mesh in the first place), this does guarantee good element quality (e.g. radii ratio). Second, because this smoother operates on the current mesh and not the original mesh, it can be regarded as a discrete version of the linear viscosity scheme discussed above. In turn, this technique may be an attractive alternative to those using the linear viscosity scheme. This smoother, practically speaking, offers very good performance (relative to Laplacian smoothing and Lineal spring smoothing) not just for external moving boundaries, but also for internal moving boundaries (e.g. turbines, pumps, etc.). A purely torsional implementation of this smoother for three-dimensional problems can be found in \citep{degand}. + diff --git a/manual/adaptivity_images/moving_mesh_schematic.pdf b/manual/adaptivity_images/moving_mesh_schematic.pdf new file mode 100644 index 0000000000..4034c652be Binary files /dev/null and b/manual/adaptivity_images/moving_mesh_schematic.pdf differ diff --git a/manual/bibliography.bib b/manual/bibliography.bib index 4daa961f3d..0a9b1ab6a2 100644 --- a/manual/bibliography.bib +++ b/manual/bibliography.bib @@ -485,6 +485,18 @@ @Article{ budd2009 bdsk-url-1 = {http://dx.doi.org/10.1017/S0962492906400015} } +@Article{ buell1973, + author = {Buell, W. R. and Bush, B. A.}, + journal = {Journal of Engineering for Industry}, + number = {1}, + pages = {332--338}, + title = {Mesh Generation - A Survey}, + volume = {95}, + year = {1973}, + doi = {10.1115/1.3438132}, + bdsk-url-1 = {http://dx.doi.org/10.1115/1.3438132} +} + @Article{ burchard1998, author = {Burchard, Hans and Petersen, Ole and Rippeth, Tom P.}, title = {Comparing the performance of the Mellor-Yamada and the @@ -847,7 +859,19 @@ @Article{ desprs_contact_2001 pages = {479--524} } -@Article{ devine2002, +@Article{ degand, + author = {Degand, C. and Farhat, C.}, + title = {A three-dimensional torsional spring analogy method for unstructured dynamic meshes}, + volume = {80}, + number = {3--4}, + journal = {Computers and Structures}, + year = {2002}, + pages = {305--316}, + doi = {10.1016/S0045-7949(02)00002-0}, + bdsk-url-1 = {http://dx.doi.org/10.1016/S0045-7949(02)00002-0} +} + +@Article{ devine2002, author = {Karen Devine and Erik Boman and Robert Heaphy and Bruce Hendrickson and Courtenay Vaughan}, title = {{Zoltan} Data Management Services for Parallel Dynamic @@ -965,7 +989,17 @@ @Book{ fannelop_94 year = {1994} } -@TechReport{ farnell1975, +@Article{ farhat1998, + author = {Farhat, C. and Degand, B. and Koobus, B. and Lesoinne, M.}, + title = {Torsional springs for two-dimensional dynamic unstructured fluid meshes}, + journal = {Computer Methods in Applied Mechanics and Engineering}, + volume = {163}, + pages = {231--245}, + year = {1998}, + doi = {10.1016/S0045-7825(98)00016-4}, + bdsk-url-1 = {http://dx.doi.org/10.1016/S0045-7825(98)00016-4} +} +@TechReportq{ farnell1975, author = {Farnell, L. and Plumb, R. A.}, title = {{Numerical integration of flow in a rotating annulus I: Axisymmetric model}}, @@ -1026,6 +1060,19 @@ @Manual{ femtools year = 2009 } + +@Article{ field, + author = {Field, D. A.}, + journal = {Communications in Applied Numerical Methods}, + number = {6}, + pages = {709--712}, + title = {Laplacian smoothing and Delaunay triangulations}, + volume = {4}, + year = {1988}, + doi = {10.1002/cnm.1630040603}, + bdsk-url-1 = {http://dx.doi.org/10.1002/cnm.1630040603} +} + @Article{ fordii, author = {R. Ford and C. C. Pain and M. D. Piggott and A. J. H. Goddard and C. R. E. de Oliveira and A. P. Umbleby}, @@ -3185,4 +3232,25 @@ @article{mccoymadras2003analytical pages={3049--3051}, year={2003}, publisher={Pergamon} +} + +@incollection{doneahuerta2004, + title = {Arbitrary Lagrangian Eulerian methods}, + author = {Donea, Jean and Huerta, Antonio and Ponthot, Jean-Philippe and others}, + booktitle = {Encyclopedia of Computational Mechanics}, + chapter = 14, + year ={2004}, + editor = {E. Stein, R. de Borst and TJR Hughes, + publisher= {Wiley \& sons} +} + +@article{alauzet2014changing, + title={A changing-topology moving mesh technique for large displacements}, + author={Alauzet, Fr{\'e}d{\'e}ric}, + journal={Engineering with Computers}, + volume={30}, + number={2}, + pages={175--200}, + year={2014}, + publisher={Springer} } \ No newline at end of file diff --git a/population_balance/Makefile.dependencies b/population_balance/Makefile.dependencies index a4f50cde3b..770c52b14c 100644 --- a/population_balance/Makefile.dependencies +++ b/population_balance/Makefile.dependencies @@ -2,9 +2,21 @@ ../include/dqmom.mod: DQMOM.o @true -DQMOM.o ../include/dqmom.mod: DQMOM.F90 ../include/fdebug.h \ - ../include/fields.mod ../include/global_parameters.mod \ - ../include/initialise_fields_module.mod ../include/spud.mod \ +DQMOM.o ../include/dqmom.mod: DQMOM.F90 ../include/elements.mod \ + ../include/fdebug.h ../include/fetools.mod ../include/fields.mod \ + ../include/fldebug.mod ../include/futils.mod \ + ../include/global_parameters.mod ../include/initialise_fields_module.mod \ + ../include/solvers.mod ../include/sparse_tools.mod \ ../include/state_fields_module.mod ../include/state_module.mod \ ../include/vector_tools.mod +../include/dqmom.mod: DQMOM_node_source.o + @true + +DQMOM_node_source.o ../include/dqmom.mod: DQMOM_node_source.F90 \ + ../include/fdebug.h ../include/fields.mod ../include/fldebug.mod \ + ../include/futils.mod ../include/global_parameters.mod \ + ../include/initialise_fields_module.mod ../include/solvers.mod \ + ../include/sparse_tools.mod ../include/state_fields_module.mod \ + ../include/state_module.mod ../include/vector_tools.mod + diff --git a/preprocessor/Boundary_Conditions_From_Options.F90 b/preprocessor/Boundary_Conditions_From_Options.F90 index 9cacb3cc61..5fe8d779a1 100644 --- a/preprocessor/Boundary_Conditions_From_Options.F90 +++ b/preprocessor/Boundary_Conditions_From_Options.F90 @@ -151,6 +151,9 @@ subroutine populate_boundary_conditions(states, suppress_warnings) vfield => extract_vector_field(states(p+1), f) field_path=vfield%option_path + call populate_vector_boundary_conditions(states(p+1),vfield, & + trim(field_path)//'/diagnostic/algorithm/boundary_conditions', position, suppress_warnings=suppress_warnings) + if (.not. have_option(trim(field_path)//'/prognostic')) cycle ! only prognostic fields from here: @@ -390,7 +393,7 @@ subroutine populate_vector_boundary_conditions(state, field, bc_path, position, end if select case(trim(bc_type)) - case("dirichlet", "neumann", "weakdirichlet", "flux") + case("dirichlet", "neumann", "weakdirichlet", "flux", "solid_body", "weaksolid_body") if(have_option(trim(bc_path_i)//"/type[0]/align_bc_with_cartesian")) then aligned_components=cartesian_aligned_components @@ -674,6 +677,13 @@ subroutine set_boundary_conditions_values(states, shift_time) do f = 1, nfields vfield => extract_vector_field(states(p+1), f) field_path=vfield%option_path + + + call set_vector_boundary_conditions_values(states(p+1), vfield, & + trim(field_path)//'/diagnostic/algorithm/boundary_conditions', & + position, shift_time=shift_time) + + if (.not. have_option(trim(field_path)//'/prognostic')) cycle ! only prognostic fields from here: diff --git a/schemas/adaptivity_options.rnc b/schemas/adaptivity_options.rnc index 2f426d9dae..edd5f59969 100644 --- a/schemas/adaptivity_options.rnc +++ b/schemas/adaptivity_options.rnc @@ -1864,6 +1864,32 @@ mesh_adaptivity_options = element imposed_grid_velocity { empty }| + ## Enable movement of mesh to satisfy a moving mesh PDE. + ## The boundary condition on the GridVelocity field is imposed + ## on the calculated solution. + element laplacian_smoothing { + empty + }| + ## Enable Lineal Spring. + ## The boundary condition on the GridVelocity field is imposed + ## on the calculated solution. + element lineal_smoothing { + empty + }| + ## Enable Lineal Torsional Spring. + ## The boundary condition on the GridVelocity field is imposed + ## on the calculated solution. + element lineal_torsional_smoothing { + empty + }| + ## Enable Centroid Relaxer. This relaxer is an adaptation of Johnson's + ## "Using Parameterization and Springs to Determine Aneurysm Wall Thickess". + ## DOI: 10.1007/978-3-642-04319-2_23. + ## The boundary condition on the GridVelocity field is imposed + ## on the calculated solution. + element centroid_relaxer { + empty + }| ## Enable movement of mesh by the Velocity. ## Note that this is only pseudo-lagrangian as the ## Velocity and GridVelocity only match at the beginning of @@ -1933,13 +1959,35 @@ mesh_adaptivity_options = ) } ), + ## Specify the solver options for methods which solve a linear equation using an internal algorithm. Default optionds are used otherwise. + element solver { + pc_ksp_solver_options + }?, ## The velocity of the mesh. element vector_field { attribute name { "GridVelocity" }, attribute rank { "1" }, ( - element diagnostic { - internal_algorithm, + element diagnostic { + ( + ## This diagnostic is internal - i.e. it is calculated + ## somewhere within the main code, and is not wrapped by + ## the automatic diagnostics wrappers. + element algorithm { + attribute name { "Internal" }, + attribute material_phase_support { "multiple" }, + element boundary_conditions { + attribute name { string }, + ## Surface id: + element surface_ids { + integer_vector + }, + velocity_boundary_conditions + }* + } + | + stress_free_field_algorithm + ), element mesh { attribute name { "CoordinateMesh" } }, diff --git a/schemas/adaptivity_options.rng b/schemas/adaptivity_options.rng index 260df54297..54d804da0f 100644 --- a/schemas/adaptivity_options.rng +++ b/schemas/adaptivity_options.rng @@ -2205,6 +2205,32 @@ on the cost of a worse internal representation of the wetting and drying interfa Requires a prescribed GridVelocity field (see below). + + Enable movement of mesh to satisfy a moving mesh PDE. +The boundary condition on the GridVelocity field is imposed +on the calculated solution. + + + + Enable Lineal Spring. +The boundary condition on the GridVelocity field is imposed +on the calculated solution. + + + + Enable Lineal Torsional Spring. +The boundary condition on the GridVelocity field is imposed +on the calculated solution. + + + + Enable Centroid Relaxer. This relaxer is an adaptation of Johnson's +"Using Parameterization and Springs to Determine Aneurysm Wall Thickess". +DOI: 10.1007/978-3-642-04319-2_23. +The boundary condition on the GridVelocity field is imposed +on the calculated solution. + + Enable movement of mesh by the Velocity. Note that this is only pseudo-lagrangian as the @@ -2342,6 +2368,12 @@ Please specify a scalar field + + + Specify the solver options for methods which solve a linear equation using an internal algorithm. Default optionds are used otherwise. + + + The velocity of the mesh. @@ -2352,7 +2384,32 @@ Please specify a scalar field - + + + This diagnostic is internal - i.e. it is calculated +somewhere within the main code, and is not wrapped by +the automatic diagnostics wrappers. + + Internal + + + multiple + + + + + + + + Surface id: + + + + + + + + CoordinateMesh diff --git a/schemas/diagnostic_algorithms.rnc b/schemas/diagnostic_algorithms.rnc index ef24bc49e0..6797f17b6e 100644 --- a/schemas/diagnostic_algorithms.rnc +++ b/schemas/diagnostic_algorithms.rnc @@ -221,6 +221,43 @@ tensor_copy_algorithm = }? } ) +stress_free_field_algorithm = + ( + ## Finds the vector field which solves the equation + ## + ## div sigma (v) - C1*v +C2*u_s = 0 + ## + ## where sigma is a stress-like tensor operator, + ## + ## sigma(v) = lambda div v I + 2 * mu * (nabla v + nabla v^T). + + element algorithm { + attribute name { "stress_free_field" }, + attribute material_phase_support { "single" }, + element solver { + linear_solver_options_sym + }, + element coordinate_field { + attribute name { string } + }?, + element lambda_parameter { real }?, + element mu_parameter { real }?, + element C1 { real }?, + element C2 { real }?, + vector_source_field, + element chi_parameter { real }?, + element depends { string }?, + element boundary_conditions { + attribute name { string }, + ## Surface id: + element surface_ids { + integer_vector + }, + velocity_boundary_conditions + }* + } + ) + scalar_galerkin_projection_algorithm = ( ## Galerkin projects the source field. @@ -1597,6 +1634,7 @@ vector_diagnostic_algorithms |= vector_galerkin_projection_algorithm vector_diagnostic_algorithms |= vector_difference_algorithm vector_diagnostic_algorithms |= vector_copy_algorithm vector_diagnostic_algorithms |= vector_advection_algorithm +vector_diagnostic_algorithms |= stress_free_field_algorithm vector_diagnostic_algorithms |= perp_algorithm vector_diagnostic_algorithms |= lumped_mass_smoothed_vector_algorithm vector_diagnostic_algorithms |= helmholtz_smoothed_vector_algorithm diff --git a/schemas/diagnostic_algorithms.rng b/schemas/diagnostic_algorithms.rng index c4793dde19..2f35d7deb8 100644 --- a/schemas/diagnostic_algorithms.rng +++ b/schemas/diagnostic_algorithms.rng @@ -310,6 +310,76 @@ this option to these these internal algorithms. + + + Finds the vector field which solves the equation + +div sigma (v) - C1*v +C2*u_s = 0 + +where sigma is a stress-like tensor operator, + +sigma(v) = lambda div v I + 2 * mu * (nabla v + nabla v^T). + + stress_free_field + + + single + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Surface id: + + + + + + + Galerkin projects the source field. @@ -2249,6 +2319,9 @@ a horizontal vector living on the specified surface. + + + diff --git a/schemas/fluidity_options.rnc b/schemas/fluidity_options.rnc index 893f13ea61..8198caf218 100644 --- a/schemas/fluidity_options.rnc +++ b/schemas/fluidity_options.rnc @@ -6153,6 +6153,20 @@ vector_field_choice = diagnostic_vector_field } ) + }| + ## Store the initial mesh coordinates, for solids/mesh movement + element vector_field { + attribute rank { "1" }, + attribute name { "InitialCoordinate" }, + ( + element diagnostic { + internal_algorithm, + element mesh { + attribute name {"CoordinateMesh"} + }, + diagnostic_vector_field + } + ) } # Insert new diagnostic vector field here using the template: diff --git a/schemas/fluidity_options.rng b/schemas/fluidity_options.rng index 3b7487a698..777f02eb86 100644 --- a/schemas/fluidity_options.rng +++ b/schemas/fluidity_options.rng @@ -7912,6 +7912,24 @@ algorithm option). + + Store the initial mesh coordinates, for solids/mesh movement + + 1 + + + InitialCoordinate + + + + + + CoordinateMesh + + + + +