← Back to branch summary

~fluidity-core/fluidity/excise-fldecomp

« back to all changes in this revision

Viewing changes to assemble/Advection_Diffusion_CG.F90

  • Committer: Mark Filipiak
  • Date: 2012-08-13 11:42:30 UTC
  • mfrom: (4003.1.23 dev-trunk)
  • Revision ID: mjf@staffmail.ed.ac.uk-20120813114230-wzoyf2gi4p4oxeh4
Merge in of the latest trunk.  To try to cure non-flredecomp tests that are passing at EPCC but failing in buildbot.

Show diffs side-by-side

added added

removed removed

Lines of Context:
assemble/Advection_Diffusion_CG.F90
47
47
  use upwind_stabilisation
48
48
  use sparsity_patterns_meshes
49
49
  use porous_media
 
50
  use multiphase_module
50
51
  use sparse_tools_petsc
51
52
  use colouring
52
53
#ifdef _OPENMP
117
118
  logical :: move_mesh
118
119
  ! Include porosity?
119
120
  logical :: include_porosity
 
121
  ! Are we running a multiphase flow simulation?
 
122
  logical :: multiphase
120
123
 
121
124
contains
122
125
 
123
 
  subroutine solve_field_equation_cg(field_name, state, dt, velocity_name, iterations_taken)
 
126
  subroutine solve_field_equation_cg(field_name, state, istate, dt, velocity_name, iterations_taken)
124
127
    !!< Construct and solve the advection-diffusion equation for the given
125
128
    !!< field using a continuous Galerkin discretisation. Based on
126
129
    !!< Advection_Diffusion_DG and Momentum_CG.
127
130
    
128
131
    character(len = *), intent(in) :: field_name
129
 
    type(state_type), intent(inout) :: state
 
132
    type(state_type), dimension(:), intent(inout) :: state
 
133
    integer, intent(in) :: istate
130
134
    real, intent(in) :: dt
131
135
    character(len = *), optional, intent(in) :: velocity_name
132
136
    integer, intent(out), optional :: iterations_taken
138
142
    ewrite(1, *) "In solve_field_equation_cg"
139
143
    
140
144
    ewrite(2, *) "Solving advection-diffusion equation for field " // &
141
 
      & trim(field_name) // " in state " // trim(state%name)
 
145
      & trim(field_name) // " in state " // trim(state(istate)%name)
142
146
 
143
 
    call initialise_advection_diffusion_cg(field_name, t, delta_t, matrix, rhs, state)
 
147
    call initialise_advection_diffusion_cg(field_name, t, delta_t, matrix, rhs, state(istate))
144
148
    
145
149
    call profiler_tic(t, "assembly")
146
 
    call assemble_advection_diffusion_cg(t, matrix, rhs, state, dt, velocity_name = velocity_name)    
 
150
    call assemble_advection_diffusion_cg(t, matrix, rhs, state(istate), dt, velocity_name = velocity_name)    
 
151
    
 
152
    ! Note: the assembly of the heat transfer term is done here to avoid 
 
153
    ! passing in the whole state array to assemble_advection_diffusion_cg.
 
154
    if(have_option("/multiphase_interaction/heat_transfer")) then
 
155
       call add_heat_transfer(state, istate, t, matrix, rhs)
 
156
    end if
147
157
    call profiler_toc(t, "assembly")
148
 
 
 
158
    
149
159
    call profiler_tic(t, "solve_total")
150
 
    call solve_advection_diffusion_cg(t, delta_t, matrix, rhs, state, &
 
160
    call solve_advection_diffusion_cg(t, delta_t, matrix, rhs, state(istate), &
151
161
                                      iterations_taken = iterations_taken)
152
162
    call profiler_toc(t, "solve_total")
153
163
 
234
244
    type(scalar_field), pointer :: density, olddensity
235
245
    character(len = FIELD_NAME_LEN) :: density_name
236
246
    type(scalar_field), pointer :: pressure
 
247
    
 
248
    ! Volume fraction fields for multiphase flow simulation
 
249
    type(scalar_field), pointer :: vfrac
 
250
    type(scalar_field) :: nvfrac ! Non-linear version
 
251
 
237
252
    ! Porosity field
238
253
    type(scalar_field) :: porosity_theta
239
254
        
373
388
       call allocate(porosity_theta, t%mesh, field_type=FIELD_TYPE_CONSTANT)
374
389
       call set(porosity_theta, 1.0)
375
390
    end if
 
391
 
376
392
    
377
393
    ! Step 2: Pull options out of the options tree
378
394
    
462
478
      stabilisation_scheme = STABILISATION_NONE
463
479
    end if
464
480
    
 
481
    ! PhaseVolumeFraction for multiphase flow simulations
 
482
    if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then
 
483
       multiphase = .true.
 
484
       vfrac => extract_scalar_field(state, "PhaseVolumeFraction")
 
485
       call allocate(nvfrac, vfrac%mesh, "NonlinearPhaseVolumeFraction")
 
486
       call zero(nvfrac)
 
487
       call get_nonlinear_volume_fraction(state, nvfrac)
 
488
      
 
489
       ewrite_minmax(nvfrac)
 
490
    else
 
491
       multiphase = .false.
 
492
       call allocate(nvfrac, t%mesh, "DummyNonlinearPhaseVolumeFraction", field_type=FIELD_TYPE_CONSTANT)
 
493
       call set(nvfrac, 1.0)
 
494
    end if
 
495
      
465
496
    call allocate(dummydensity, t%mesh, "DummyDensity", field_type=FIELD_TYPE_CONSTANT)
466
497
    call set(dummydensity, 1.0)
467
498
    ! find out equation type and hence if density is needed or not
479
510
      if(move_mesh) then
480
511
        FLExit("Haven't implemented a moving mesh energy equation yet.")
481
512
      end if
 
513
      
 
514
      ! Get old and current densities
482
515
      call get_option(trim(t%option_path)//'/prognostic/equation[0]/density[0]/name', &
483
516
                      density_name)
484
517
      density=>extract_scalar_field(state, trim(density_name))
485
518
      ewrite_minmax(density)
486
 
      
487
519
      olddensity=>extract_scalar_field(state, "Old"//trim(density_name))
488
520
      ewrite_minmax(olddensity)
489
521
      
490
 
      call get_option(trim(density%option_path)//"/prognostic/temporal_discretisation/theta", &
491
 
                      density_theta)
 
522
      if(.not.multiphase .or. have_option('/material_phase::'//trim(state%name)//&
 
523
         &'/equation_of_state/compressible')) then         
 
524
         call get_option(trim(density%option_path)//"/prognostic/temporal_discretisation/theta", &
 
525
                        density_theta)
 
526
      else
 
527
         ! Since the particle phase is always incompressible then its Density
 
528
         ! will not be prognostic. Just use a fixed theta value of 1.0.
 
529
         density_theta = 1.0
 
530
      end if
492
531
                      
493
532
      pressure=>extract_scalar_field(state, "Pressure")
494
533
      ewrite_minmax(pressure)
 
534
 
495
535
    case default
496
536
      FLExit("Unknown field equation type for cg advection diffusion.")
497
537
    end select
531
571
              positions, old_positions, new_positions, &
532
572
              velocity, grid_velocity, &
533
573
              source, absorption, diffusivity, &
534
 
              density, olddensity, pressure, porosity_theta, &
 
574
              density, olddensity, pressure, porosity_theta, nvfrac, &
535
575
              supg_element(thread_num+1))
536
576
      end do element_loop
537
577
      !$OMP END DO
545
585
    ! which must be included before dirichlet BC's.
546
586
    if (add_src_directly_to_rhs) call addto(rhs, source)
547
587
    
 
588
    
548
589
    ! Step 4: Boundary conditions
549
590
    
550
591
    if( &
571
612
        call assemble_advection_diffusion_face_cg(i, t_bc_types(i), t, t_bc, t_bc_2, &
572
613
                                                  matrix, rhs, &
573
614
                                                  positions, velocity, grid_velocity, &
574
 
                                                  density, olddensity)
 
615
                                                  density, olddensity, nvfrac)
575
616
      end do
576
617
    
577
618
      call deallocate(t_bc)
586
627
    ewrite_minmax(rhs)
587
628
    
588
629
    call deallocate(velocity)
 
630
    call deallocate(nvfrac)
589
631
    call deallocate(dummydensity)
590
632
    if (stabilisation_scheme == STABILISATION_SUPG) then
591
633
       do i = 1, num_threads
645
687
                                      positions, old_positions, new_positions, &
646
688
                                      velocity, grid_velocity, &
647
689
                                      source, absorption, diffusivity, &
648
 
                                      density, olddensity, pressure, porosity_theta, supg_shape)
 
690
                                      density, olddensity, pressure, porosity_theta, nvfrac, supg_shape)
649
691
    integer, intent(in) :: ele
650
692
    type(scalar_field), intent(in) :: t
651
693
    type(csr_matrix), intent(inout) :: matrix
661
703
    type(scalar_field), intent(in) :: olddensity
662
704
    type(scalar_field), intent(in) :: pressure
663
705
    type(scalar_field), intent(in) :: porosity_theta
 
706
    type(scalar_field), intent(in) :: nvfrac
664
707
    type(element_type), intent(inout) :: supg_shape
665
708
 
666
709
    integer, dimension(:), pointer :: element_nodes
669
712
    real, dimension(ele_loc(density, ele), ele_ngi(density, ele), mesh_dim(density)) :: drho_t
670
713
    real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: du_t 
671
714
    real, dimension(ele_loc(positions, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: dug_t 
 
715
    ! Derivative of shape function for nvfrac field
 
716
    real, dimension(ele_loc(nvfrac, ele), ele_ngi(nvfrac, ele), mesh_dim(nvfrac)) :: dnvfrac_t
 
717
    
672
718
    real, dimension(mesh_dim(t), mesh_dim(t), ele_ngi(t, ele)) :: j_mat 
673
719
    type(element_type) :: test_function
674
720
    type(element_type), pointer :: t_shape
741
787
      end if
742
788
    end if
743
789
    
 
790
    if(have_advection .and. multiphase .and. (equation_type==FIELD_EQUATION_INTERNALENERGY)) then
 
791
      ! If the field and nvfrac meshes are different, then we need to
 
792
      ! compute the derivatives of the nvfrac shape functions.    
 
793
      if(ele_shape(nvfrac, ele) == t_shape) then
 
794
         dnvfrac_t = dt_t
 
795
      else
 
796
         call transform_to_physical(positions, ele, ele_shape(nvfrac, ele), dshape=dnvfrac_t)
 
797
      end if
 
798
    end if
 
799
 
744
800
    ! Step 2: Set up test function
745
801
    
746
802
    select case(stabilisation_scheme)
763
819
    ! Step 3: Assemble contributions
764
820
    
765
821
    ! Mass
766
 
    if(have_mass) call add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto)
 
822
    if(have_mass) call add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, nvfrac, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto)
767
823
    
768
824
    ! Advection
769
825
    if(have_advection) call add_advection_element_cg(ele, test_function, t, &
770
826
                                        velocity, grid_velocity, diffusivity, &
771
 
                                        density, olddensity, &
772
 
                                        dt_t, du_t, dug_t, drho_t, detwei, j_mat, matrix_addto, rhs_addto)
 
827
                                        density, olddensity, nvfrac, &
 
828
                                        dt_t, du_t, dug_t, drho_t, dnvfrac_t, detwei, j_mat, matrix_addto, rhs_addto)
773
829
        
774
830
    ! Absorption
775
831
    if(have_absorption) call add_absorption_element_cg(ele, test_function, t, absorption, detwei, matrix_addto, rhs_addto)
784
840
    
785
841
    ! Pressure
786
842
    if(equation_type==FIELD_EQUATION_INTERNALENERGY) call add_pressurediv_element_cg(ele, test_function, t, &
787
 
                                                                                  velocity, pressure, &
788
 
                                                                                  du_t, detwei, rhs_addto)
 
843
                                                                                  velocity, pressure, nvfrac, &
 
844
                                                                                  du_t, dnvfrac_t, detwei, rhs_addto)
789
845
    
790
846
    ! Step 4: Insertion
791
847
            
795
851
 
796
852
  end subroutine assemble_advection_diffusion_element_cg
797
853
  
798
 
  subroutine add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto)
 
854
  subroutine add_mass_element_cg(ele, test_function, t, density, olddensity, porosity_theta, nvfrac, detwei, detwei_old, detwei_new, matrix_addto, rhs_addto)
799
855
    integer, intent(in) :: ele
800
856
    type(element_type), intent(in) :: test_function
801
857
    type(scalar_field), intent(in) :: t
802
858
    type(scalar_field), intent(in) :: density, olddensity
803
859
    type(scalar_field), intent(in) :: porosity_theta    
 
860
    type(scalar_field), intent(in) :: nvfrac    
804
861
    real, dimension(ele_ngi(t, ele)), intent(in) :: detwei, detwei_old, detwei_new
805
862
    real, dimension(ele_loc(t, ele), ele_loc(t, ele)), intent(inout) :: matrix_addto
806
863
    real, dimension(ele_loc(t, ele)), intent(inout) :: rhs_addto
827
884
      else
828
885
        if (include_porosity) then
829
886
          mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad*porosity_theta_at_quad)        
 
887
        else if(multiphase) then
 
888
          mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad*ele_val_at_quad(nvfrac, ele))
830
889
        else
831
890
          mass_matrix = shape_shape(test_function, ele_shape(t, ele), detwei*density_at_quad)
832
891
        end if
877
936
  
878
937
  subroutine add_advection_element_cg(ele, test_function, t, &
879
938
                                velocity, grid_velocity, diffusivity, &
880
 
                                density, olddensity, &
881
 
                                dt_t, du_t, dug_t, drho_t, detwei, j_mat, matrix_addto, rhs_addto)
 
939
                                density, olddensity, nvfrac, &
 
940
                                dt_t, du_t, dug_t, drho_t, dnvfrac_t, detwei, j_mat, matrix_addto, rhs_addto)
882
941
    integer, intent(in) :: ele
883
942
    type(element_type), intent(in) :: test_function
884
943
    type(scalar_field), intent(in) :: t
886
945
    type(vector_field), pointer :: grid_velocity
887
946
    type(tensor_field), intent(in) :: diffusivity
888
947
    type(scalar_field), intent(in) :: density, olddensity
 
948
    type(scalar_field), intent(in) :: nvfrac
889
949
    real, dimension(ele_loc(t, ele), ele_ngi(t, ele), mesh_dim(t)), intent(in) :: dt_t
890
950
    real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: du_t
891
951
    real, dimension(:, :, :) :: dug_t
892
952
    real, dimension(ele_loc(density, ele), ele_ngi(density, ele), mesh_dim(density)), intent(in) :: drho_t
 
953
    real, dimension(:, :, :), intent(in) :: dnvfrac_t
893
954
    real, dimension(ele_ngi(t, ele)), intent(in) :: detwei
894
955
    real, dimension(mesh_dim(t), mesh_dim(t), ele_ngi(t, ele)), intent(in) :: j_mat 
895
956
    real, dimension(ele_loc(t, ele), ele_loc(t, ele)), intent(inout) :: matrix_addto
904
965
    real, dimension(velocity%dim, ele_ngi(density, ele)) :: densitygrad_at_quad
905
966
    real, dimension(ele_ngi(density, ele)) :: udotgradrho_at_quad
906
967
    
 
968
    real, dimension(ele_ngi(t, ele)) :: nvfrac_at_quad
 
969
    real, dimension(velocity%dim, ele_ngi(t, ele)) :: nvfracgrad_at_quad
 
970
    real, dimension(ele_ngi(t, ele)) :: udotgradnvfrac_at_quad
 
971
    
907
972
    assert(have_advection)
908
973
    
909
974
    t_shape => ele_shape(t, ele)
922
987
      densitygrad_at_quad = density_theta*ele_grad_at_quad(density, ele, drho_t) &
923
988
                           +(1.-density_theta)*ele_grad_at_quad(olddensity, ele, drho_t)
924
989
      udotgradrho_at_quad = sum(densitygrad_at_quad*velocity_at_quad, 1)
 
990
      
 
991
      if(multiphase) then
 
992
         nvfrac_at_quad = ele_val_at_quad(nvfrac, ele)
 
993
      end if
925
994
    end select
926
995
                
927
996
    if(integrate_advection_by_parts) then
931
1000
      !    /                                        /
932
1001
      select case(equation_type)
933
1002
      case(FIELD_EQUATION_INTERNALENERGY)
934
 
        advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei*density_at_quad)
 
1003
        if(multiphase) then      
 
1004
           advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei*density_at_quad*nvfrac_at_quad)
 
1005
        else
 
1006
           advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei*density_at_quad)
 
1007
        end if
 
1008
           
935
1009
        if(abs(1.0 - beta) > epsilon(0.0)) then
936
1010
          velocity_div_at_quad = ele_div_at_quad(velocity, ele, du_t)
937
 
          advection_mat = advection_mat &
938
 
                    - (1.0-beta) * shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad &
939
 
                                                                      +udotgradrho_at_quad)* detwei)
 
1011
          
 
1012
          if(multiphase) then
 
1013
          
 
1014
             nvfracgrad_at_quad = ele_grad_at_quad(nvfrac, ele, dnvfrac_t)
 
1015
             udotgradnvfrac_at_quad = sum(nvfracgrad_at_quad*velocity_at_quad, 1)
 
1016
             
 
1017
             advection_mat = advection_mat - (1.0-beta) * ( shape_shape(test_function, t_shape, detwei*velocity_div_at_quad*density_at_quad*nvfrac_at_quad) &
 
1018
                                                           + shape_shape(test_function, t_shape, detwei*nvfrac_at_quad*udotgradrho_at_quad) &
 
1019
                                                           + shape_shape(test_function, t_shape, detwei*density_at_quad*udotgradnvfrac_at_quad) )
 
1020
          else
 
1021
             advection_mat = advection_mat - (1.0-beta) * shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad + udotgradrho_at_quad)*detwei)
 
1022
          end if
940
1023
        end if
941
1024
      case default
942
1025
        advection_mat = -dshape_dot_vector_shape(dt_t, velocity_at_quad, t_shape, detwei)
953
1036
      !  /                                 /
954
1037
      select case(equation_type)
955
1038
      case(FIELD_EQUATION_INTERNALENERGY)
956
 
        advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei*density_at_quad)
 
1039
        
 
1040
        if(multiphase) then     
 
1041
           ! vfrac*rho*nu*grad(internalenergy)
 
1042
           advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei*density_at_quad*nvfrac_at_quad)
 
1043
        else
 
1044
           advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei*density_at_quad)
 
1045
        end if
 
1046
        
957
1047
        if(abs(beta) > epsilon(0.0)) then
958
1048
          velocity_div_at_quad = ele_div_at_quad(velocity, ele, du_t)
959
 
          advection_mat = advection_mat &
960
 
                    + beta*shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad &
961
 
                                                                +udotgradrho_at_quad)*detwei)
 
1049
          
 
1050
          if(multiphase) then
 
1051
             ! advection_mat + internalenergy*div(vfrac*rho*nu)
 
1052
             ! Split up div(vfrac*rho*nu) = vfrac*rho*div(nu) + nu*grad(vfrac*rho) = vfrac*rho*div(nu) + nu*(vfrac*grad(rho) + rho*grad(nvfrac))
 
1053
             
 
1054
             nvfracgrad_at_quad = ele_grad_at_quad(nvfrac, ele, dnvfrac_t)
 
1055
             udotgradnvfrac_at_quad = sum(nvfracgrad_at_quad*velocity_at_quad, 1)
 
1056
             
 
1057
             advection_mat = advection_mat + beta * ( shape_shape(test_function, t_shape, detwei*velocity_div_at_quad*density_at_quad*nvfrac_at_quad) &
 
1058
                                                     + shape_shape(test_function, t_shape, detwei*nvfrac_at_quad*udotgradrho_at_quad) &
 
1059
                                                     + shape_shape(test_function, t_shape, detwei*density_at_quad*udotgradnvfrac_at_quad) )
 
1060
          else
 
1061
             advection_mat = advection_mat + beta*shape_shape(test_function, t_shape, (velocity_div_at_quad*density_at_quad &
 
1062
                                                              +udotgradrho_at_quad)*detwei)
 
1063
          end if
962
1064
        end if
963
1065
      case default
964
1066
        advection_mat = shape_vector_dot_dshape(test_function, velocity_at_quad, dt_t, detwei)
1060
1162
    
1061
1163
  end subroutine add_diffusivity_element_cg
1062
1164
  
1063
 
  subroutine add_pressurediv_element_cg(ele, test_function, t, velocity, pressure, du_t, detwei, rhs_addto)
 
1165
  subroutine add_pressurediv_element_cg(ele, test_function, t, velocity, pressure, nvfrac, du_t, dnvfrac_t, detwei, rhs_addto)
1064
1166
  
1065
1167
    integer, intent(in) :: ele
1066
1168
    type(element_type), intent(in) :: test_function
1067
1169
    type(scalar_field), intent(in) :: t
1068
1170
    type(vector_field), intent(in) :: velocity
1069
1171
    type(scalar_field), intent(in) :: pressure
1070
 
    real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)) :: du_t
 
1172
    type(scalar_field), intent(in) :: nvfrac
 
1173
    real, dimension(ele_loc(velocity, ele), ele_ngi(velocity, ele), mesh_dim(t)), intent(in) :: du_t
 
1174
    real, dimension(:, :, :), intent(in) :: dnvfrac_t
1071
1175
    real, dimension(ele_ngi(t, ele)), intent(in) :: detwei
1072
1176
    real, dimension(ele_loc(t, ele)), intent(inout) :: rhs_addto
1073
1177
    
 
1178
    real, dimension(velocity%dim, ele_ngi(t, ele)) :: nvfracgrad_at_quad
 
1179
    real, dimension(ele_ngi(t, ele)) :: udotgradnvfrac_at_quad
 
1180
    
1074
1181
    assert(equation_type==FIELD_EQUATION_INTERNALENERGY)
1075
1182
    assert(ele_ngi(pressure, ele)==ele_ngi(t, ele))
1076
1183
    
1077
 
    rhs_addto = rhs_addto - &
1078
 
                shape_rhs(test_function, ele_div_at_quad(velocity, ele, du_t) * ele_val_at_quad(pressure, ele) * detwei)
 
1184
    if(multiphase) then
 
1185
       ! -p * (div(vfrac*nu)) 
 
1186
       nvfracgrad_at_quad = ele_grad_at_quad(nvfrac, ele, dnvfrac_t)
 
1187
       udotgradnvfrac_at_quad = sum(nvfracgrad_at_quad*ele_val_at_quad(velocity, ele), 1)
 
1188
             
 
1189
       rhs_addto = rhs_addto - ( shape_rhs(test_function, ele_div_at_quad(velocity, ele, du_t) * ele_val_at_quad(nvfrac, ele) * ele_val_at_quad(pressure, ele) * detwei) &
 
1190
                                 + shape_rhs(test_function, udotgradnvfrac_at_quad * ele_val_at_quad(pressure, ele) * detwei) )
 
1191
    else
 
1192
       rhs_addto = rhs_addto - shape_rhs(test_function, ele_div_at_quad(velocity, ele, du_t) * ele_val_at_quad(pressure, ele) * detwei)
 
1193
    end if
1079
1194
    
1080
1195
  end subroutine add_pressurediv_element_cg
1081
1196
  
1082
 
  subroutine assemble_advection_diffusion_face_cg(face, bc_type, t, t_bc, t_bc_2, matrix, rhs, positions, velocity, grid_velocity, density, olddensity)
 
1197
  subroutine assemble_advection_diffusion_face_cg(face, bc_type, t, t_bc, t_bc_2, matrix, rhs, positions, velocity, grid_velocity, density, olddensity, nvfrac)
1083
1198
    integer, intent(in) :: face
1084
1199
    integer, intent(in) :: bc_type
1085
1200
    type(scalar_field), intent(in) :: t
1092
1207
    type(vector_field), pointer :: grid_velocity
1093
1208
    type(scalar_field), intent(in) :: density
1094
1209
    type(scalar_field), intent(in) :: olddensity
 
1210
    type(scalar_field), intent(in) :: nvfrac
1095
1211
    
1096
1212
    integer, dimension(face_loc(t, face)) :: face_nodes
1097
1213
    real, dimension(face_ngi(t, face)) :: detwei
1125
1241
    
1126
1242
    ! Advection
1127
1243
    if(have_advection .and. integrate_advection_by_parts) &
1128
 
      call add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, detwei, normal, matrix_addto, rhs_addto)
 
1244
      call add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, nvfrac, detwei, normal, matrix_addto, rhs_addto)
1129
1245
    
1130
1246
    ! Diffusivity
1131
1247
    if(have_diffusivity) call add_diffusivity_face_cg(face, bc_type, t, t_bc, t_bc_2, detwei, matrix_addto, rhs_addto)
1138
1254
    
1139
1255
  end subroutine assemble_advection_diffusion_face_cg
1140
1256
  
1141
 
  subroutine add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, detwei, normal, matrix_addto, rhs_addto)
 
1257
  subroutine add_advection_face_cg(face, bc_type, t, t_bc, velocity, grid_velocity, density, olddensity, nvfrac, detwei, normal, matrix_addto, rhs_addto)
1142
1258
    integer, intent(in) :: face
1143
1259
    integer, intent(in) :: bc_type
1144
1260
    type(scalar_field), intent(in) :: t
1147
1263
    type(vector_field), pointer :: grid_velocity
1148
1264
    type(scalar_field), intent(in) :: density
1149
1265
    type(scalar_field), intent(in) :: olddensity
 
1266
    type(scalar_field), intent(in) :: nvfrac
1150
1267
    real, dimension(face_ngi(t, face)), intent(in) :: detwei
1151
1268
    real, dimension(mesh_dim(t), face_ngi(t, face)), intent(in) :: normal
1152
1269
    real, dimension(face_loc(t, face), face_loc(t, face)), intent(inout) :: matrix_addto
1171
1288
    case(FIELD_EQUATION_INTERNALENERGY)
1172
1289
      density_at_quad = density_theta*face_val_at_quad(density, face) &
1173
1290
                       +(1.0-density_theta)*face_val_at_quad(olddensity, face)
1174
 
 
1175
 
      advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1) * density_at_quad)
 
1291
                       
 
1292
      if(multiphase) then
 
1293
         advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1) * density_at_quad * face_val_at_quad(nvfrac, face))
 
1294
      else
 
1295
         advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1) * density_at_quad)
 
1296
      end if
1176
1297
    case default
1177
1298
      
1178
1299
      advection_mat = shape_shape(t_shape, t_shape, detwei * sum(velocity_at_quad * normal, 1))