356
358
end subroutine assemble_divergence_matrix_cg
358
subroutine assemble_compressible_divergence_matrix_cg(ctp_m, state, ct_rhs)
360
subroutine assemble_compressible_divergence_matrix_cg(ctp_m, state, istate, ct_rhs, div_mass)
361
363
! bucket full of fields
362
364
type(state_type), dimension(:), intent(inout) :: state
365
integer, intent(in) :: istate
364
366
! the compressible divergence matrix
365
367
type(block_csr_matrix), pointer :: ctp_m
367
369
type(scalar_field), intent(inout), optional :: ct_rhs
369
if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then
371
call assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state(1), ct_rhs)
371
type(csr_matrix), intent(inout), optional :: div_mass
373
if(option_count("/material_phase/vector_field::Velocity/prognostic") > 1) then
375
call assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, istate, ct_rhs, div_mass)
375
FLExit("Multimaterial compressible continuous_galerkin pressure not possible.")
379
if((size(state)==1).and.(.not.has_scalar_field(state(1), "MaterialVolumeFraction"))) then
381
call assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, 1, ct_rhs, div_mass)
385
FLExit("Multimaterial compressible continuous_galerkin pressure not possible.")
379
393
end subroutine assemble_compressible_divergence_matrix_cg
381
subroutine assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, ct_rhs)
395
subroutine assemble_1mat_compressible_divergence_matrix_cg(ctp_m, state, istate, ct_rhs, div_mass)
384
398
! bucket full of fields
385
type(state_type), intent(inout) :: state
399
type(state_type), dimension(:), intent(inout) :: state
400
integer, intent(in) :: istate
387
402
! the compressible divergence matrix
388
403
type(block_csr_matrix), pointer :: ctp_m
390
405
type(scalar_field), intent(inout), optional :: ct_rhs
407
type(csr_matrix), intent(inout), optional :: div_mass
393
410
type(mesh_type), pointer :: test_mesh
408
425
real, dimension(:), allocatable :: density_at_quad, olddensity_at_quad
409
426
real, dimension(:,:), allocatable :: density_grad_at_quad, nlvelocity_at_quad
428
real, dimension(:,:), allocatable :: div_mass_mat
412
431
integer :: ele, sele, dim
425
444
integer, dimension(:), allocatable :: density_bc_type
426
445
type(scalar_field) :: density_bc
449
!! Multiphase variables
450
! Volume fraction fields
451
type(scalar_field), pointer :: vfrac
452
type(scalar_field) :: nvfrac
453
type(element_type), pointer :: nvfrac_shape
454
! Transformed gradient function for the non-linear PhaseVolumeFraction.
455
real, dimension(:, :, :), allocatable :: dnvfrac_t
456
logical :: is_compressible_phase ! Is this the (single) compressible phase in a multiphase simulation?
430
458
! =============================================================
431
459
! Subroutine to construct the matrix ctp_m (a.k.a. C1/2/3TP).
434
462
ewrite(2,*) 'In assemble_1mat_compressible_divergence_matrix_cg'
436
coordinate=> extract_vector_field(state, "Coordinate")
438
density => extract_scalar_field(state, "Density")
439
olddensity => extract_scalar_field(state, "OldDensity")
441
pressure => extract_scalar_field(state, "Pressure")
443
velocity=>extract_vector_field(state, "Velocity")
444
nonlinearvelocity=>extract_vector_field(state, "NonlinearVelocity") ! maybe this should be updated after the velocity solve?
464
coordinate=> extract_vector_field(state(istate), "Coordinate")
466
density => extract_scalar_field(state(istate), "Density")
467
olddensity => extract_scalar_field(state(istate), "OldDensity")
469
if(have_option(trim(state(istate)%option_path)//"/equation_of_state/compressible")) then
470
is_compressible_phase = .true.
472
is_compressible_phase = .false.
474
! Find the compressible phase and extract it's density to form rho_c * div(vfrac_i*u_i), where _c and _i represent
475
! the compressible and incompressible phases respectively.
476
do i = 1, size(state)
477
if(have_option(trim(state(i)%option_path)//"/equation_of_state/compressible")) then
478
density => extract_scalar_field(state(i), "Density")
479
olddensity => extract_scalar_field(state(i), "OldDensity")
484
pressure => extract_scalar_field(state(istate), "Pressure")
486
velocity=>extract_vector_field(state(istate), "Velocity")
487
nonlinearvelocity=>extract_vector_field(state(istate), "NonlinearVelocity") ! maybe this should be updated after the velocity solve?
489
! Get the non-linear PhaseVolumeFraction field if multiphase
491
vfrac => extract_scalar_field(state(istate), "PhaseVolumeFraction")
492
call allocate(nvfrac, vfrac%mesh, "NonlinearPhaseVolumeFraction")
494
call get_nonlinear_volume_fraction(state(istate), nvfrac)
495
ewrite_minmax(nvfrac)
446
500
integrate_by_parts=have_option(trim(complete_field_path(density%option_path, stat=stat))//&
447
501
&"/spatial_discretisation/continuous_galerkin/advection_terms/integrate_advection_by_parts")&
448
502
.or. have_option(trim(complete_field_path(velocity%option_path, stat=stat))//&
494
548
olddensity_at_quad(ele_ngi(density, 1)), &
495
549
nlvelocity_at_quad(nonlinearvelocity%dim, ele_ngi(nonlinearvelocity, 1)), &
496
550
density_grad_at_quad(field%dim, ele_ngi(density,1)), &
497
j_mat(field%dim, field%dim, ele_ngi(density, 1)))
551
j_mat(field%dim, field%dim, ele_ngi(density, 1)), &
552
div_mass_mat(ele_loc(test_mesh, 1), ele_loc(test_mesh, 1)))
555
! We will need grad(nvfrac) if we are not integrating by parts below
556
allocate(dnvfrac_t(ele_loc(nvfrac,1), ele_ngi(nvfrac,1), field%dim))
499
559
do ele=1, element_count(test_mesh)
501
561
test_nodes=>ele_nodes(test_mesh, ele)
509
569
olddensity_at_quad = ele_val_at_quad(olddensity, ele)
511
571
nlvelocity_at_quad = ele_val_at_quad(nonlinearvelocity, ele)
513
573
if(any(stabilisation_scheme == (/STABILISATION_STREAMLINE_UPWIND, STABILISATION_SUPG/))) then
514
574
call transform_to_physical(coordinate, ele, test_shape_ptr, dshape = dtest_t, &
515
575
detwei = detwei, j = j_mat)
517
577
call transform_to_physical(coordinate, ele, test_shape_ptr, dshape = dtest_t, detwei=detwei)
520
if(.not.integrate_by_parts) then
580
if(.not.integrate_by_parts .or. (multiphase .and. integrate_by_parts .and. .not.is_compressible_phase)) then
521
581
! transform the field (velocity) derivatives into physical space
522
582
call transform_to_physical(coordinate, ele, field_shape, dshape=dfield_t)
547
607
if(integrate_by_parts) then
548
609
! if SUPG is fixed for P>1 then this dtest_t should be updated
549
ele_mat = -dshape_shape(dtest_t, field_shape, &
550
detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad))
610
if(multiphase .and. .not.is_compressible_phase) then
611
density_grad_at_quad = theta*(ele_grad_at_quad(density, ele, ddensity_t))+&
612
(1-theta)*(ele_grad_at_quad(olddensity, ele, ddensity_t))
614
ele_mat = -dshape_shape(dtest_t, field_shape, detwei*ele_val_at_quad(nvfrac, ele)*(theta*density_at_quad + (1-theta)*olddensity_at_quad)) - shape_shape_vector(test_shape, field_shape, detwei*ele_val_at_quad(nvfrac, ele), density_grad_at_quad)
616
else if(multiphase .and. is_compressible_phase) then
617
ele_mat = -dshape_shape(dtest_t, field_shape, detwei*ele_val_at_quad(nvfrac, ele)*(theta*density_at_quad + (1-theta)*olddensity_at_quad))
621
ele_mat = -dshape_shape(dtest_t, field_shape, detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad))
552
624
density_grad_at_quad = theta*(ele_grad_at_quad(density, ele, ddensity_t))+&
553
625
(1-theta)*(ele_grad_at_quad(olddensity, ele, ddensity_t))
555
ele_mat = shape_dshape(test_shape, dfield_t, &
629
! If the field and nvfrac meshes are different, then we need to
630
! compute the derivatives of the nvfrac shape functions.
631
if(.not.(nvfrac%mesh == field%mesh)) then
632
nvfrac_shape => ele_shape(nvfrac%mesh, ele)
633
call transform_to_physical(coordinate, ele, nvfrac_shape, dshape=dnvfrac_t)
638
! Split up the divergence term div(rho*vfrac*u) = rho*div(u*vfrac) + vfrac*u*grad(rho)
639
! = (rho*vfrac*div(u) + rho*u*grad(vfrac)) + u*vfrac*grad(rho)
641
! First assemble rho*div(u*vfrac). This is the incompressible phase's divergence matrix.
642
ele_mat = shape_dshape(test_shape, dfield_t, detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad)*ele_val_at_quad(nvfrac, ele)) + shape_shape_vector(test_shape, field_shape, detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad), ele_grad_at_quad(nvfrac, ele, dnvfrac_t))
644
! If the phase is compressible, then we now complete the assembly of div(rho*vfrac*u) below.
645
if(is_compressible_phase) then
646
ele_mat = ele_mat + shape_shape_vector(test_shape, field_shape, detwei*ele_val_at_quad(nvfrac, ele), density_grad_at_quad)
650
ele_mat = shape_dshape(test_shape, dfield_t, &
556
651
detwei*(theta*density_at_quad + (1-theta)*olddensity_at_quad)) + &
557
652
shape_shape_vector(test_shape, field_shape, detwei, density_grad_at_quad)
560
657
! Stabilisation does not return the right shape for this operator!
568
665
do dim = 1, field%dim
569
666
call addto(ctp_m, 1, dim, test_nodes, field_nodes, ele_mat(dim,:,:))
669
if(present(div_mass)) then
670
div_mass_mat = shape_shape(test_shape, test_shape, detwei)
671
call addto(div_mass, test_nodes, test_nodes, div_mass_mat)
572
674
call deallocate(test_shape)
618
720
& detwei_f=detwei_bdy,&
619
721
& normal=normal_bdy)
621
ele_mat_bdy = shape_shape_vector(test_shape_ptr, field_shape, &
622
detwei_bdy*density_bdy, normal_bdy)
725
ele_mat_bdy = shape_shape_vector(test_shape_ptr, field_shape, &
726
detwei_bdy*density_bdy*face_val_at_quad(nvfrac, sele), normal_bdy)
728
ele_mat_bdy = shape_shape_vector(test_shape_ptr, field_shape, &
729
detwei_bdy*density_bdy, normal_bdy)
624
732
do dim = 1, field%dim
625
733
if((field_bc_type(dim, sele)==1).and.present(ct_rhs)) then
added
removed
assemble/Divergence_Matrix_CG.F90
