riemann2d.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: riemann2d.f90 #
5!# #
6!# Copyright (C) 2006-2024 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# #
9!# This program is free software; you can redistribute it and/or modify #
10!# it under the terms of the GNU General Public License as published by #
11!# the Free Software Foundation; either version 2 of the License, or (at #
12!# your option) any later version. #
13!# #
14!# This program is distributed in the hope that it will be useful, but #
15!# WITHOUT ANY WARRANTY; without even the implied warranty of #
16!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
17!# NON INFRINGEMENT. See the GNU General Public License for more #
18!# details. #
19!# #
20!# You should have received a copy of the GNU General Public License #
21!# along with this program; if not, write to the Free Software #
22!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
23!# #
24!#############################################################################
466 USE fosite_mod
467#include "tap.h"
468 IMPLICIT NONE
469 !--------------------------------------------------------------------------!
470 ! simulation parameter
471 INTEGER, PARAMETER :: select_test = 6 ! run this test, should be in 1..19
472 LOGICAL, PARAMETER :: run_all = .true. ! set this if you want to run all 19 tests
473 REAL, PARAMETER :: gamma = 1.4 ! ratio of specific heats
474 ! mesh settings
475 INTEGER, PARAMETER :: mgeo = cartesian ! geometry of the mesh
476! INTEGER, PARAMETER :: MGEO = CYLINDRICAL
477! INTEGER, PARAMETER :: MGEO = LOGCYLINDRICAL
478 INTEGER, PARAMETER :: xres = 100 ! resolution
479 INTEGER, PARAMETER :: yres = 100
480 INTEGER, PARAMETER :: zres = 1
481 REAL, PARAMETER :: rmin = 1.0e-4 ! inner radius for polar grids
482 ! output file parameter
483 INTEGER, PARAMETER :: onum = 1 ! number of output data sets
484 CHARACTER(LEN=256), PARAMETER & ! output data dir
485 :: odir = './'
486 CHARACTER(LEN=256), PARAMETER & ! output data file name
487 :: ofname = 'riemann2d'
488 ! do not change these parameters
489 INTEGER, PARAMETER :: num_tests = 19 ! total number of 2D setups
490 REAL, PARAMETER, DIMENSION(NUM_TESTS) :: & ! stop time for each test
491 test_stoptime = (/ 0.2,0.2,0.3,0.25,0.23,0.3,0.25, &
492 0.25,0.3,0.15,0.3,0.25,0.3,0.1, &
493 0.2,0.2,0.3,0.2,0.3/)
494 !--------------------------------------------------------------------------!
495 CLASS(fosite), ALLOCATABLE :: sim
496 LOGICAL :: ok(num_tests), mpifinalize = .false.
497 INTEGER :: ic,err
498 !--------------------------------------------------------------------------!
499
500 IF (run_all) THEN
501tap_plan(num_tests)
502 ELSE
503tap_plan(1)
504 END IF
505
506 DO ic=1,num_tests
507 IF (.NOT.run_all.AND.ic.NE.select_test) cycle ! skip if not selected
508
509 IF (ALLOCATED(sim)) DEALLOCATE(sim)
510 ALLOCATE(sim,stat=err)
511 IF (err.NE.0) &
512 CALL sim%Error("riemann2d","cannot allocate memory for Sim")
513
514 CALL sim%InitFosite()
515 CALL makeconfig(sim, sim%config, ic)
516
517! CALL PrintDict(config)
518
519 CALL sim%Setup()
520
521 ! set initial condition
522 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc, ic)
523 CALL sim%Run()
524 ok(ic) = .NOT.sim%aborted
525
526 IF ((run_all.AND.ic.EQ.num_tests).OR.(.NOT.run_all)) mpifinalize=.true.
527 CALL sim%Finalize(mpifinalize)
528
529 END DO
530
531 DO ic=1,num_tests
532 IF (.NOT.run_all.AND.ic.NE.select_test) cycle ! skip if not selected
533tap_check(ok(ic),"Simulation finished")
534 END DO
535
536 IF (ALLOCATED(sim)) DEALLOCATE(sim)
537
538tap_done
539
540CONTAINS
541
542 SUBROUTINE makeconfig(Sim, config, ic)
543 USE functions, ONLY : asinh
544 IMPLICIT NONE
545 !------------------------------------------------------------------------!
546 CLASS(fosite) :: Sim
547 TYPE(dict_typ),POINTER :: config
548 INTEGER :: ic
549 !------------------------------------------------------------------------!
550 ! Local variable declaration
551 INTEGER :: bc(6)
552 TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, &
553 timedisc, fluxes
554 REAL :: x1,x2,y1,y2,z1,z2,sc
555 CHARACTER(LEN=16) :: fext
556 !------------------------------------------------------------------------!
557 INTENT(INOUT) :: sim
558 INTENT(IN) :: ic
559 CHARACTER(LEN=8) :: geo_str
560 !------------------------------------------------------------------------!
561 IF (ic.LT.1.OR.ic.GT.num_tests) &
562 CALL sim%Error("riemann2d::MakeConfig","invalid test number selected")
563
564 ! mesh settings and boundary conditions
565 SELECT CASE(mgeo)
566 CASE(cartesian)
567 sc = 1.0
568 x1 = -0.5
569 x2 = 0.5
570 y1 = -0.5
571 y2 = 0.5
572 z1 = -0.0
573 z2 = 0.0
574 bc(west) = no_gradients
575 bc(east) = no_gradients
576 bc(south) = no_gradients
577 bc(north) = no_gradients
578 bc(bottom)= no_gradients
579 bc(top) = no_gradients
580 geo_str = "-cart"
581 CASE(cylindrical)
582 sc = 1.0
583 x1 = 0.0
584 x2 = 0.5*sqrt(2.0)
585 y1 = 0.0
586 y2 = 2*pi
587 z1 = -0.0
588 z2 = 0.0
589 bc(west) = absorbing
590 bc(east) = no_gradients
591 bc(south) = periodic
592 bc(north) = periodic
593 bc(bottom) = reflecting
594 bc(top) = reflecting
595 geo_str = "-cyl"
596 CASE(logcylindrical)
597 sc = 0.3
598 x1 = log(rmin/sc)
599 x2 = log(0.5*sqrt(2.0)/sc)
600 y1 = -pi
601 y2 = pi
602 z1 = 0.0
603 z2 = 0.0
604 bc(west) = axis
605 bc(east) = no_gradients
606 bc(south) = periodic
607 bc(north) = periodic
608 bc(bottom) = no_gradients
609 bc(top) = no_gradients
610 geo_str = "-logcyl"
611 CASE DEFAULT
612 CALL sim%Error("riemann2d::MakeConfig","geometry not supported in this test")
613 END SELECT
614
615 ! mesh settings
616 mesh => dict("meshtype" / midpoint, &
617 "geometry" / mgeo, &
618 "inum" / xres, &
619 "jnum" / yres, &
620 "knum" / zres, &
621 "xmin" / x1, &
622 "xmax" / x2, &
623 "ymin" / y1, &
624 "ymax" / y2, &
625 "zmin" / z1, &
626 "zmax" / z2, &
627 "output/commutator" / 1,&
628 "gparam" / sc)
629
630 ! boundary conditions
631 boundary => dict("western" / bc(west), &
632 "eastern" / bc(east), &
633 "southern" / bc(south), &
634 "northern" / bc(north), &
635 "bottomer" / bc(bottom),&
636 "topper" / bc(top))
637
638 ! physics settings
639 physics => dict("problem" / euler, &
640 "gamma" / gamma) ! ratio of specific heats !
641
642 ! flux calculation and reconstruction method
643 fluxes => dict("order" / linear, &
644 "fluxtype" / kt, &
645 "variables" / conservative, & ! vars. to use for reconstruction!
646 "limiter" / monocent, & ! one of: minmod, monocent,... !
647 "theta" / 1.2) ! optional parameter for limiter !
648
649 ! time discretization settings
650 timedisc => dict( &
651 "method" / modified_euler, &
652 "order" / 3, &
653 "cfl" / 0.4, &
654 "stoptime" / test_stoptime(ic), &
655 "dtlimit" / 1.0e-10, &
656 "maxiter" / 10000000, &
657 "output/geometrical_sources" / 1)
658
659 ! initialize data input/output
660 ! append test number to file names
661 WRITE (fext, '(A,I2.2)') trim(geo_str) // "_", ic
662 datafile => dict( &
663! "fileformat" / GNUPLOT, "filecycles" / 0, &
664 "fileformat" / vtk, &
665! "fileformat" / XDMF, &
666 "filename" / (trim(odir) // trim(ofname) // trim(fext)), &
667 "count" / onum)
668
669 config => dict("mesh" / mesh, &
670 "physics" / physics, &
671 "boundary" / boundary, &
672 "fluxes" / fluxes, &
673 "timedisc" / timedisc, &
674! "logfile" / logfile, &
675 "datafile" / datafile)
676 END SUBROUTINE makeconfig
677
679 SUBROUTINE initdata(Mesh,Physics,Timedisc,ic)
680 IMPLICIT NONE
681 !------------------------------------------------------------------------!
682 CLASS(physics_base) :: Physics
683 CLASS(mesh_base) :: Mesh
684 CLASS(timedisc_base):: Timedisc
685 INTEGER :: ic
686 !------------------------------------------------------------------------!
687 ! Local variable declaration
688 TYPE(marray_base) :: vcart,vcurv
689 REAL :: xmin,ymin,xmax,ymax,x0,y0
690#ifdef PARALLEL
691 REAL :: xmin_all,xmax_all,ymin_all,ymax_all
692 INTEGER :: ierr
693#endif
694 REAL, DIMENSION(NUM_TESTS,4):: den,pre,vx,vy
695 CHARACTER(LEN=64) :: teststr
696 !------------------------------------------------------------------------!
697 INTENT(IN) :: mesh,physics,ic
698 INTENT(INOUT) :: timedisc
699 !------------------------------------------------------------------------!
700 ! minima and maxima of _cartesian_ coordinates
701 xmin = minval(mesh%bccart(:,:,:,1))
702 ymin = minval(mesh%bccart(:,:,:,2))
703 xmax = maxval(mesh%bccart(:,:,:,1))
704 ymax = maxval(mesh%bccart(:,:,:,2))
705#ifdef PARALLEL
706 CALL mpi_allreduce(xmin,xmin_all,1,default_mpi_real,mpi_min,mesh%comm_cart,ierr)
707 xmin = xmin_all
708 CALL mpi_allreduce(ymin,ymin_all,1,default_mpi_real,mpi_min,mesh%comm_cart,ierr)
709 ymin = ymin_all
710 CALL mpi_allreduce(xmax,xmax_all,1,default_mpi_real,mpi_max,mesh%comm_cart,ierr)
711 xmax = xmax_all
712 CALL mpi_allreduce(ymax,ymax_all,1,default_mpi_real,mpi_max,mesh%comm_cart,ierr)
713 ymax = ymax_all
714#endif
715 x0 = xmin + 0.5*abs(xmax-xmin)
716 y0 = ymin + 0.5*abs(ymax-ymin)
717
718 vcart = marray_base(3)
719 vcurv = marray_base(3)
720
721 ! set defaults
722 den(:,:) = 1.
723 vx(:,:) = 0.
724 vy(:,:) = 0.
725 pre(:,:) = 1.
726
727 IF (ic.LT.1.OR.ic.GT.num_tests) &
728 CALL mesh%Error("InitData","Sorry, this 2D Riemann problem is currently not supported!")
729
730 WRITE (teststr,'(A,I2)') "2D Riemann problem no. ", ic
731
732 ! initial conditions
733
734 ! Test configuration no. 1
735 ! 1st quadrant
736 den(1,1) = 1.
737 pre(1,1) = 1.
738 ! 2nd quadrant
739 den(1,2) = .5197
740 vx(1,2) = -.7259
741 pre(1,2) = .4
742 ! 3rd quadrant
743 den(1,3) = .1072
744 vx(1,3) = -.7259
745 vy(1,3) = -1.4045
746 pre(1,3) = .0439
747 ! 4th quadrant
748 den(1,4) = .2579
749 vy(1,4) = -1.4045
750 pre(1,4) = .15
751
752 ! Test configuration no. 2
753 ! 1st quadrant
754 den(2,1) = 1.
755 pre(2,1) = 1.
756 ! 2nd quadrant
757 den(2,2) = .5197
758 vx(2,2) = -.7259
759 pre(2,2) = .4
760 ! 3rd quadrant
761 den(2,3) = 1.
762 vx(2,3) = -.7259
763 vy(2,3) = -.7259
764 pre(2,3) = 1.
765 ! 4th quadrant
766 den(2,4) = .5197
767 vy(2,4) = -.7259
768 pre(2,4) = .4
769
770 ! Test configuration no. 3
771 ! 1st quadrant
772 den(3,1) = 1.5
773 pre(3,1) = 1.5
774 ! 2nd quadrant
775 den(3,2) = .5323
776 vx(3,2) = 1.206
777 pre(3,2) = .3
778 ! 3rd quadrant
779 den(3,3) = 0.138
780 vx(3,3) = 1.206
781 vy(3,3) = 1.206
782 pre(3,3) = 0.029
783 ! 4th quadrant
784 den(3,4) = .5323
785 vy(3,4) = 1.206
786 pre(3,4) = .3
787
788 ! Test configuration no. 4
789 ! 1st quadrant
790 den(4,1) = 1.1
791 pre(4,1) = 1.1
792 ! 2nd quadrant
793 den(4,2) = .5065
794 vx(4,2) = .8939
795 pre(4,2) = .35
796 ! 3rd quadrant
797 den(4,3) = 1.1
798 vx(4,3) = .8939
799 vy(4,3) = .8939
800 pre(4,3) = 1.1
801 ! 4th quadrant
802 den(4,4) = .5065
803 vy(4,4) = .8939
804 pre(4,4) = .35
805
806 ! 1st quadrant
807 ! Test configuration no. 5
808 den(5,1) = 1.
809 vx(5,1) = -.75
810 vy(5,1) = -.5
811 pre(5,1) = 1.
812 ! 2nd quadrant
813 den(5,2) = 2.
814 vx(5,2) = -.75
815 vy(5,2) = .5
816 pre(5,2) = 1.
817 ! 3rd quadrant
818 den(5,3) = 1.
819 vx(5,3) = .75
820 vy(5,3) = .5
821 pre(5,3) = 1.
822 ! 4th quadrant
823 den(5,4) = 3.
824 vx(5,4) = .75
825 vy(5,4) = -.5
826 pre(5,4) = 1.
827
828 ! Test configuration no. 6
829 ! 1st quadrant
830 den(6,1) = 1.
831 vx(6,1) = 0.75
832 vy(6,1) = -0.5
833 pre(6,1) = 1.
834 ! 2nd quadrant
835 den(6,2) = 2.
836 vx(6,2) = 0.75
837 vy(6,2) = 0.5
838 pre(6,2) = 1.
839 ! 3rd quadrant
840 den(6,3) = 1.
841 vx(6,3) = -0.75
842 vy(6,3) = 0.5
843 pre(6,3) = 1.
844 ! 4th quadrant
845 den(6,4) = 3.
846 vx(6,4) = -0.75
847 vy(6,4) = -0.5
848 pre(6,4) = 1.
849
850 ! Test configuration no. 7
851 ! 1st quadrant
852 den(7,1) = 1.
853 vx(7,1) = 0.1
854 vy(7,1) = 0.1
855 pre(7,1) = 1.
856 ! 2nd quadrant
857 den(7,2) = 0.5197
858 vx(7,2) = -0.6259
859 vy(7,2) = 0.1
860 pre(7,2) = 0.4
861 ! 3rd quadrant
862 den(7,3) = 0.8
863 vx(7,3) = 0.1
864 vy(7,3) = 0.1
865 pre(7,3) = 0.4
866 ! 4th quadrant
867 den(7,4) = 0.5197
868 vx(7,4) = 0.1
869 vy(7,4) = -0.6259
870 pre(7,4) = 0.4
871
872 ! Test configuration no. 8
873 ! 1st quadrant
874 den(8,1) = 0.5197
875 vx(8,1) = 0.1
876 vy(8,1) = 0.1
877 pre(8,1) = 0.4
878 ! 2nd quadrant
879 den(8,2) = 1.
880 vx(8,2) = -0.6259
881 vy(8,2) = 0.1
882 pre(8,2) = 1.
883 ! 3rd quadrant
884 den(8,3) = 0.8
885 vx(8,3) = 0.1
886 vy(8,3) = 0.1
887 pre(8,3) = 1.
888 ! 4th quadrant
889 den(8,4) = 1.
890 vx(8,4) = 0.1
891 vy(8,4) = -0.6259
892 pre(8,4) = 1.
893
894 ! Test configuration no. 9
895 ! 1st quadrant
896 den(9,1) = 1.
897 vx(9,1) = 0.
898 vy(9,1) = 0.3
899 pre(9,1) = 1.
900 ! 2nd quadrant
901 den(9,2) = 2.
902 vx(9,2) = 0.
903 vy(9,2) = -0.3
904 pre(9,2) = 1.
905 ! 3rd quadrant
906 den(9,3) = 1.039
907 vx(9,3) = 0.
908 vy(9,3) = -0.8133
909 pre(9,3) = 0.4
910 ! 4th quadrant
911 den(9,4) = 0.5197
912 vx(9,4) = 0.
913 vy(9,4) = -0.4259
914 pre(9,4) = 0.4
915
916 ! Test configuration no. 10
917 ! 1st quadrant
918 den(10,1) = 1.
919 vx(10,1) = 0.
920 vy(10,1) = 0.4297
921 pre(10,1) = 1.
922 ! 2nd quadrant
923 den(10,2) = 0.5
924 vx(10,2) = 0.
925 vy(10,2) = 0.6076
926 pre(10,2) = 1.
927 ! 3rd quadrant
928 den(10,3) = 0.2281
929 vx(10,3) = 0.
930 vy(10,3) = -0.6076
931 pre(10,3) = 0.3333
932 ! 4th quadrant
933 den(10,4) = 0.4562
934 vx(10,4) = 0.
935 vy(10,4) = -0.4297
936 pre(10,4) = 0.3333
937
938 ! Test configuration no. 11
939 ! 1st quadrant
940 den(11,1) = 1.
941 vx(11,1) = 0.1
942 vy(11,1) = 0.
943 pre(11,1) = 1.
944 ! 2nd quadrant
945 den(11,2) = 0.5313
946 vx(11,2) = 0.8276
947 vy(11,2) = 0.
948 pre(11,2) = 0.4
949 ! 3rd quadrant
950 den(11,3) = 0.8
951 vx(11,3) = 0.1
952 vy(11,3) = 0.
953 pre(11,3) = 0.4
954 ! 4th quadrant
955 den(11,4) = 0.5313
956 vx(11,4) = 0.1
957 vy(11,4) = 0.7276
958 pre(11,4) = 0.4
959
960 ! Test configuration no. 12
961 ! 1st quadrant
962 den(12,1) = 0.5313
963 vx(12,1) = 0.
964 vy(12,1) = 0.
965 pre(12,1) = 0.4
966 ! 2nd quadrant
967 den(12,2) = 1.
968 vx(12,2) = 0.7276
969 vy(12,2) = 0.
970 pre(12,2) = 1.
971 ! 3rd quadrant
972 den(12,3) = 0.8
973 vx(12,3) = 0.
974 vy(12,3) = 0.
975 pre(12,3) = 1.
976 ! 4th quadrant
977 den(12,4) = 1.
978 vx(12,4) = 0.
979 vy(12,4) = 0.7276
980 pre(12,4) = 1.
981
982 ! Test configuration no. 13
983 ! 1st quadrant
984 den(13,1) = 1.
985 vx(13,1) = 0.
986 vy(13,1) = -0.3
987 pre(13,1) = 1.
988 ! 2nd quadrant
989 den(13,2) = 2.
990 vx(13,2) = 0.
991 vy(13,2) = 0.3
992 pre(13,2) = 1.
993 ! 3rd quadrant
994 den(13,3) = 1.0625
995 vx(13,3) = 0.
996 vy(13,3) = 0.8145
997 pre(13,3) = 0.4
998 ! 4th quadrant
999 den(13,4) = 0.5313
1000 vx(13,4) = 0.
1001 vy(13,4) = 0.4276
1002 pre(13,4) = 0.4
1003
1004 ! Test configuration no. 14
1005 ! 1st quadrant
1006 den(14,1) = 2.
1007 vx(14,1) = 0.
1008 vy(14,1) = -0.5606
1009 pre(14,1) = 8.
1010 ! 2nd quadrant
1011 den(14,2) = 1.
1012 vx(14,2) = 0.
1013 vy(14,2) = -1.2172
1014 pre(14,2) = 8.
1015 ! 3rd quadrant
1016 den(14,3) = 0.4736
1017 vx(14,3) = 0.
1018 vy(14,3) = 1.2172
1019 pre(14,3) = 2.6667
1020 ! 4th quadrant
1021 den(14,4) = 0.9474
1022 vx(14,4) = 0.
1023 vy(14,4) = 1.1606
1024 pre(14,4) = 2.6667
1025
1026 ! Test configuration no. 15
1027 ! 1st quadrant
1028 den(15,1) = 1.
1029 vx(15,1) = 0.1
1030 vy(15,1) = -0.3
1031 pre(15,1) = 1.
1032 ! 2nd quadrant
1033 den(15,2) = 0.5197
1034 vx(15,2) = -0.6259
1035 vy(15,2) = -0.3
1036 pre(15,2) = 0.4
1037 ! 3rd quadrant
1038 den(15,3) = 0.8
1039 vx(15,3) = 0.1
1040 vy(15,3) = -0.3
1041 pre(15,3) = 0.4
1042 ! 4th quadrant
1043 den(15,4) = 0.5313
1044 vx(15,4) = 0.1
1045 vy(15,4) = 0.4276
1046 pre(15,4) = 0.4
1047
1048 ! Test configuration no. 16
1049 ! 1st quadrant
1050 den(16,1) = 0.5313
1051 vx(16,1) = 0.1
1052 vy(16,1) = 0.1
1053 pre(16,1) = 0.4
1054 ! 2nd quadrant
1055 den(16,2) = 1.02222
1056 vx(16,2) = -0.6179
1057 vy(16,2) = 0.1
1058 pre(16,2) = 1.
1059 ! 3rd quadrant
1060 den(16,3) = 0.8
1061 vx(16,3) = 0.1
1062 vy(16,3) = 0.1
1063 pre(16,3) = 1.
1064 ! 4th quadrant
1065 den(16,4) = 1.
1066 vx(16,4) = 0.1
1067 vy(16,4) = 0.8276
1068 pre(16,4) = 1.
1069
1070 ! Test configuration no. 17
1071 ! 1st quadrant
1072 den(17,1) = 1.
1073 vx(17,1) = 0.
1074 vy(17,1) = -0.4
1075 pre(17,1) = 1.
1076 ! 2nd quadrant
1077 den(17,2) = 2.
1078 vx(17,2) = 0.
1079 vy(17,2) = -0.3
1080 pre(17,2) = 1.
1081 ! 3rd quadrant
1082 den(17,3) = 1.0625
1083 vx(17,3) = 0.
1084 vy(17,3) = 0.2145
1085 pre(17,3) = 0.4
1086 ! 4th quadrant
1087 den(17,4) = 0.5197
1088 vx(17,4) = 0.
1089 vy(17,4) = -1.1259
1090 pre(17,4) = 0.4
1091
1092 ! Test configuration no. 18
1093 ! 1st quadrant
1094 den(18,1) = 1.
1095 vx(18,1) = 0.
1096 vy(18,1) = 1.
1097 pre(18,1) = 1.
1098 ! 2nd quadrant
1099 den(18,2) = 2.
1100 vx(18,2) = 0.
1101 vy(18,2) = -0.3
1102 pre(18,2) = 1.
1103 ! 3rd quadrant
1104 den(18,3) = 1.0625
1105 vx(18,3) = 0.
1106 vy(18,3) = 0.2145
1107 pre(18,3) = 0.4
1108 ! 4th quadrant
1109 den(18,4) = 0.5197
1110 vx(18,4) = 0.
1111 vy(18,4) = 0.2741
1112 pre(18,4) = 0.4
1113
1114 ! Test configuration no. 19
1115 ! 1st quadrant
1116 den(19,1) = 1.
1117 vx(19,1) = 0.
1118 vy(19,1) = 0.3
1119 pre(19,1) = 1.
1120 ! 2nd quadrant
1121 den(19,2) = 2.
1122 vx(19,2) = 0.
1123 vy(19,2) = -0.3
1124 pre(19,2) = 1.
1125 ! 3rd quadrant
1126 den(19,3) = 1.0625
1127 vx(19,3) = 0.
1128 vy(19,3) = 0.2145
1129 pre(19,3) = 0.4
1130 ! 4th quadrant
1131 den(19,4) = 0.5197
1132 vx(19,4) = 0.
1133 vy(19,4) = -0.4259
1134 pre(19,4) = 0.4
1135
1136 SELECT TYPE(pvar => timedisc%pvar)
1137 TYPE IS(statevector_euler) ! non-isothermal HD
1138 WHERE ( (mesh%bccart(:,:,:,1).GE.x0).AND.(mesh%bccart(:,:,:,2).GE.y0) )
1139 pvar%density%data3d(:,:,:) = den(ic,1)
1140 pvar%pressure%data3d(:,:,:) = pre(ic,1)
1141 vcart%data4d(:,:,:,1) = vx(ic,1)
1142 vcart%data4d(:,:,:,2) = vy(ic,1)
1143 ELSEWHERE ( (mesh%bccart(:,:,:,1).LT.x0).AND.(mesh%bccart(:,:,:,2).GE.y0) )
1144 pvar%density%data3d(:,:,:) = den(ic,2)
1145 pvar%pressure%data3d(:,:,:) = pre(ic,2)
1146 vcart%data4d(:,:,:,1) = vx(ic,2)
1147 vcart%data4d(:,:,:,2) = vy(ic,2)
1148 ELSEWHERE ( (mesh%bccart(:,:,:,1).LT.x0).AND.(mesh%bccart(:,:,:,2).LT.y0) )
1149 pvar%density%data3d(:,:,:) = den(ic,3)
1150 pvar%pressure%data3d(:,:,:) = pre(ic,3)
1151 vcart%data4d(:,:,:,1) = vx(ic,3)
1152 vcart%data4d(:,:,:,2) = vy(ic,3)
1153 ELSEWHERE ( (mesh%bccart(:,:,:,1).GE.x0).AND.(mesh%bccart(:,:,:,2).LT.y0) )
1154 pvar%density%data3d(:,:,:) = den(ic,4)
1155 pvar%pressure%data3d(:,:,:) = pre(ic,4)
1156 vcart%data4d(:,:,:,1) = vx(ic,4)
1157 vcart%data4d(:,:,:,2) = vy(ic,4)
1158 END WHERE
1159 vcart%data2d(:,3) = 0.0 ! z-velocity
1160 CALL mesh%geometry%Convert2Curvilinear(mesh%bcenter,vcart%data4d,vcurv%data4d)
1161 pvar%velocity%data2d(:,1:2) = vcurv%data2d(:,1:2)
1162 END SELECT
1163
1164 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
1165 CALL mesh%Info(" DATA-----> initial condition: " // trim(teststr))
1166
1167 END SUBROUTINE initdata
1168
1169END PROGRAM riemann2d
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
Mathematical functions.
Definition: functions.f90:48
elemental real function, public gamma(x)
Compute the Gamma function for arguments != 0,-1,-2,..
Definition: functions.f90:1133
real, parameter pi
Definition: functions.f90:52
elemental real function, public asinh(x)
inverse hyperbolic sine function
Definition: functions.f90:100
program riemann2d
Definition: riemann2d.f90:465
main fosite class
Definition: fosite.f90:71