← Back to branch summary

~f-milthaler/fluidity/fsi-model-stationary-solid-with-velocity

« back to all changes in this revision

Viewing changes to assemble/Free_Surface.F90

  • Committer: f.milthaler10 at uk
  • Date: 2013-11-06 13:43:56 UTC
  • mfrom: (3463.184.85 fluidity)
  • Revision ID: f.milthaler10@imperial.ac.ic.uk.-20131106134356-v3lw1dheesckywj0
merge from trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
assemble/Free_Surface.F90
297
297
      end if
298
298
      
299
299
      ! reference density
300
 
      call get_reference_density_from_options(rho0, state%option_path)
 
300
      call get_fs_reference_density_from_options(rho0, state%option_path)
301
301
      density => extract_scalar_field(state, "Density", stat=dens_stat)
302
302
      have_density = (dens_stat==0)
303
303
 
671
671
    gravity_normal => extract_vector_field(state, "GravityDirection")
672
672
    density => extract_scalar_field(state, "Density", stat=stat)
673
673
    have_density = (stat==0)
674
 
    call get_reference_density_from_options(rho0, state%option_path)
 
674
    call get_fs_reference_density_from_options(rho0, state%option_path)
675
675
    call get_option('/timestepping/timestep', dt)
676
676
 
677
677
    surface_sparsity => get_csr_sparsity_firstorder(state, surface_mesh, surface_mesh)
852
852
    end if
853
853
    
854
854
    ! reference density
855
 
    call get_reference_density_from_options(rho0, state%option_path)
 
855
    call get_fs_reference_density_from_options(rho0, state%option_path)
856
856
    density => extract_scalar_field(state, "Density", stat=dens_stat)
857
857
    have_density = (dens_stat==0)
858
858
 
1113
1113
    end if
1114
1114
    density => extract_scalar_field(state, "Density", stat=stat)
1115
1115
    have_density = (stat==0)
1116
 
    call get_reference_density_from_options(rho0, state%option_path)
 
1116
    call get_fs_reference_density_from_options(rho0, state%option_path)
1117
1117
 
1118
1118
    surface_sparsity => get_csr_sparsity_firstorder(state, scaled_fs%mesh, scaled_fs%mesh)
1119
1119
    call allocate(fs_matrix, surface_sparsity, name="FSMatrix")
1324
1324
      call allocate(extended_mesh, node_count(pressure_mesh)+node_count(fs_mesh), &
1325
1325
          element_count(pressure_mesh), pressure_mesh%shape, &
1326
1326
          "Extended"//trim(pressure_mesh%name))
 
1327
      extended_mesh%periodic = pressure_mesh%periodic
1327
1328
 
1328
1329
      do ele=1, element_count(pressure_mesh)
1329
1330
        nodes => ele_nodes(pressure_mesh,ele)
1795
1796
   move_mesh = have_option("/mesh_adaptivity/mesh_movement/free_surface")
1796
1797
   x => extract_vector_field(state, "Coordinate")
1797
1798
   ! reference density
1798
 
   call get_reference_density_from_options(rho0, state%option_path)
 
1799
   call get_fs_reference_density_from_options(rho0, state%option_path)
1799
1800
   density => extract_scalar_field(state, "Density", stat=dens_stat)
1800
1801
   have_density = (dens_stat==0)
1801
1802
   call get_option('/physical_parameters/gravity/magnitude', g, stat=stat)
2378
2379
     have_wd=have_option("/mesh_adaptivity/mesh_movement/free_surface/wetting_and_drying")
2379
2380
 
2380
2381
     u => extract_vector_field(state, "Velocity")
2381
 
     call get_reference_density_from_options(rho0, state%option_path)
 
2382
     call get_fs_reference_density_from_options(rho0, state%option_path)
 
2383
 
 
2384
     if (have_option(trim(u%option_path)//'/prognostic/equation::ShallowWater')) then
 
2385
       ! For the shallow water equations we only have a 2D, horizontal mesh
 
2386
       ! and f.s. is simply p/g everywhere:
 
2387
       call set(free_surface, p)
 
2388
       call scale(free_surface, 1./g)
 
2389
       return
 
2390
     end if
2382
2391
 
2383
2392
     !
2384
2393
     ! first we compute the right free surface values at the free surface
2490
2499
       if (bctype=="free_surface" .and. has_scalar_surface_field(u, i, "WettingDryingAlpha")) then
2491
2500
             scalar_surface_field => extract_scalar_surface_field(u, i, "WettingDryingAlpha")
2492
2501
             ! Update WettingDryingAlpha
2493
 
             call get_reference_density_from_options(rho0, state%option_path)
 
2502
             call get_fs_reference_density_from_options(rho0, state%option_path)
2494
2503
             original_bottomdist_remap => extract_scalar_field(state, "OriginalDistanceToBottomPressureMesh") 
2495
2504
             call get_option('/physical_parameters/gravity/magnitude', g)
2496
2505
             call get_option("/mesh_adaptivity/mesh_movement/free_surface/wetting_and_drying/d0", d0)
2605
2614
      end if
2606
2615
      
2607
2616
      ! reference density
2608
 
      call get_reference_density_from_options(rho0, state%option_path)
 
2617
      call get_fs_reference_density_from_options(rho0, state%option_path)
2609
2618
      ! Timestep
2610
2619
      call get_option('/timestepping/timestep', dt)
2611
2620
      move_mesh = have_option("/mesh_adaptivity/mesh_movement/free_surface")
2691
2700
    
2692
2701
    character(len=OPTION_PATH_LEN):: option_path, phase_path, pressure_path, pade_path
2693
2702
    character(len=FIELD_NAME_LEN):: fs_meshname, p_meshname, bctype
2694
 
    logical:: have_free_surface, have_explicit_free_surface, have_viscous_free_surface, have_wd
 
2703
    logical:: have_free_surface, have_explicit_free_surface, have_viscous_free_surface
 
2704
    logical:: have_wd, have_swe
2695
2705
    integer i, p
2696
2706
 
2697
2707
    do p=1, option_count('/material_phase')
2727
2737
        have_viscous_free_surface = .false.
2728
2738
      end if
2729
2739
      have_wd=have_option("/mesh_adaptivity/mesh_movement/free_surface/wetting_and_drying")
 
2740
 
 
2741
      have_swe=have_option(trim(option_path)//'/equation::ShallowWater')
2730
2742
      
2731
2743
      if (have_free_surface) then
2732
2744
         ewrite(2,*) "You have a free surface boundary condition, checking its options"
2759
2771
      if (have_option(trim(option_path))) then
2760
2772
        call get_option(trim(option_path)//'/mesh[0]/name', fs_meshname)
2761
2773
        call get_option(trim(pressure_path)//'/mesh[0]/name', p_meshname)
2762
 
        if (.not. have_free_surface) then
 
2774
        if (.not. (have_free_surface .or. have_swe)) then
2763
2775
          ewrite(-1,*) "The diagnostic FreeSurface field has to be used in combination " // &
2764
 
            "with the free_surface boundary condition under Velocity."
 
2776
            "with the free_surface boundary condition under Velocity, or with " // &
 
2777
            "equation type ShallowWater for Velocity."
2765
2778
          FLExit("Exit")
2766
2779
        end if
2767
2780
        if (.not. fs_meshname==p_meshname) then
2768
2781
          FLExit("The diagnostic FreeSurface field and the Pressure field have to be on the same mesh")
2769
2782
        end if
2770
 
        if (.not. have_option('/geometry/ocean_boundaries')) then
 
2783
        if (.not. have_option('/geometry/ocean_boundaries') .and. .not. have_swe) then
2771
2784
          ewrite(0,*) "Warning: your diagnostic free surface will only be " // &
2772
2785
            "defined at the free surface nodes and not extrapolated downwards, " // &
2773
2786
            "because you didn't specify geometry/ocean_boundaries."