diff --git a/demo/richards/richards.cfg b/demo/richards/richards.cfg index d251ddb6..46756e7f 100644 --- a/demo/richards/richards.cfg +++ b/demo/richards/richards.cfg @@ -15,25 +15,25 @@ standard_parallel= pressure = 1.0e-12 relative [richards-driver-snes-prob1] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES [richards-driver-ts-prob1] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-ts-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method TS +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-ts-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method TS [richards-driver-snes-prob1-np4] np=4 -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method SNES +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method SNES [richards-driver-ts-prob1-np4] np=4 timeout=300. -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-ts-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method TS +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-ts-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method TS [richards-driver-snes-initial-cond-isotropic-k] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 5,4,3 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-initial-cond-isotropic-k -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_time_integration_method SNES -init_permeability_file initial.h5 -ic_file initial.h5 +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 5,4,3 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-initial-cond-isotropic-k -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_time_integration_method SNES -init_permeability_file initial.h5 -ic_file initial.h5 [richards-driver-snes-initial-cond-anisotropic-k] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 5,4,3 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-initial-cond-anisotropic-k -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_time_integration_method SNES -init_porosity_file initial.h5 -init_permeability_file initial.h5 -anisotropic_perm -ic_file initial.h5 -ic_dataset fields/IC +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 5,4,3 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename richards-driver-snes-initial-cond-anisotropic-k -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_time_integration_method SNES -init_porosity_file initial.h5 -init_permeability_file initial.h5 -anisotropic_perm -ic_file initial.h5 -ic_dataset fields/IC [richards-driver-64xy-3z-wedge] input_arguments=-problem 4 -dm_plex_simplex 0 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-64xy-3z-wedge -dm_plex_filename ../../share/meshes/64xy_3z_wedge.exo @@ -48,7 +48,7 @@ input_arguments=-snes_linesearch_basic -ts_dt 100 -ts_type beuler -ts_max_time 1 input_arguments=-snes_linesearch_type basic -ts_dt 100 -ts_type bdf -ts_adapt_type none -ts_max_time 1000 -ts_max_steps 10 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-transientvar -tdy_method MPFA_O_TRANSIENTVAR -pc_type lu [richards-driver-snes-checkpoint-write] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-checkpoint-write -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES -enable_checkpoint +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-checkpoint-write -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES -enable_checkpoint [richards-driver-snes-checkpoint-read] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-checkpoint-read -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES -ic_file 3.15360e+03_chk.h5 -ic_dataset fields/IC \ No newline at end of file +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename richards-driver-snes-checkpoint-read -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_file richards_driver_snes_prob1_3x3x3_random.vec -tdy_time_integration_method SNES -ic_file 3.15360e+03_chk.h5 -ic_dataset fields/IC diff --git a/demo/richards/richards_driver.c b/demo/richards/richards_driver.c index 35c2fa98..67efcbc3 100644 --- a/demo/richards/richards_driver.c +++ b/demo/richards/richards_driver.c @@ -11,14 +11,15 @@ int main(int argc, char **argv) { TDyIOFormat format = HDF5Format; ierr = TDyInit(argc, argv); CHKERRQ(ierr); - ierr = TDyCreate(&tdy); CHKERRQ(ierr); + MPI_Comm comm = PETSC_COMM_WORLD; + ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); ierr = TDySetMode(tdy,RICHARDS); CHKERRQ(ierr); - ierr = TDySetDiscretizationMethod(tdy,MPFA_O); CHKERRQ(ierr); + ierr = TDySetDiscretization(tdy,MPFA_O); CHKERRQ(ierr); - ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr); - ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size); CHKERRQ(ierr); - PetscPrintf(PETSC_COMM_WORLD,"Beginning Richards Driver simulation.\n"); - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Sample Options",""); + ierr = MPI_Comm_rank(comm,&rank); CHKERRQ(ierr); + ierr = MPI_Comm_size(comm,&size); CHKERRQ(ierr); + PetscPrintf(comm,"Beginning Richards Driver simulation.\n"); + ierr = PetscOptionsBegin(comm,NULL,"Sample Options",""); CHKERRQ(ierr); ierr = PetscOptionsInt("-successful_exit_code", "Code passed on successful completion","", @@ -33,7 +34,7 @@ int main(int argc, char **argv) { if (!rank) { ierr = TDyIOSetIOProcess(tdy->io, PETSC_TRUE); CHKERRQ(ierr); } - PetscPrintf(PETSC_COMM_WORLD,"--\n"); + PetscPrintf(comm,"--\n"); if (size == 1) { ierr = TDyIOSetMode(tdy,format);CHKERRQ(ierr); ierr = TDyIOWriteVec(tdy); CHKERRQ(ierr); @@ -44,8 +45,8 @@ int main(int argc, char **argv) { ierr = TDyOutputRegression(tdy,tdy->solution); CHKERRQ(ierr); ierr = TDyDestroy(&tdy); CHKERRQ(ierr); - PetscPrintf(PETSC_COMM_WORLD,"--\n"); - PetscPrintf(PETSC_COMM_WORLD,"Simulation complete.\n"); + PetscPrintf(comm,"--\n"); + PetscPrintf(comm,"Simulation complete.\n"); ierr = TDyFinalize(); CHKERRQ(ierr); return(successful_exit_code); } diff --git a/demo/steady/steady.c b/demo/steady/steady.c index f67a4bea..fab4c1c4 100644 --- a/demo/steady/steady.c +++ b/demo/steady/steady.c @@ -264,7 +264,6 @@ PetscErrorCode Forcing3(TDy tdy,PetscReal *x,PetscReal *f,void *ctx) { PetscFunctionReturn(0); } - PetscErrorCode PerturbVerticesRandom(DM dm,PetscReal h) { PetscErrorCode ierr; PetscFunctionBegin; @@ -274,8 +273,7 @@ PetscErrorCode PerturbVerticesRandom(DM dm,PetscReal h) { PetscScalar *coords; PetscInt v,vStart,vEnd,offset,value,dim; ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - /* this is the 'marker' label which marks boundary entities */ - ierr = DMGetLabelByNum(dm,2,&label); CHKERRQ(ierr); + ierr = DMGetLabel(dm, "boundary", &label); CHKERRQ(ierr); ierr = DMGetCoordinateSection(dm, &coordSection); CHKERRQ(ierr); ierr = DMGetCoordinatesLocal(dm, &coordinates); CHKERRQ(ierr); ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd); CHKERRQ(ierr); @@ -292,7 +290,7 @@ PetscErrorCode PerturbVerticesRandom(DM dm,PetscReal h) { coords[offset+1] += r*PetscSinReal(t); } } else { - /* This is because 'marker' is broken in 3D */ + /* This is because 'boundary' is broken in 3D */ if(coords[offset] > 0 && coords[offset] < 1 && coords[offset+1] > 0 && coords[offset+1] < 1 && coords[offset+2] > 0 && coords[offset+2] < 1) { @@ -303,6 +301,7 @@ PetscErrorCode PerturbVerticesRandom(DM dm,PetscReal h) { ierr = VecRestoreArray(coordinates,&coords); CHKERRQ(ierr); PetscFunctionReturn(0); } + PetscErrorCode PerturbVerticesSmooth(DM dm) { PetscErrorCode ierr; PetscFunctionBegin; @@ -393,7 +392,9 @@ PetscErrorCode SaveVertices(DM dm, char filename[256]){ Vec coordinates; PetscViewer viewer; - ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer); CHKERRQ(ierr); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_WRITE,&viewer); CHKERRQ(ierr); ierr = DMGetCoordinatesLocal(dm, &coordinates); CHKERRQ(ierr); ierr = VecView(coordinates,viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); @@ -414,7 +415,9 @@ PetscErrorCode SaveCentroids(DM dm, char filename[256]){ ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr); - ierr = VecCreateMPI(PETSC_COMM_WORLD,(cEnd-cStart)*dim,PETSC_DETERMINE,¢roids); CHKERRQ(ierr); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + ierr = VecCreateMPI(comm,(cEnd-cStart)*dim,PETSC_DETERMINE,¢roids); CHKERRQ(ierr); ierr = VecGetArray(centroids,¢roids_p); CHKERRQ(ierr); for(c=cStart; cN; + if(options->exo) { + ierr = DMPlexCreateExodusFromFile(PETSC_COMM_WORLD, options->exofile, + PETSC_TRUE,dm); CHKERRQ(ierr); + } else { + PetscInt Nx=N,Ny=N,Nz=N; + PetscReal Lx=1,Ly=1,Lz=1; + if(options->column){ + Nx = 1; Ny = 1; Nz = N; + Lx = 10; Ly = 10; Lz = 1; + } + const PetscInt faces[3] = {Nx ,Ny ,Nz }; + const PetscReal lower[3] = {0.0,0.0,0.0}; + const PetscReal upper[3] = {Lx ,Ly ,Lz }; + ierr = DMPlexCreateBoxMesh(comm, options->dim, PETSC_FALSE, + faces,lower,upper,NULL,PETSC_TRUE,dm); CHKERRQ(ierr); + } + PetscFunctionReturn(0); +} + int main(int argc, char **argv) { /* Initialize */ PetscErrorCode ierr; - PetscInt N = 4, dim = 2, problem = 2; PetscInt successful_exit_code=0; - PetscBool perturb = PETSC_FALSE; - char exofile[256],paper[256]="none"; char vertices_filename[256]="none"; char centroids_filename[256]="none"; char true_pres_filename[256]="none"; char forcing_filename[256]="none"; - PetscBool exo = PETSC_FALSE; + char paper[256]; + PetscBool wheeler2006, wheeler2012; + + // DM-related options. + PetscInt N = 4, dim = 2, problem = 2; + PetscBool perturb = PETSC_FALSE, exo = PETSC_FALSE, column; + char exofile[256]; ierr = TDyInit(argc, argv); CHKERRQ(ierr); - TDy tdy; - ierr = TDyCreate(&tdy); CHKERRQ(ierr); - ierr = TDySetDiscretizationMethod(tdy,MPFA_O); CHKERRQ(ierr); + MPI_Comm comm = PETSC_COMM_WORLD; + TDy tdy; + ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); + ierr = TDySetMode(tdy,RICHARDS); CHKERRQ(ierr); + ierr = TDySetDiscretization(tdy,MPFA_O); CHKERRQ(ierr); - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Sample Options",""); CHKERRQ(ierr); + ierr = PetscOptionsBegin(comm,NULL,"Sample Options",""); CHKERRQ(ierr); ierr = PetscOptionsInt("-dim","Problem dimension","",dim,&dim,NULL); CHKERRQ(ierr); ierr = PetscOptionsInt("-N","Number of elements in 1D","",N,&N,NULL); CHKERRQ(ierr); ierr = PetscOptionsInt("-problem","Problem number","",problem,&problem,NULL); CHKERRQ(ierr); - ierr = PetscOptionsBool("-perturb","Perturb interior vertices","",perturb,&perturb,NULL); CHKERRQ(ierr); ierr = PetscOptionsReal("-alpha","Permeability scaling","",alpha,&alpha,NULL); CHKERRQ(ierr); ierr = PetscOptionsInt("-successful_exit_code","Code passed on successful completion","",successful_exit_code,&successful_exit_code,NULL); ierr = PetscOptionsString("-exo","Mesh file in exodus format","",exofile,exofile,256,&exo); CHKERRQ(ierr); @@ -580,59 +626,22 @@ int main(int argc, char **argv) { ierr = PetscOptionsString("-view_forcing","Filename to save forcing","",forcing_filename,forcing_filename,256,NULL); CHKERRQ(ierr); ierr = PetscOptionsEnd(); CHKERRQ(ierr); - PetscBool wheeler2006,wheeler2012,column; ierr = PetscStrcasecmp(paper,"wheeler2006",&wheeler2006); CHKERRQ(ierr); ierr = PetscStrcasecmp(paper,"wheeler2012",&wheeler2012); CHKERRQ(ierr); ierr = PetscStrcasecmp(paper,"column",&column); CHKERRQ(ierr); - /* Create and distribute the mesh */ - DM dm, dmDist = NULL; - DMLabel marker; - if(exo){ - ierr = DMPlexCreateExodusFromFile(PETSC_COMM_WORLD,exofile, - PETSC_TRUE,&dm); CHKERRQ(ierr); - ierr = DMSetFromOptions(dm); CHKERRQ(ierr); - ierr = DMCreateLabel(dm,"marker"); CHKERRQ(ierr); - ierr = DMGetLabel(dm,"marker",&marker); CHKERRQ(ierr); - ierr = DMPlexMarkBoundaryFaces(dm,1,marker); CHKERRQ(ierr); - }else{ - PetscInt Nx=N,Ny=N,Nz=N; - PetscReal Lx=1,Ly=1,Lz=1; - if(column){ - Nx = 1; Ny = 1; Nz = N; - Lx = 10; Ly = 10; Lz = 1; - } - const PetscInt faces[3] = {Nx ,Ny ,Nz }; - const PetscReal lower[3] = {0.0,0.0,0.0}; - const PetscReal upper[3] = {Lx ,Ly ,Lz }; - ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,dim,PETSC_FALSE,faces,lower,upper, - NULL,PETSC_TRUE,&dm); CHKERRQ(ierr); - if(wheeler2006 && (problem==2 || problem==0)){ - ierr = GeometryWheeler2006_2(dm,N); CHKERRQ(ierr); - } - if(column) { - ierr = GeometryColumn(dm); CHKERRQ(ierr); - } - if(perturb) { - ierr = PerturbVerticesRandom(dm,1./N); CHKERRQ(ierr); - }else{ - if(wheeler2012) { - ierr = PerturbVerticesSmooth(dm); CHKERRQ(ierr); - } - } - } - ierr = DMPlexDistribute(dm, 1, NULL, &dmDist); - if (dmDist) {DMDestroy(&dm); dm = dmDist;} - ierr = DMSetFromOptions(dm); CHKERRQ(ierr); - ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); + // Specify a special DM to be constructed for this demo, and pass it the + // relevant options. + DMOptions dm_options = {.N = N, .dim = dim, .exo = exo, .exofile = exofile}; + ierr = TDySetDMConstructor(tdy, &dm_options, CreateDM); CHKERRQ(ierr); - ierr = TDySetDM(tdy,dm); CHKERRQ(ierr); + // Apply overrides. ierr = TDySetFromOptions(tdy); CHKERRQ(ierr); - /* Setup problem parameters */ + // Setup problem parameters if(wheeler2006){ if(dim != 2){ - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"-paper wheeler2006 is only for -dim 2 problems"); + SETERRQ(comm,PETSC_ERR_USER,"-paper wheeler2006 is only for -dim 2 problems"); } switch(problem) { case 0: // not a problem in the paper, but want to check constants on the geometry @@ -654,13 +663,13 @@ int main(int argc, char **argv) { ierr = TDySetBoundaryVelocityFn(tdy,VelocityWheeler2006_2,NULL); CHKERRQ(ierr); break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"-paper wheeler2006 only valid for -problem {0,1,2}"); + SETERRQ(comm,PETSC_ERR_USER,"-paper wheeler2006 only valid for -problem {0,1,2}"); } }else if(wheeler2012){ switch(problem) { case 1: if(dim != 2){ - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"-paper wheeler2012 -problem 1 is only for -dim 2"); + SETERRQ(comm,PETSC_ERR_USER,"-paper wheeler2012 -problem 1 is only for -dim 2"); } ierr = TDySetPermeabilityTensor(tdy,PermWheeler2012_1); CHKERRQ(ierr); ierr = TDySetForcingFunction(tdy,ForcingWheeler2012_1,NULL); CHKERRQ(ierr); @@ -669,7 +678,7 @@ int main(int argc, char **argv) { break; case 2: if(dim != 3){ - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"-paper wheeler2012 -problem 2 is only for -dim 3"); + SETERRQ(comm,PETSC_ERR_USER,"-paper wheeler2012 -problem 2 is only for -dim 3"); } ierr = TDySetPermeabilityTensor(tdy,PermWheeler2012_2); CHKERRQ(ierr); ierr = TDySetForcingFunction(tdy,ForcingWheeler2012_2,NULL); CHKERRQ(ierr); @@ -677,7 +686,7 @@ int main(int argc, char **argv) { ierr = TDySetBoundaryVelocityFn(tdy,VelocityWheeler2012_2,NULL); CHKERRQ(ierr); break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"-paper wheeler2012 only valid for -problem {1,2}"); + SETERRQ(comm,PETSC_ERR_USER,"-paper wheeler2012 only valid for -problem {1,2}"); } }else{ switch(problem) { @@ -736,6 +745,26 @@ int main(int argc, char **argv) { ierr = TDySetup(tdy); CHKERRQ(ierr); + // Make adjustments to our DM based on the problem. + DM dm; + TDyGetDM(tdy, &dm); + if(wheeler2006 && (problem==2 || problem==0)){ + ierr = GeometryWheeler2006_2(dm,N); CHKERRQ(ierr); + } + if(column) { + ierr = GeometryColumn(dm); CHKERRQ(ierr); + } + if(perturb) { + ierr = PerturbVerticesRandom(dm,1./N); CHKERRQ(ierr); + }else{ + if(wheeler2012) { + ierr = PerturbVerticesSmooth(dm); CHKERRQ(ierr); + } + } + + // View the configured DM. + ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); + /* Compute system */ Mat K; Vec U,Ue,F; @@ -747,7 +776,7 @@ int main(int argc, char **argv) { /* Solve system */ KSP ksp; - ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr); + ierr = KSPCreate(comm,&ksp); CHKERRQ(ierr); ierr = KSPSetOperators(ksp,K,K); CHKERRQ(ierr); ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); ierr = KSPSetUp(ksp); CHKERRQ(ierr); @@ -766,7 +795,7 @@ int main(int argc, char **argv) { /* Evaluate error norms */ PetscReal normp,normv; ierr = TDyComputeErrorNorms(tdy,U,&normp,&normv); - ierr = PetscPrintf(PETSC_COMM_WORLD,"%e %e\n",normp,normv); CHKERRQ(ierr); + ierr = PetscPrintf(comm,"%e %e\n",normp,normv); CHKERRQ(ierr); /* Save vertex coordinates */ PetscBool file_not_specified; diff --git a/demo/steady/steady.cfg b/demo/steady/steady.cfg index c203d220..a463e9df 100644 --- a/demo/steady/steady.cfg +++ b/demo/steady/steady.cfg @@ -7,10 +7,6 @@ standard= steady-wy-prob1-3D steady-wy-prob2-3D steady-wy-prob3-3D - steady-tpf-prob1 - steady-tpf-prob2 - steady-tpf-prob3 - steady-tpf-prob4 steady-bdm-prob1 steady-bdm-prob2 steady-bdm-prob3 @@ -20,48 +16,36 @@ standard= pressure = 1.0e-12 relative [steady-wy-prob1] -input_arguments=-tdy_method wy -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1 +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1 [steady-wy-prob2] -input_arguments=-tdy_method wy -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2 +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2 [steady-wy-prob3] -input_arguments=-tdy_method wy -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3 +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3 [steady-wy-prob4] -input_arguments=-tdy_method wy -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob4 +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob4 [steady-wy-prob1-3D] -input_arguments=-tdy_method wy -tdy_regression_test -problem 1 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1-3D +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 1 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob1-3D [steady-wy-prob2-3D] -input_arguments=-tdy_method wy -tdy_regression_test -problem 2 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2-3D +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 2 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob2-3D [steady-wy-prob3-3D] -input_arguments=-tdy_method wy -tdy_regression_test -problem 3 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3-3D - -[steady-tpf-prob1] -input_arguments=-tdy_method tpf -tdy_regression_test -problem 1 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob1 - -[steady-tpf-prob2] -input_arguments=-tdy_method tpf -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob2 - -[steady-tpf-prob3] -input_arguments=-tdy_method tpf -tdy_regression_test -problem 3 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob3 - -[steady-tpf-prob4] -input_arguments=-tdy_method tpf -tdy_regression_test -problem 4 -N 8 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-tpf-prob4 +input_arguments=-tdy_discretization wy -tdy_regression_test -problem 3 -dim 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-wy-prob3-3D [steady-bdm-prob1] pressure = 1.0e-12 absolute -input_arguments=-tdy_method bdm -tdy_regression_test -problem 1 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob1 -ksp_type preonly -pc_type lu +input_arguments=-tdy_discretization bdm -tdy_regression_test -problem 1 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob1 -ksp_type preonly -pc_type lu [steady-bdm-prob2] -input_arguments=-tdy_method bdm -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob2 -ksp_type preonly -pc_type lu +input_arguments=-tdy_discretization bdm -tdy_regression_test -problem 2 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob2 -ksp_type preonly -pc_type lu [steady-bdm-prob3] -input_arguments=-tdy_method bdm -tdy_regression_test -problem 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob3 -ksp_type preonly -pc_type lu +input_arguments=-tdy_discretization bdm -tdy_regression_test -problem 3 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob3 -ksp_type preonly -pc_type lu [steady-bdm-prob4] -input_arguments=-tdy_method bdm -tdy_regression_test -problem 4 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob4 -ksp_type preonly -pc_type lu +input_arguments=-tdy_discretization bdm -tdy_regression_test -problem 4 -N 4 -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename steady-bdm-prob4 -ksp_type preonly -pc_type lu diff --git a/demo/th/th.cfg b/demo/th/th.cfg index de037b6a..98c3f2a0 100644 --- a/demo/th/th.cfg +++ b/demo/th/th.cfg @@ -9,17 +9,17 @@ standard_parallel= pressure = 1.0e-12 relative #[th-driver-snes-prob1] -#input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename th-driver-snes-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method SNES +#input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename th-driver-snes-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method SNES [th-driver-ts-prob1] -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename th-driver-ts-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method TS +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename th-driver-ts-prob1 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method TS #[th-driver-snes-prob1-np4] #np=4 -#input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename th-driver-snes-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method SNES +#input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename th-driver-snes-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_init_with_random_field -tdy_time_integration_method SNES [th-driver-ts-prob1-np4] np=4 timeout=300. -input_arguments=-tdy_generate_mesh -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename th-driver-ts-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method TS +input_arguments=-dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1,1,1 -tdy_water_density exponential -tdy_regression_test -tdy_regression_test_num_cells_per_process 1 -tdy_regression_test_filename th-driver-ts-prob1-np4 -tdy_final_time 3.1536e3 -tdy_dt_max 600. -tdy_dt_growth_factor 1.5 -tdy_timers -tdy_init_with_random_field -tdy_time_integration_method TS diff --git a/demo/th/th_driver.c b/demo/th/th_driver.c index 239232b0..5b4b271d 100644 --- a/demo/th/th_driver.c +++ b/demo/th/th_driver.c @@ -8,17 +8,18 @@ int main(int argc, char **argv) { PetscInt successful_exit_code; PetscMPIInt rank, size; TDy tdy = PETSC_NULL; - TDyIOFormat format = PetscViewerASCIIFormat; + TDyIOFormat format = PetscViewerASCIIFormat; ierr = TDyInit(argc, argv); CHKERRQ(ierr); - ierr = TDyCreate(&tdy); CHKERRQ(ierr); + MPI_Comm comm = PETSC_COMM_WORLD; + ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); ierr = TDySetMode(tdy,TH); CHKERRQ(ierr); - ierr = TDySetDiscretizationMethod(tdy,MPFA_O); CHKERRQ(ierr); + ierr = TDySetDiscretization(tdy,MPFA_O); CHKERRQ(ierr); - ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr); - ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size); CHKERRQ(ierr); - PetscPrintf(PETSC_COMM_WORLD,"Beginning TH Driver simulation.\n"); - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Sample Options",""); + ierr = MPI_Comm_rank(comm,&rank); CHKERRQ(ierr); + ierr = MPI_Comm_size(comm,&size); CHKERRQ(ierr); + PetscPrintf(comm,"Beginning TH Driver simulation.\n"); + ierr = PetscOptionsBegin(comm,NULL,"Sample Options",""); CHKERRQ(ierr); ierr = PetscOptionsInt("-successful_exit_code", "Code passed on successful completion","", @@ -33,19 +34,19 @@ int main(int argc, char **argv) { if (!rank) { ierr = TDyIOSetIOProcess(tdy->io, PETSC_TRUE); CHKERRQ(ierr); } - PetscPrintf(PETSC_COMM_WORLD,"--\n"); + PetscPrintf(comm,"--\n"); if (size == 1) { ierr = TDyIOSetMode(tdy,format);CHKERRQ(ierr); ierr = TDyIOWriteVec(tdy); CHKERRQ(ierr); } - ierr = TDyTimeIntegratorRunToTime(tdy,tdy->ti->final_time); + ierr = TDyTimeIntegratorRunToTime(tdy,tdy->ti->final_time); CHKERRQ(ierr); if (size == 1) {ierr = TDyIOWriteVec(tdy); CHKERRQ(ierr);} ierr = TDyOutputRegression(tdy,tdy->solution); CHKERRQ(ierr); ierr = TDyDestroy(&tdy); CHKERRQ(ierr); - PetscPrintf(PETSC_COMM_WORLD,"--\n"); - PetscPrintf(PETSC_COMM_WORLD,"Simulation complete.\n"); + PetscPrintf(comm,"--\n"); + PetscPrintf(comm,"Simulation complete.\n"); ierr = TDyFinalize(); CHKERRQ(ierr); return(successful_exit_code); } diff --git a/demo/transient/transient.c b/demo/transient/transient.c index ac70aa4c..b0c7febb 100644 --- a/demo/transient/transient.c +++ b/demo/transient/transient.c @@ -58,6 +58,34 @@ PetscErrorCode PerturbInteriorVertices(DM dm,PetscReal h) { PetscFunctionReturn(0); } +// This data is used by CreateDM below to create a DM for this demo. +typedef struct DMOptions { + PetscInt dim; // Dimension of DM (2 or 3) + PetscInt N; // Number of cells on a side + PetscBool exo; // whether to load a named exodus file + const char* exofile; // name of the exodus file to load +} DMOptions; + +// This function creates a DM specifically for this demo. Overrides are applied +// to the resulting DM with TDySetFromOptions. +PetscErrorCode CreateDM(MPI_Comm comm, void* context, DM* dm) { + int ierr; + DMOptions* options = context; + + PetscInt N = options->N; + if(options->exo){ + ierr = DMPlexCreateExodusFromFile(PETSC_COMM_WORLD,options->exofile, + PETSC_TRUE,dm); CHKERRQ(ierr); + } else { + const PetscInt faces[3] = {N,N,N }; + const PetscReal lower[3] = {-0.5,-0.5,-0.5}; + const PetscReal upper[3] = {+0.5,+0.5,+0.5}; + ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,options->dim,PETSC_FALSE, + faces,lower,upper,NULL,PETSC_TRUE,dm); CHKERRQ(ierr); + ierr = PerturbInteriorVertices(*dm,1./N); CHKERRQ(ierr); + } +} + int main(int argc, char **argv) { /* Initialize */ PetscErrorCode ierr; @@ -67,9 +95,11 @@ int main(int argc, char **argv) { PetscBool exo = PETSC_FALSE; ierr = TDyInit(argc, argv); CHKERRQ(ierr); + MPI_Comm comm = PETSC_COMM_WORLD; TDy tdy; - ierr = TDyCreate(&tdy); CHKERRQ(ierr); - ierr = TDySetDiscretizationMethod(tdy,WY); CHKERRQ(ierr); + ierr = TDyCreate(comm, &tdy); CHKERRQ(ierr); + ierr = TDySetMode(tdy, RICHARDS); CHKERRQ(ierr); + ierr = TDySetDiscretization(tdy,WY); CHKERRQ(ierr); ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL, "Transient Options",""); CHKERRQ(ierr); @@ -81,33 +111,19 @@ int main(int argc, char **argv) { exofile,exofile,256,&exo); CHKERRQ(ierr); ierr = PetscOptionsEnd(); CHKERRQ(ierr); - /* Create and distribute the mesh */ - DM dm, dmDist = NULL; - DMLabel marker; - if(exo){ - ierr = DMPlexCreateExodusFromFile(PETSC_COMM_WORLD,exofile, - PETSC_TRUE,&dm); CHKERRQ(ierr); - //ierr = DMPlexOrient(dm); CHKERRQ(ierr); - ierr = DMSetFromOptions(dm); CHKERRQ(ierr); - ierr = DMCreateLabel(dm,"marker"); CHKERRQ(ierr); - ierr = DMGetLabel(dm,"marker",&marker); CHKERRQ(ierr); - ierr = DMPlexMarkBoundaryFaces(dm,1,marker); CHKERRQ(ierr); - }else{ - const PetscInt faces[3] = {N,N,N }; - const PetscReal lower[3] = {-0.5,-0.5,-0.5}; - const PetscReal upper[3] = {+0.5,+0.5,+0.5}; - ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,dim,PETSC_FALSE,faces,lower,upper, - NULL,PETSC_TRUE,&dm); CHKERRQ(ierr); - ierr = PerturbInteriorVertices(dm,1./N); CHKERRQ(ierr); - } - ierr = DMPlexDistribute(dm, 1, NULL, &dmDist); - if (dmDist) {DMDestroy(&dm); dm = dmDist;} - ierr = DMSetFromOptions(dm); CHKERRQ(ierr); - ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); + // Specify a special DM to be constructed for this demo, and pass it the + // relevant options. + DMOptions dm_options = {.N = N, .dim = dim, .exo = exo, .exofile = exofile}; + ierr = TDySetDMConstructor(tdy, &dm_options, CreateDM); CHKERRQ(ierr); - ierr = TDySetDM(tdy,dm); CHKERRQ(ierr); + // Apply overrides. ierr = TDySetFromOptions(tdy); CHKERRQ(ierr); + // View the configured DM. + DM dm; + TDyGetDM(tdy, &dm); + ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); + /* Setup problem parameters */ ierr = TDySetPorosity(tdy,Porosity); CHKERRQ(ierr); ierr = TDySetPermeabilityScalar(tdy,Permeability); CHKERRQ(ierr); diff --git a/demo/transient/transient_mpfaof90.F90 b/demo/transient/transient_mpfaof90.F90 index 7adf2790..343c9b99 100644 --- a/demo/transient/transient_mpfaof90.F90 +++ b/demo/transient/transient_mpfaof90.F90 @@ -4,6 +4,9 @@ module mpfaof90mod #include #include + ! Module variables + PetscInt :: dim = 3 + contains subroutine PorosityFunction(tdy,x,theta,dummy,ierr) @@ -53,6 +56,26 @@ subroutine PressureFunction(tdy,x,pressure,dummy,ierr) ierr = 0 end subroutine PressureFunction + subroutine CreateDM(comm, dm, ierr) + use petscdm + implicit none + + MPI_Comm :: comm + DM :: dm + PetscErrorCode :: ierr + + PetscInt :: faces(3), nx, ny, nz + PetscReal :: lower(3), upper(3) + + nx = 3; ny = 3; nz = 3; + faces(1) = nx; faces(2) = ny; faces(3) = nz; + lower(:) = 0.d0; + upper(:) = 1.d0; + + call DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, faces, lower, upper, & + PETSC_NULL_INTEGER, PETSC_TRUE, dm, ierr); + CHKERRA(ierr) + end subroutine CreateDM end module mpfaof90mod program main @@ -71,15 +94,12 @@ program main implicit none TDy :: tdy - DM :: dm, dmDist + DM :: dm Vec :: U TS :: ts PetscInt :: rank, successful_exit_code PetscBool :: flg - PetscInt :: dim, faces(3) - PetscReal :: lower(3), upper(3) PetscErrorCode :: ierr - PetscInt :: nx, ny, nz PetscInt, pointer :: index(:) PetscReal, pointer :: residualSat(:), blockPerm(:), liquid_sat(:), liquid_mass(:) PetscReal, pointer :: p_loc(:) @@ -90,11 +110,11 @@ program main CHKERRA(ierr); call TDyCreate(tdy, ierr); CHKERRA(ierr); - call TDySetDiscretizationMethod(tdy,MPFA_O,ierr); + call TDySetMode(tdy,RICHARDS,ierr); + CHKERRA(ierr); + call TDySetDiscretization(tdy,MPFA_O,ierr); CHKERRA(ierr); - nx = 3; ny = 3; nz = 3; - dim = 3 successful_exit_code= 0 call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-successful_exit_code',successful_exit_code,flg,ierr); @@ -102,23 +122,14 @@ program main call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr); CHKERRA(ierr) - faces(1) = nx; faces(2) = ny; faces(3) = nz; - lower(:) = 0.d0; - upper(:) = 1.d0; - - call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, faces, lower, upper, & - PETSC_NULL_INTEGER, PETSC_TRUE, dm, ierr); + ! Set a constructor for a DM. + call TDySetDMConstructor(tdy, CreateDM,ierr); CHKERRA(ierr); - call DMPlexDistribute(dm, 1, PETSC_NULL_SF, dmDist, ierr); + + call TDySetFromOptions(tdy,ierr); CHKERRA(ierr); - if (dmDist /= PETSC_NULL_DM) then - call DMDestroy(dm, ierr); - CHKERRA(ierr); - dm = dmDist; - end if - call DMSetUp(dm,ierr); - CHKERRA(ierr) + call TDyGetDM(tdy, dm, ierr); call DMPlexGetHeightStratum(dm,0,cStart,cEnd,ierr); CHKERRA(ierr); @@ -138,14 +149,6 @@ program main enddo enddo - call DMSetFromOptions(dm, ierr); - CHKERRA(ierr); - - call TDySetDM(tdy, dm, ierr); - CHKERRA(ierr); - call TDySetFromOptions(tdy,ierr); - CHKERRA(ierr); - call TDySetPorosityFunction(tdy,PorosityFunction,0,ierr); CHKERRA(ierr); diff --git a/demo/transient/transient_snes_mpfaof90.F90 b/demo/transient/transient_snes_mpfaof90.F90 index 14113b83..6d74f0f3 100644 --- a/demo/transient/transient_snes_mpfaof90.F90 +++ b/demo/transient/transient_snes_mpfaof90.F90 @@ -4,6 +4,13 @@ module snes_mpfaof90mod #include #include + ! DM-related variables for CreateDM. + character (len=256) :: mesh_filename + PetscBool :: mesh_file_flg + PetscInt :: nx, ny, nz + PetscInt :: dim + PetscInt :: dm_plex_extrude_layers + contains subroutine PorosityFunction(tdy,x,theta,dummy,ierr) @@ -113,6 +120,40 @@ subroutine PressureFunction(tdy,x,pressure,dummy,ierr) ierr = 0 end subroutine PressureFunction + subroutine CreateDM(comm, dm, ierr) + use petscdm + implicit none + + MPI_Comm :: comm + DM :: dm, edm + PetscErrorCode :: ierr + + PetscInt :: faces(3) + PetscReal :: lower(3), upper(3) + + if (.not.mesh_file_flg) then + faces(1) = nx; faces(2) = ny; faces(3) = nz; + lower(:) = 0.d0; + upper(:) = 1.d0; + + call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, faces, lower, upper, & + PETSC_NULL_INTEGER, PETSC_TRUE, dm, ierr); + CHKERRA(ierr); + else + call DMPlexCreateFromFile(PETSC_COMM_WORLD, mesh_filename, PETSC_TRUE, dm, ierr); + CHKERRA(ierr); + call DMGetDimension(dm, dim, ierr); + CHKERRA(ierr); + if (dm_plex_extrude_layers > 0) then + call DMPlexExtrude(dm, PETSC_DETERMINE, -1.d0, PETSC_TRUE, PETSC_NULL_REAL, PETSC_TRUE, edm, ierr); + CHKERRA(ierr); + call DMDestroy(dm ,ierr); + dm = edm + end if + endif + + end subroutine + end module snes_mpfaof90mod program main @@ -134,30 +175,31 @@ program main implicit none TDy :: tdy - DM :: dm, edm, dmDist + DM :: dm Vec :: U !TS :: ts SNES :: snes PetscInt :: rank, successful_exit_code PetscBool :: flg - PetscInt :: dim, faces(3) - PetscReal :: lower(3), upper(3) PetscErrorCode :: ierr - PetscInt :: nx, ny, nz, ncell, bc_type, dm_plex_extrude_layers + PetscInt :: ncell, bc_type PetscInt , pointer :: index(:) PetscReal , pointer :: residualSat(:), blockPerm(:), liquid_sat(:), liquid_mass(:) PetscReal , pointer :: alpha(:), m(:) PetscReal :: perm(9), resSat PetscInt :: c, cStart, cEnd, j, nvalues,g, max_steps, step PetscReal :: dtime, mass_pre, mass_post, ic_value - character (len=256) :: mesh_filename, ic_filename + character (len=256) :: ic_filename character(len=256) :: string, bc_type_name - PetscBool :: mesh_file_flg, ic_file_flg, pflotran_consistent, use_tdydriver + PetscBool :: ic_file_flg, pflotran_consistent, use_tdydriver PetscViewer :: viewer PetscInt :: step_mod PetscFE :: fe SNESConvergedReason :: reason + dim = 3 + nx = 1; ny = 1; nz = 15; + call TDyInit(ierr); CHKERRA(ierr); @@ -169,11 +211,11 @@ program main call TDyCreate(tdy, ierr); CHKERRA(ierr); - call TDySetDiscretizationMethod(tdy,MPFA_O,ierr); + call TDySetMode(tdy,RICHARDS,ierr); + CHKERRA(ierr); + call TDySetDiscretization(tdy,MPFA_O,ierr); CHKERRA(ierr); - nx = 1; ny = 1; nz = 15; - dim = 3 successful_exit_code= 0 max_steps = 2 dtime = 1800.d0 @@ -224,27 +266,17 @@ program main bc_type = MPFAO_SEEPAGE_BC endif - if (.not.mesh_file_flg) then - faces(1) = nx; faces(2) = ny; faces(3) = nz; - lower(:) = 0.d0; - upper(:) = 1.d0; + ! Set a constructor for a DM. + call TDySetDMConstructor(tdy, CreateDM,ierr); + CHKERRA(ierr); - call DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, faces, lower, upper, & - PETSC_NULL_INTEGER, PETSC_TRUE, dm, ierr); - CHKERRA(ierr); - else - call DMPlexCreateFromFile(PETSC_COMM_WORLD, mesh_filename, PETSC_TRUE, dm, ierr); - CHKERRA(ierr); - call DMGetDimension(dm, dim, ierr); - CHKERRA(ierr); - if (dm_plex_extrude_layers > 0) then - call DMPlexExtrude(dm, PETSC_DETERMINE, -1.d0, PETSC_TRUE, PETSC_NULL_REAL, PETSC_TRUE, edm, ierr); - CHKERRA(ierr); - call DMDestroy(dm ,ierr); - dm = edm - end if - endif + ! Apply overrides. + call TDySetFromOptions(tdy,ierr); + CHKERRA(ierr); + ! Set up the discretization. + call TDyGetDM(tdy, dm, ierr); + CHKERRA(ierr); call DMGetDimension(dm, dim, ierr); CHKERRA(ierr) call PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, "p_", -1, fe, ierr); @@ -259,25 +291,6 @@ program main CHKERRA(ierr) call DMSetUseNatural(dm, PETSC_TRUE, ierr); CHKERRA(ierr) - - call DMPlexDistribute(dm, 1, PETSC_NULL_SF, dmDist, ierr); - CHKERRA(ierr); - if (dmDist /= PETSC_NULL_DM) then - call DMDestroy(dm, ierr); - CHKERRA(ierr); - dm = dmDist; - end if - call DMSetUp(dm,ierr); - CHKERRA(ierr) - - - call DMSetFromOptions(dm, ierr); - CHKERRA(ierr); - - - call TDySetWaterDensityType(tdy,WATER_DENSITY_EXPONENTIAL,ierr); - CHKERRA(ierr) - call DMPlexGetHeightStratum(dm,0,cStart,cEnd,ierr); CHKERRA(ierr); @@ -306,10 +319,8 @@ program main enddo enddo - call TDySetDM(tdy, dm, ierr); - CHKERRA(ierr); - call TDySetFromOptions(tdy,ierr); - CHKERRA(ierr); + call TDySetWaterDensityType(tdy,WATER_DENSITY_EXPONENTIAL,ierr); + CHKERRA(ierr) if (pflotran_consistent) then ! call TDySetPorosityFunction(tdy,PorosityFunctionPFLOTRAN,0,ierr); diff --git a/demo/transient/transient_snes_mpfaof90.cfg b/demo/transient/transient_snes_mpfaof90.cfg index 4339405f..424d79d8 100644 --- a/demo/transient/transient_snes_mpfaof90.cfg +++ b/demo/transient/transient_snes_mpfaof90.cfg @@ -26,7 +26,7 @@ input_arguments=-max_steps 48 -ic_value 90325 -nx 5 -ny 5 -nz 5 -tdy_water_densi input_arguments=-max_steps 10 -use_tdydriver -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename transient-snes-mpfaof90-tdydriver -pc_type lu [transient-snes-mpfaof90-transect-neumann] -input_arguments=-mesh_filename ../../share/meshes/transect_44x1x3_hex_uniform_dz_mesh.exo -ic_value 57180 -dtime 1800 -tdy_mpfao_boundary_condition_type MPFAO_NEUMANN_BC -max_steps 10 -snes_rtol 1e-30 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename transient-snes-mpfaof90-transect-neumann +input_arguments=-tdy_read_mesh ../../share/meshes/transect_44x1x3_hex_uniform_dz_mesh.exo -ic_value 57180 -dtime 1800 -tdy_mpfao_boundary_condition_type MPFAO_NEUMANN_BC -max_steps 10 -snes_rtol 1e-30 -tdy_regression_test -tdy_regression_test_num_cells_per_process 2 -tdy_regression_test_filename transient-snes-mpfaof90-transect-neumann [transient-snes-mpfaof90-region] input_arguments=-max_steps 10 -tdy_regression_test -tdy_regression_test_num_cells_per_process 10 -tdy_regression_test_filename transient-snes-mpfaof90-region -pc_type lu -tdy_connected_region ../../share/meshes/transect_44x1x3_hex_uniform_dz_mesh_regions.bin -mesh_filename ../../share/meshes/transect_44x1x3_hex_uniform_dz_mesh.exo -tdy_mpfao_gmatrix_method MPFAO_GMATRIX_TPF -tdy_mpfao_boundary_condition_type MPFAO_DIRICHLET_BC -tdy_water_density EXPONENTIAL -pflotran_consistent diff --git a/include/finclude/tdycore.h b/include/finclude/tdycore.h index f4b01478..fdc9d1e7 100644 --- a/include/finclude/tdycore.h +++ b/include/finclude/tdycore.h @@ -2,7 +2,7 @@ #define TDYCORECOREDEF_H #define TDy type(tTDy) -#define TDyMethod PetscEnum +#define TDyDiscretization PetscEnum #define MPFAO_DIRICHLET_BC 0 #define MPFAO_NEUMANN_BC 1 diff --git a/include/private/tdybdmimpl.h b/include/private/tdybdmimpl.h new file mode 100644 index 00000000..39e8328c --- /dev/null +++ b/include/private/tdybdmimpl.h @@ -0,0 +1,9 @@ +#if !defined(TDYBDMIMPL_H) +#define TDYBDMIMPL_H + +#include +#include + +PETSC_INTERN PetscErrorCode TDyBDM_Setup(TDy, DM); + +#endif diff --git a/include/private/tdycoreimpl.h b/include/private/tdycoreimpl.h index 89f08d74..1c32167b 100644 --- a/include/private/tdycoreimpl.h +++ b/include/private/tdycoreimpl.h @@ -24,23 +24,37 @@ struct _TDyOps { // Called by TDyDestroy to free implementation-specific resources. PetscErrorCode (*destroy)(void*); + // Creates a DM for a particular simulation (optional). + // Arguments: + // 1. A valid MPI communicator + // 2. A context pointer storing data specifically for the constructor + // function + // 3. A location for the newly created DM. + // We pass the dycore as the first argument here because we don't expect + // the caller to know implementation details. + PetscErrorCode (*create_dm)(MPI_Comm, void*, DM*); + // Implements the view operation for the TDy implementation with the given // viewer. PetscErrorCode (*view)(void*, PetscViewer); - //PetscErrorCode (*setup)(TDy); - // Called by TDySetFromOptions -- sets implementation-specific options // from command-line arguments. - PetscErrorCode (*set_from_options)(void*); + // FIXME: convert the arg here to void* when we've moved specific data + // FIXME: out of TDy. + PetscErrorCode (*set_from_options)(TDy); - // Called by TDySetup -- configures the DM for the dycore. - PetscErrorCode (*config_dm)(void*, DM); + // Called by TDySetup -- configures the DM for solvers. + // FIXME: we should convert the first argument here to void* when we've moved + // FIXME: all discretization-specific data out of TDy itself. + PetscErrorCode (*setup)(TDy, DM); // Called by TDyComputeErrorNorms -- computes error norms given a solution // vector. PetscErrorCode (*compute_error_norms)(void*,Vec,PetscReal*,PetscReal*); + // Functions used to define solver behavior. + // Material and boundary condition functions--we'll sort these out later. PetscErrorCode (*computeporosity)(TDy,PetscReal*,PetscReal*,void*); PetscErrorCode (*computepermeability)(TDy,PetscReal*,PetscReal*,void*); @@ -68,6 +82,9 @@ struct _p_TDy { // configured DM dm; + // Contextual information passed to create_dm (if given). + void* create_dm_context; + // We'll likely get rid of this. TDyTimeIntegrator ti; diff --git a/include/private/tdympfaoimpl.h b/include/private/tdympfaoimpl.h index c87e247a..1f70eacb 100644 --- a/include/private/tdympfaoimpl.h +++ b/include/private/tdympfaoimpl.h @@ -3,7 +3,9 @@ #include -PETSC_INTERN PetscErrorCode TDyMPFAOInitialize(TDy); +PETSC_INTERN PetscErrorCode TDyRichards_MPFA_O_Setup(TDy, DM); +PETSC_INTERN PetscErrorCode TDyRichards_MPFA_O_DAE_Setup(TDy, DM); +PETSC_INTERN PetscErrorCode TDyRichards_MPFA_O_TRANSIENTVAR_Setup(TDy, DM); PETSC_INTERN PetscErrorCode TDyComputeGMatrixMPFAO(TDy); PETSC_INTERN PetscErrorCode TDyComputeGMatrixTPF(TDy); PETSC_INTERN PetscErrorCode TDyUpdateTransmissibilityMatrix(TDy); @@ -21,7 +23,7 @@ PETSC_INTERN PetscErrorCode TDyMPFAOTransientVariable(TS,Vec,Vec,void*); PETSC_INTERN PetscErrorCode TDyMPFAOIFunction_TransientVariable(TS,PetscReal,Vec,Vec,Vec,void*); PETSC_INTERN PetscErrorCode TDyMPFAOSNESFunction(SNES,Vec,Vec,void*); PETSC_INTERN PetscErrorCode TDyMPFAOSNESJacobian(SNES,Vec,Mat,Mat,void*); -PETSC_INTERN PetscErrorCode TDyMPFAOSetFromOptions(TDy); +PETSC_INTERN PetscErrorCode TDyMPFAO_SetFromOptions(TDy); PETSC_INTERN PetscErrorCode TDyMPFAOSNESPreSolve(TDy); #endif diff --git a/include/private/tdyoptions.h b/include/private/tdyoptions.h index f6844a0d..2dbb688c 100644 --- a/include/private/tdyoptions.h +++ b/include/private/tdyoptions.h @@ -13,7 +13,7 @@ typedef struct { PetscInt enthalpy_type; // Numerics settings - TDyMethod method; + TDyDiscretization discretization; TDyQuadratureType qtype; PetscInt mpfao_gmatrix_method; PetscInt mpfao_bc_type; @@ -43,7 +43,6 @@ typedef struct { char init_file[PETSC_MAX_PATH_LEN]; // Mesh-related options - PetscBool generate_mesh; PetscBool read_mesh; char mesh_file[PETSC_MAX_PATH_LEN]; diff --git a/include/private/tdyrichardsimpl.h b/include/private/tdyrichardsimpl.h index 67608fe7..753716f7 100644 --- a/include/private/tdyrichardsimpl.h +++ b/include/private/tdyrichardsimpl.h @@ -4,7 +4,7 @@ #include #include -PETSC_INTERN PetscErrorCode TDyRichardsInitialize(TDy); +PETSC_INTERN PetscErrorCode TDyRichardsInitialize(TDy); // TDyDriver setup fn PETSC_INTERN PetscErrorCode TDyRichardsTSPostStep(TS); PETSC_INTERN PetscErrorCode TDyRichardsSNESPostCheck(SNESLineSearch,Vec,Vec,Vec, PetscBool*,PetscBool*, diff --git a/include/private/tdythimpl.h b/include/private/tdythimpl.h index d0171237..fb8ac240 100644 --- a/include/private/tdythimpl.h +++ b/include/private/tdythimpl.h @@ -4,7 +4,8 @@ #include #include -PETSC_INTERN PetscErrorCode TDyTHInitialize(TDy); +PETSC_INTERN PetscErrorCode TDyTH_MPFA_O_Setup(TDy tdy, DM dm); +PETSC_INTERN PetscErrorCode TDyTHInitialize(TDy); // TDyDriver setup fn PETSC_INTERN PetscErrorCode TDyTHTSPostStep(TS); PETSC_INTERN PetscErrorCode TDyTHSNESPostCheck(SNESLineSearch,Vec,Vec,Vec, PetscBool*,PetscBool*, diff --git a/include/private/tdywyimpl.h b/include/private/tdywyimpl.h new file mode 100644 index 00000000..b24f2c52 --- /dev/null +++ b/include/private/tdywyimpl.h @@ -0,0 +1,9 @@ +#if !defined(TDYWYIMPL_H) +#define TDYWYIMPL_H + +#include +#include + +PETSC_INTERN PetscErrorCode TDyWY_Setup(TDy, DM); + +#endif diff --git a/include/tdycore.h b/include/tdycore.h index 6d909c0b..39dc62b9 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -6,16 +6,34 @@ /* ---------------------------------------------------------------- */ +/// This type lists modes that identify the set of governing equations solved +/// by the dycore. typedef enum { - TPF=0, /* two point flux, classic finite volumes */ - MPFA_O, /* multipoint flux approximation - O method */ - MPFA_O_DAE, /* multipoint flux approximation - O method using DAE */ - MPFA_O_TRANSIENTVAR, /* multipoint flux approximation - O method using TS transient (conservative) approach */ - BDM, /* P0,BDM1 spaces, standard approach */ - WY /* P0,BDM1 spaces, vertex quadrature, statically condensed */ -} TDyMethod; + /// Richards equation + RICHARDS=0, + /// Non-isothermal flows + TH +} TDyMode; + +PETSC_EXTERN const char *const TDyModes[]; -PETSC_EXTERN const char *const TDyMethods[]; +/// This type enumerates discretizations supported by the dycore. +typedef enum { + /// multi-point flux approximation - O method + MPFA_O=0, + /// multi-point flux approximation - O method using DAE + MPFA_O_DAE, + /// multipoint flux approximation - O method using TS transient (conservative) + /// approach + MPFA_O_TRANSIENTVAR, + /// finite element using P0, BDM1 spaces, standard approach + BDM, + /// finite element using P0,BDM1 spaces, vertex quadrature, statically + /// condensed + WY +} TDyDiscretization; + +PETSC_EXTERN const char *const TDyDiscretizations[]; typedef enum { MPFAO_GMATRIX_DEFAULT=0, /* default method to compute gmatrix for MPFA-O method */ @@ -39,11 +57,6 @@ typedef enum { PETSC_EXTERN const char *const TDyQuadratureTypes[]; -typedef enum { - RICHARDS=0, - TH -} TDyMode; - typedef enum { TDySNES=0, TDyTS @@ -52,12 +65,12 @@ typedef enum { typedef enum { TDyCreated=0x0, TDyParametersInitialized=0x1, - TDyOptionsSet=0x2, - TDySetupFinished=0x4, + TDyModeSet=0x1<<1, + TDyDiscretizationSet=0x1<<2, + TDyOptionsSet=0x1<<3, + TDySetupFinished=0x1<<4, } TDySetupFlags; -PETSC_EXTERN const char *const TDyModes[]; - typedef void (*SpatialFunction)(PetscReal *x,PetscReal *f); /* returns f(x) */ typedef struct _p_TDy *TDy; @@ -79,13 +92,14 @@ PETSC_EXTERN PetscErrorCode TDyInitNoArguments(void); PETSC_EXTERN PetscErrorCode TDyOnFinalize(void (*)(void)); PETSC_EXTERN PetscErrorCode TDyFinalize(void); -PETSC_EXTERN PetscErrorCode TDyCreate(TDy*); +PETSC_EXTERN PetscErrorCode TDyCreate(MPI_Comm, TDy*); PETSC_EXTERN PetscErrorCode TDySetMode(TDy,TDyMode); -PETSC_EXTERN PetscErrorCode TDySetDiscretizationMethod(TDy,TDyMethod); +PETSC_EXTERN PetscErrorCode TDySetDiscretization(TDy,TDyDiscretization); +PETSC_EXTERN PetscErrorCode TDySetDMConstructor(TDy,void*,PetscErrorCode(*)(MPI_Comm, void*, DM*)); PETSC_EXTERN PetscErrorCode TDySetFromOptions(TDy); PETSC_EXTERN PetscErrorCode TDySetup(TDy); -PETSC_EXTERN PetscErrorCode TDyDestroy(TDy *tdy); -PETSC_EXTERN PetscErrorCode TDyView(TDy,PetscViewer viewer); +PETSC_EXTERN PetscErrorCode TDyDestroy(TDy*); +PETSC_EXTERN PetscErrorCode TDyView(TDy,PetscViewer); PETSC_EXTERN PetscErrorCode TDyGetDimension(TDy,PetscInt*); PETSC_EXTERN PetscErrorCode TDyGetDM(TDy,DM*); @@ -146,8 +160,7 @@ PETSC_EXTERN PetscErrorCode TDyGetNumCellsLocal(TDy,PetscInt*); PETSC_EXTERN PetscErrorCode TDyGetCellNaturalIDsLocal(TDy,PetscInt*,PetscInt[]); PETSC_EXTERN PetscErrorCode TDyGetCellIsLocal(TDy,PetscInt*,PetscInt[]); - -PETSC_EXTERN PetscErrorCode TDyResetDiscretizationMethod(TDy); +PETSC_EXTERN PetscErrorCode TDyResetDiscretization(TDy); PETSC_EXTERN PetscErrorCode TDySetQuadratureType(TDy,TDyQuadratureType); PETSC_EXTERN PetscErrorCode TDySetWaterDensityType(TDy,TDyWaterDensityType); @@ -155,7 +168,6 @@ PETSC_EXTERN PetscErrorCode TDySetMPFAOGmatrixMethod(TDy,TDyMPFAOGmatrixMethod); PETSC_EXTERN PetscErrorCode TDySetMPFAOBoundaryConditionType(TDy,TDyMPFAOBoundaryConditionType); // We will probably remove the following functions. -PETSC_EXTERN PetscErrorCode TDySetDM(TDy,DM); PETSC_EXTERN PetscErrorCode TDyComputeSystem(TDy,Mat,Vec); PETSC_EXTERN PetscErrorCode TDySetIFunction(TS,TDy); PETSC_EXTERN PetscErrorCode TDySetIJacobian(TS,TDy); @@ -172,16 +184,8 @@ PETSC_EXTERN PetscErrorCode TDyPostSolveSNESSolver(TDy,Vec); PETSC_EXTERN PetscErrorCode TDyOutputRegression(TDy,Vec); -PETSC_EXTERN PetscErrorCode TDyTPFInitialize(TDy); -PETSC_EXTERN PetscErrorCode TDyTPFComputeSystem(TDy,Mat,Vec); -PETSC_EXTERN PetscReal TDyTPFPressureNorm(TDy,Vec); -PETSC_EXTERN PetscReal TDyTPFVelocityNorm(TDy,Vec); -PETSC_EXTERN PetscErrorCode TDyTPFCheckMeshSuitability(TDy); - -PETSC_EXTERN PetscErrorCode TDyWYInitialize(TDy); PETSC_EXTERN PetscErrorCode TDyWYComputeSystem(TDy,Mat,Vec); -PETSC_EXTERN PetscErrorCode TDyBDMInitialize(TDy); PETSC_EXTERN PetscErrorCode TDyBDMComputeSystem(TDy,Mat,Vec); PETSC_EXTERN PetscReal TDyBDMPressureNorm(TDy,Vec); PETSC_EXTERN PetscReal TDyBDMVelocityNorm(TDy,Vec); diff --git a/src/f90-mod/tdycore.h b/src/f90-mod/tdycore.h index 46f62e7c..9fc37a81 100644 --- a/src/f90-mod/tdycore.h +++ b/src/f90-mod/tdycore.h @@ -10,7 +10,6 @@ ! ! Types of TDy methods ! - PetscEnum TPF PetscEnum MPFA_O PetscEnum MPFA_O_DAE PetscEnum MPFA_O_TRANSIENT_VAR @@ -19,8 +18,8 @@ ! The parameters values need to match those defined in ! the C code (i.e. include/tdycore.h) - parameter (TPF=0,MPFA_O=1,MPFA_O_DAE=2,MPFA_O_TRANSIENT_VAR=3) - parameter (BDM=4,WY=5) + parameter (MPFA_O=0,MPFA_O_DAE=1,MPFA_O_TRANSIENT_VAR=2) + parameter (BDM=3,WY=4) ! ! Types of TDy modes diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index f217c9da..6a9df75b 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -31,21 +31,12 @@ subroutine TDySetMode(a,b,z) end subroutine TDySetMode end interface interface - subroutine TDySetDiscretizationMethod(a,b,z) + subroutine TDySetDiscretization(a,b,z) use tdycoredef TDy a PetscInt b integer z - end subroutine TDySetDiscretizationMethod - end interface - interface - subroutine TDySetDM(a,b,z) - use petscdm - use tdycoredef - TDy a - DM b - integer z - end subroutine TDySetDM + end subroutine TDySetDiscretization end interface interface subroutine TDySetFromOptions(a,z) @@ -461,6 +452,17 @@ function GetRegFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) end function end interface + abstract interface + subroutine TDyDMConstructor(comm, dm, ierr) + use, intrinsic :: iso_c_binding + use tdycoredef + use petscdm + MPI_Comm :: comm + DM :: dm + PetscErrorCode :: ierr + end subroutine + end interface + contains subroutine TDyInit(ierr) @@ -494,6 +496,27 @@ function RegisterFn(name, func) bind (c, name="TDyRegisterFunction") result(ierr ierr = RegisterFn(FtoCString(name), c_funloc(func)) end subroutine + subroutine TDySetDMConstructor(tdy, dm_ctor, ierr) + use, intrinsic :: iso_c_binding + use tdycoredef + implicit none + TDy, target :: tdy + procedure(TDyDMConstructor) :: dm_ctor + PetscErrorCode :: ierr + + interface + function SetDMConstructor(tdy, dm_ctor) bind (c, name="TDySetDMConstructorF90") result(ierr) + use, intrinsic :: iso_c_binding + implicit none + type(c_ptr), value :: tdy + type(c_funptr), value :: dm_ctor + integer(c_int) :: ierr + end function + end interface + + ierr = SetDMConstructor(c_loc(tdy), c_funloc(dm_ctor)) + end subroutine + subroutine TDySelectPorosityFunction(tdy, name, ierr) use, intrinsic :: iso_c_binding use tdycoredef diff --git a/src/interface/ftn/tdycoref.c b/src/interface/ftn/tdycoref.c index 4058835b..b250b3a1 100644 --- a/src/interface/ftn/tdycoref.c +++ b/src/interface/ftn/tdycoref.c @@ -10,15 +10,15 @@ #define tdyinitnoarguments_ TDYINITNOARGUMENTS #define tdyfinalize_ TDYFINALIZE #define tdycreate_ TDYCREATE -#define tdycreatesetdm_ TDYSETDM #define tdysetmode_ TDYSETMODE -#define tdysetdiscretizationmethod_ TDYSETDISCRETIZATIONMETHOD +#define tdysetdiscretization_ TDYSETDISCRETIZATION #define tdysetfromoptions_ TDYSETFROMOPTIONS #define tdydriverinitializettdy_ TDYDRIVERINITIALIZETDY #define tdydtimeintegratorruntotime_ TDYTIMEINTEGRATORRUNTOTIME #define tdydtimeintegratorsettimestep_ TDYTIMEINTEGRATORSETTIMESTEP #define tdydtimeintegratoroutputregression_ TDYTIMEINTEGRATOROUTPUTREGRESSION #define tdysetup_ TDYSETUP +#define tdygetdm_ TDYGETDM #define tdysetwaterdensitytype_ TDYSETWATERDENSITYTYPE #define tdysetmpfaogmatrixmethod_ TDYSETMPFAOGMATRIXMETHOD #define tdysetmpfaoboundaryconditiontype_ TDYSETMPFAOGBOUNDARYCONDITIONTYPE @@ -72,15 +72,15 @@ #define tdyinitnoarguments_ tdyinitnoarguments #define tdyfinalize_ tdyfinalize #define tdycreate_ tdycreate -#define tdysetdm_ tdysetdm #define tdysetmode_ tdysetmode -#define tdysetdiscretizationmethod_ tdysetdiscretizationmethod +#define tdysetdiscretization_ tdysetdiscretization #define tdysetfromoptions_ tdysetfromoptions #define tdydriverinitializettdy_ tdydriverinitializetdy #define tdydtimeintegratorruntotime_ tdydtimeintegratorruntotime #define tdydtimeintegratorsettimestep_ tdydtimeintegratorsettimestep #define tdydtimeintegratoroutputregression_ tdydtimeintegratoroutputregression #define tdysetup_ tdysetup +#define tdygetdm_ tdygetdm #define tdysetwaterdensitytype_ tdysetwaterdensitytype #define tdysetmpfaogmatrixmethod_ tdysetmpfaogmatrixmethod #define tdysetmpfaoboundaryconditiontype_ tdysetmpfaoboundaryconditiontype @@ -167,7 +167,7 @@ PETSC_EXTERN void tdyfinalize_(int *__ierr){ extern "C" { #endif PETSC_EXTERN void tdycreate_(TDy *_tdy, int *__ierr){ -*__ierr = TDyCreate(_tdy); +*__ierr = TDyCreate(PETSC_COMM_WORLD,_tdy); } #if defined(__cplusplus) } @@ -185,18 +185,8 @@ PETSC_EXTERN void tdysetmode_(TDy _tdy, PetscInt *mode, int *__ierr){ #if defined(__cplusplus) extern "C" { #endif -PETSC_EXTERN void tdysetdiscretizationmethod_(TDy _tdy, PetscInt *method, int *__ierr){ -*__ierr = TDySetDiscretizationMethod((TDy)PetscToPointer((_tdy)), *method); -} -#if defined(__cplusplus) -} -#endif - -#if defined(__cplusplus) -extern "C" { -#endif -PETSC_EXTERN void tdysetdm_(TDy _tdy, DM dm, int *__ierr){ -*__ierr = TDySetDM((TDy)PetscToPointer((_tdy)), (DM)PetscToPointer((dm))); +PETSC_EXTERN void tdysetdiscretization_(TDy _tdy, PetscInt *discretization, int *__ierr){ +*__ierr = TDySetDiscretization((TDy)PetscToPointer((_tdy)), *discretization); } #if defined(__cplusplus) } @@ -262,6 +252,16 @@ PETSC_EXTERN void tdysetup_(TDy _tdy, int *__ierr){ } #endif +#if defined(__cplusplus) +extern "C" { +#endif +PETSC_EXTERN void tdygetdm_(TDy _tdy, DM *dm, int *__ierr){ +*__ierr = TDyGetDM((TDy)PetscToPointer((_tdy)), dm); +} +#if defined(__cplusplus) +} +#endif + #if defined(__cplusplus) extern "C" { #endif diff --git a/src/mpfao/tdympfao.c b/src/mpfao/tdympfao.c index 1169e75d..6d81c2bb 100644 --- a/src/mpfao/tdympfao.c +++ b/src/mpfao/tdympfao.c @@ -307,15 +307,206 @@ PetscErrorCode TDyMPFAO_AllocateMemoryForEnergySourceSinkValues(TDy tdy) { TDY_STOP_FUNCTION_TIMER() PetscFunctionReturn(0); } - /* -------------------------------------------------------------------------- */ -PetscErrorCode TDyMPFAOInitialize(TDy tdy) { +//----------------- +// Setup functions +//----------------- + +// There's a lot of duplicated code here for the moment. We'll factor out the +// repeated stuff when we move the data out of the TDy struct and into +// discretization-specific types. +// Setup function for Richards + MPFA_O +PetscErrorCode TDyRichards_MPFA_O_Setup(TDy tdy, DM dm) { + PetscErrorCode ierr; + PetscFunctionBegin; + + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + + PetscInt dim; + ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); + if (dim == 2) { + SETERRQ(comm,PETSC_ERR_USER,"MPFA-O method supports only 3D calculations."); + } + + // Allocate mesh data. + tdy->mesh = (TDyMesh *) malloc(sizeof(TDyMesh)); + ierr = TDyAllocateMemoryForMesh(tdy); CHKERRQ(ierr); + ierr = TDyAllocate_RealArray_3D(&tdy->Trans, tdy->mesh->num_vertices, tdy->nfv, tdy->nfv + tdy->ncv); CHKERRQ(ierr); + ierr = PetscMalloc(tdy->mesh->num_faces*sizeof(PetscReal), + &(tdy->vel )); CHKERRQ(ierr); + ierr = TDyInitialize_RealArray_1D(tdy->vel, tdy->mesh->num_faces, 0.0); CHKERRQ(ierr); + ierr = PetscMalloc(tdy->mesh->num_faces*sizeof(PetscInt), + &(tdy->vel_count)); CHKERRQ(ierr); + ierr = TDyInitialize_IntegerArray_1D(tdy->vel_count, tdy->mesh->num_faces, 0); CHKERRQ(ierr); + PetscInt nsubcells = 8; + PetscInt nrow = 3; + PetscInt ncol = 3; + ierr = TDyAllocate_RealArray_4D(&tdy->subc_Gmatrix, tdy->mesh->num_cells, + nsubcells, nrow, ncol); CHKERRQ(ierr); + + // Set up the section, 1 dof per cell + PetscSection sec; + PetscInt p, pStart, pEnd; + ierr = PetscSectionCreate(comm, &sec); CHKERRQ(ierr); + ierr = PetscSectionSetNumFields(sec, 1); CHKERRQ(ierr); + ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); + ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); + + ierr = DMPlexGetHeightStratum(dm,0,&pStart,&pEnd); CHKERRQ(ierr); + ierr = PetscSectionSetChart(sec,pStart,pEnd); CHKERRQ(ierr); + for(p=pStart; pmesh->num_faces; + nLocalCells = TDyMeshGetNumberOfLocalCells(tdy->mesh); + nNonLocalFaces = TDyMeshGetNumberOfNonLocalFacess(tdy->mesh); + nNonInternalFaces = TDyMeshGetNumberOfNonInternalFacess(tdy->mesh); + + nrow = 4*nFaces; + ncol = nLocalCells + nNonLocalFaces + nNonInternalFaces; + nz = tdy->nfv; + ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,nrow,ncol,nz,NULL,&tdy->Trans_mat); CHKERRQ(ierr); + ierr = MatSetOption(tdy->Trans_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + ierr = VecCreateSeq(PETSC_COMM_SELF,ncol,&tdy->P_vec); + ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->TtimesP_vec); + ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->GravDisVec); + ierr = VecZeroEntries(tdy->GravDisVec); + } + + // Set up material models. + if (tdy->ops->computeporosity) { ierr = SetPorosityFromFunction(tdy); CHKERRQ(ierr); } + if (tdy->ops->computepermeability) {ierr = SetPermeabilityFromFunction(tdy); CHKERRQ(ierr);} + + // why must these be placed after SetPermeabilityFromFunction()? + ierr = TDyComputeGMatrix(tdy); CHKERRQ(ierr); + ierr = TDyComputeTransmissibilityMatrix(tdy); CHKERRQ(ierr); + ierr = TDyComputeGravityDiscretization(tdy); CHKERRQ(ierr); + + ierr = TDyMPFAO_AllocateMemoryForBoundaryValues(tdy); CHKERRQ(ierr); + ierr = TDyMPFAO_AllocateMemoryForSourceSinkValues(tdy); CHKERRQ(ierr); + + PetscFunctionReturn(0); +} + +// Setup function for Richards + MPFA_O_DAE +PetscErrorCode TDyRichards_MPFA_O_DAE_Setup(TDy tdy, DM dm) { + PetscErrorCode ierr; + PetscFunctionBegin; + PetscFunctionReturn(0); + + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + + PetscInt dim; + ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); + if (dim == 2) { + SETERRQ(comm,PETSC_ERR_USER,"MPFA-O method supports only 3D calculations."); + } + + // Allocate mesh data. + tdy->mesh = (TDyMesh *) malloc(sizeof(TDyMesh)); + ierr = TDyAllocateMemoryForMesh(tdy); CHKERRQ(ierr); + ierr = TDyAllocate_RealArray_3D(&tdy->Trans, tdy->mesh->num_vertices, tdy->nfv, tdy->nfv + tdy->ncv); CHKERRQ(ierr); + ierr = PetscMalloc(tdy->mesh->num_faces*sizeof(PetscReal), + &(tdy->vel )); CHKERRQ(ierr); + ierr = TDyInitialize_RealArray_1D(tdy->vel, tdy->mesh->num_faces, 0.0); CHKERRQ(ierr); + ierr = PetscMalloc(tdy->mesh->num_faces*sizeof(PetscInt), + &(tdy->vel_count)); CHKERRQ(ierr); + ierr = TDyInitialize_IntegerArray_1D(tdy->vel_count, tdy->mesh->num_faces, 0); CHKERRQ(ierr); + PetscInt nsubcells = 8; + PetscInt nrow = 3; + PetscInt ncol = 3; + ierr = TDyAllocate_RealArray_4D(&tdy->subc_Gmatrix, tdy->mesh->num_cells, + nsubcells, nrow, ncol); CHKERRQ(ierr); + + // Set up the section, 1 dof per cell + PetscSection sec; + PetscInt p, pStart, pEnd; + ierr = PetscSectionCreate(comm, &sec); CHKERRQ(ierr); + ierr = PetscSectionSetNumFields(sec, 2); CHKERRQ(ierr); + ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); + ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); + ierr = PetscSectionSetFieldName(sec, 1, "LiquidMass"); CHKERRQ(ierr); + ierr = PetscSectionSetFieldComponents(sec, 1, 1); CHKERRQ(ierr); + + ierr = DMPlexGetHeightStratum(dm,0,&pStart,&pEnd); CHKERRQ(ierr); + ierr = PetscSectionSetChart(sec,pStart,pEnd); CHKERRQ(ierr); + for(p=pStart; pmesh->num_faces; + nLocalCells = TDyMeshGetNumberOfLocalCells(tdy->mesh); + nNonLocalFaces = TDyMeshGetNumberOfNonLocalFacess(tdy->mesh); + nNonInternalFaces = TDyMeshGetNumberOfNonInternalFacess(tdy->mesh); + + nrow = 4*nFaces; + ncol = nLocalCells + nNonLocalFaces + nNonInternalFaces; + nz = tdy->nfv; + ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,nrow,ncol,nz,NULL,&tdy->Trans_mat); CHKERRQ(ierr); + ierr = MatSetOption(tdy->Trans_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + ierr = VecCreateSeq(PETSC_COMM_SELF,ncol,&tdy->P_vec); + ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->TtimesP_vec); + ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->GravDisVec); + ierr = VecZeroEntries(tdy->GravDisVec); + } + + // Set up material models. + if (tdy->ops->computeporosity) { ierr = SetPorosityFromFunction(tdy); CHKERRQ(ierr); } + if (tdy->ops->computepermeability) {ierr = SetPermeabilityFromFunction(tdy); CHKERRQ(ierr);} + + // why must these be placed after SetPermeabilityFromFunction()? + ierr = TDyComputeGMatrix(tdy); CHKERRQ(ierr); + ierr = TDyComputeTransmissibilityMatrix(tdy); CHKERRQ(ierr); + ierr = TDyComputeGravityDiscretization(tdy); CHKERRQ(ierr); + + ierr = TDyMPFAO_AllocateMemoryForBoundaryValues(tdy); CHKERRQ(ierr); + ierr = TDyMPFAO_AllocateMemoryForSourceSinkValues(tdy); CHKERRQ(ierr); + +} + +// Setup function for Richards + MPFA_O_TRANSIENTVAR +PetscErrorCode TDyRichards_MPFA_O_TRANSIENTVAR_Setup(TDy tdy, DM dm) { + PetscErrorCode ierr; + PetscFunctionBegin; + // This is essentially the same as the MPFA-O one. + ierr = TDyRichards_MPFA_O_Setup(tdy, dm); + PetscFunctionReturn(0); +} + +// Setup function for TH + MPFA-O +PetscErrorCode TDyTH_MPFA_O_Setup(TDy tdy, DM dm) { PetscErrorCode ierr; - MPI_Comm comm; - DM dm = tdy->dm; - PetscInt nrow,ncol,nsubcells; PetscFunctionBegin; TDY_START_FUNCTION_TIMER() @@ -326,16 +517,15 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"MPFA-O method supports only 3D calculations."); } + MPI_Comm comm; ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); + // Allocate mesh data. tdy->mesh = (TDyMesh *) malloc(sizeof(TDyMesh)); - ierr = TDyAllocateMemoryForMesh(tdy); CHKERRQ(ierr); ierr = TDyAllocate_RealArray_3D(&tdy->Trans, tdy->mesh->num_vertices, tdy->nfv, tdy->nfv + tdy->ncv); CHKERRQ(ierr); - if (tdy->options.mode == TH) { - ierr = TDyAllocate_RealArray_3D(&tdy->Temp_Trans, tdy->mesh->num_vertices, - tdy->nfv, tdy->nfv); CHKERRQ(ierr); - } + ierr = TDyAllocate_RealArray_3D(&tdy->Temp_Trans, tdy->mesh->num_vertices, + tdy->nfv, tdy->nfv); CHKERRQ(ierr); ierr = PetscMalloc(tdy->mesh->num_faces*sizeof(PetscReal), &(tdy->vel )); CHKERRQ(ierr); ierr = TDyInitialize_RealArray_1D(tdy->vel, tdy->mesh->num_faces, 0.0); CHKERRQ(ierr); @@ -343,74 +533,41 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { &(tdy->vel_count)); CHKERRQ(ierr); ierr = TDyInitialize_IntegerArray_1D(tdy->vel_count, tdy->mesh->num_faces, 0); CHKERRQ(ierr); - nsubcells = 8; - nrow = 3; - ncol = 3; + PetscInt nsubcells = 8; + PetscInt nrow = 3; + PetscInt ncol = 3; ierr = TDyAllocate_RealArray_4D(&tdy->subc_Gmatrix, tdy->mesh->num_cells, nsubcells, nrow, ncol); CHKERRQ(ierr); - if (tdy->options.mode == TH) { - ierr = TDyAllocate_RealArray_4D(&tdy->Temp_subc_Gmatrix, tdy->mesh->num_cells, - nsubcells, nrow, ncol); CHKERRQ(ierr); - } + ierr = TDyAllocate_RealArray_4D(&tdy->Temp_subc_Gmatrix, tdy->mesh->num_cells, + nsubcells, nrow, ncol); CHKERRQ(ierr); - /* Set up the section, 1 dof per cell */ + // Set up the section, 1 dof per cell PetscSection sec; PetscInt p, pStart, pEnd; - PetscBool use_dae; - - use_dae = (tdy->options.method == MPFA_O_DAE); ierr = PetscSectionCreate(comm, &sec); CHKERRQ(ierr); - if (!use_dae) { - if (tdy->options.mode == TH){ - ierr = PetscSectionSetNumFields(sec, 2); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 1, "LiquidTemperature"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 1, 1); CHKERRQ(ierr); - } else { - ierr = PetscSectionSetNumFields(sec, 1); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); - } - } else { - if (tdy->options.mode == TH) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TH unsupported with MPFA_O_DAE"); - } - ierr = PetscSectionSetNumFields(sec, 2); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec, 1, "LiquidMass"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec, 1, 1); CHKERRQ(ierr); - } + ierr = PetscSectionSetNumFields(sec, 2); CHKERRQ(ierr); + ierr = PetscSectionSetFieldName(sec, 0, "LiquidPressure"); CHKERRQ(ierr); + ierr = PetscSectionSetFieldComponents(sec, 0, 1); CHKERRQ(ierr); + ierr = PetscSectionSetFieldName(sec, 1, "LiquidTemperature"); CHKERRQ(ierr); + ierr = PetscSectionSetFieldComponents(sec, 1, 1); CHKERRQ(ierr); ierr = DMPlexGetHeightStratum(dm,0,&pStart,&pEnd); CHKERRQ(ierr); ierr = PetscSectionSetChart(sec,pStart,pEnd); CHKERRQ(ierr); for(p=pStart; poptions.mode == TH) { - ierr = PetscSectionSetFieldDof(sec,p,1,1); CHKERRQ(ierr); - ierr = PetscSectionSetDof(sec,p,2); CHKERRQ(ierr); - } - if (use_dae) { - if (tdy->options.mode == TH) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TH unsupported with MPFA_O_DAE"); - } - ierr = PetscSectionSetFieldDof(sec,p,1,1); CHKERRQ(ierr); - ierr = PetscSectionSetDof(sec,p,2); CHKERRQ(ierr); - } + ierr = PetscSectionSetFieldDof(sec,p,1,1); CHKERRQ(ierr); + ierr = PetscSectionSetDof(sec,p,2); CHKERRQ(ierr); } ierr = PetscSectionSetUp(sec); CHKERRQ(ierr); ierr = DMSetSection(dm,sec); CHKERRQ(ierr); ierr = PetscSectionViewFromOptions(sec, NULL, "-layout_view"); CHKERRQ(ierr); ierr = PetscSectionDestroy(&sec); CHKERRQ(ierr); ierr = DMSetBasicAdjacency(dm,PETSC_TRUE,PETSC_TRUE); CHKERRQ(ierr); - //ierr = DMPlexSetAdjacencyUseCone(dm,PETSC_TRUE);CHKERRQ(ierr); - //ierr = DMPlexSetAdjacencyUseClosure(dm,PETSC_TRUE);CHKERRQ(ierr); + // Build the mesh. ierr = TDyBuildMesh(tdy); CHKERRQ(ierr); - { PetscInt nLocalCells, nFaces, nNonLocalFaces, nNonInternalFaces; PetscInt nrow, ncol, nz; @@ -430,25 +587,22 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->GravDisVec); ierr = VecZeroEntries(tdy->GravDisVec); - if (tdy->options.mode == TH) { - ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,nrow,ncol,nz,NULL,&tdy->Temp_Trans_mat); CHKERRQ(ierr); - ierr = MatSetOption(tdy->Temp_Trans_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); - ierr = VecCreateSeq(PETSC_COMM_SELF,ncol,&tdy->Temp_P_vec); - ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->Temp_TtimesP_vec); - } + ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,nrow,ncol,nz,NULL,&tdy->Temp_Trans_mat); CHKERRQ(ierr); + ierr = MatSetOption(tdy->Temp_Trans_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + ierr = VecCreateSeq(PETSC_COMM_SELF,ncol,&tdy->Temp_P_vec); + ierr = VecCreateSeq(PETSC_COMM_SELF,nrow,&tdy->Temp_TtimesP_vec); } + // Set up material models. if (tdy->ops->computeporosity) { ierr = SetPorosityFromFunction(tdy); CHKERRQ(ierr); } if (tdy->ops->computepermeability) {ierr = SetPermeabilityFromFunction(tdy); CHKERRQ(ierr);} - if (tdy->options.mode == TH){ - if (tdy->ops->computethermalconductivity) {ierr = SetThermalConductivityFromFunction(tdy); CHKERRQ(ierr);} - if (tdy->ops->computesoildensity) { - ierr = SetSoilDensityFromFunction(tdy); CHKERRQ(ierr); - } - if (tdy->ops->computesoilspecificheat) { - printf("SetSoilSpecificHeatFromFunction\n"); - ierr = SetSoilSpecificHeatFromFunction(tdy); CHKERRQ(ierr); - } + if (tdy->ops->computethermalconductivity) {ierr = SetThermalConductivityFromFunction(tdy); CHKERRQ(ierr);} + if (tdy->ops->computesoildensity) { + ierr = SetSoilDensityFromFunction(tdy); CHKERRQ(ierr); + } + if (tdy->ops->computesoilspecificheat) { + printf("SetSoilSpecificHeatFromFunction\n"); + ierr = SetSoilSpecificHeatFromFunction(tdy); CHKERRQ(ierr); } // why must these be placed after SetPermeabilityFromFunction()? @@ -469,7 +623,7 @@ PetscErrorCode TDyMPFAOInitialize(TDy tdy) { } /* -------------------------------------------------------------------------- */ -PetscErrorCode TDyMPFAOSetFromOptions(TDy tdy) { +PetscErrorCode TDyMPFAO_SetFromOptions(TDy tdy) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() diff --git a/src/tdybdm.c b/src/tdybdm.c index 7b865a7f..156880f4 100644 --- a/src/tdybdm.c +++ b/src/tdybdm.c @@ -1,4 +1,5 @@ #include +#include #include #include @@ -133,7 +134,7 @@ void HdivBasisHex(const PetscReal *x,PetscReal *B,PetscReal *DF,PetscReal J) { TDY_STOP_FUNCTION_TIMER() } -PetscErrorCode TDyBDMInitialize(TDy tdy) { +PetscErrorCode TDyBDM_Setup(TDy tdy, DM dm) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() PetscErrorCode ierr; @@ -142,7 +143,6 @@ PetscErrorCode TDyBDMInitialize(TDy tdy) { PetscSection sec; PetscInt d,dim,dofs_per_face = 1; PetscBool found; - DM dm = tdy->dm; ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); /* Get plex limits */ diff --git a/src/tdycore.c b/src/tdycore.c index a0a5a1b3..9aa1acd8 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -10,6 +10,9 @@ static char help[] = "TDycore \n\ #include #include #include +#include +#include +#include #include #include #include @@ -20,15 +23,14 @@ static char help[] = "TDycore \n\ #include #include -const char *const TDyMethods[] = { - "TPF", +const char *const TDyDiscretizations[] = { "MPFA_O", "MPFA_O_DAE", "MPFA_O_TRANSIENTVAR", "BDM", "WY", /* */ - "TDyMethod","TDY_METHOD_",NULL + "TDyDiscretization","TDY_DISCRETIZATION_",NULL }; const char *const TDyMPFAOGmatrixMethods[] = { @@ -195,7 +197,7 @@ static PetscErrorCode SetDefaultOptions(TDy tdy) { TDyOptions *options = &tdy->options; options->mode = RICHARDS; - options->method = WY; + options->discretization = MPFA_O; options->gravity_constant = 9.8068; options->rho_type = WATER_DENSITY_CONSTANT; options->mu_type = WATER_VISCOSITY_CONSTANT; @@ -228,7 +230,7 @@ static PetscErrorCode SetDefaultOptions(TDy tdy) { PetscFunctionReturn(0); } -PetscErrorCode TDyCreate(TDy *_tdy) { +PetscErrorCode TDyCreate(MPI_Comm comm, TDy *_tdy) { TDy tdy; PetscErrorCode ierr; PetscFunctionBegin; @@ -238,7 +240,7 @@ PetscErrorCode TDyCreate(TDy *_tdy) { // Initialize TDycore-specific subsystems. ierr = TDyInitSubsystems(); CHKERRQ(ierr); - ierr = PetscHeaderCreate(tdy,TDY_CLASSID,"TDy","TDy","TDy",PETSC_COMM_WORLD, + ierr = PetscHeaderCreate(tdy,TDY_CLASSID,"TDy","TDy","TDy",comm, TDyDestroy,TDyView); CHKERRQ(ierr); *_tdy = tdy; tdy->setup_flags |= TDyCreated; @@ -264,21 +266,42 @@ PetscErrorCode TDyCreate(TDy *_tdy) { PetscFunctionReturn(0); } -PetscErrorCode TDySetDM(TDy tdy, DM dm) { +/// Provides a function to be used to create a DM in special cases where a +/// specific geometry is needed. After the function is executed on the DM, +/// DMSetFromOptions is called to apply overrides, and then the DM is +/// distributed appropriately. This function must be called before +/// TDySetFromOptions. +/// @param [inout] tdy the dycore +/// @param [in] context a pointer to contextual information that can be used by +/// dm_func to create a DM. This pointer is not managed +/// by the dycore. +/// @param [in] dm_func A function that, given an MPI communicator and a context +/// pointer, creates a given DM and returns an error +PetscErrorCode TDySetDMConstructor(TDy tdy, void* context, + PetscErrorCode (*dm_func)(MPI_Comm, void*, DM*)) { PetscFunctionBegin; - if (tdy->dm != NULL) { - if (tdy->options.generate_mesh) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, - "Can't call TDySetDM: A DM has already been generated from supplied options."); - } else { // we read a mesh from a file - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, - "Can't call TDySetDM: A DM has already been read from a given file."); - } + MPI_Comm comm; + int ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + if ((tdy->setup_flags & TDyDiscretizationSet) == 0) { + SETERRQ(comm,PETSC_ERR_USER, + "You must call TDySetDiscretization before TDySetDMConstructor()"); } - if (!dm) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"A DM must be created prior to TDySetDM()"); + if (tdy->setup_flags & TDyOptionsSet) { + SETERRQ(comm,PETSC_ERR_USER, + "You must call TDySetDMConstructor before TDySetFromOptions()"); } - tdy->dm = dm; + tdy->create_dm_context = context; + tdy->ops->create_dm = dm_func; + PetscFunctionReturn(0); +} + +// This function is a wrapper used to eliminate the context pointer argument +// from TDySetDMConstructor so the function can be called from Fortran. +static void (*create_dm_f90_)(MPI_Fint*, DM*, PetscErrorCode*) = NULL; +PetscErrorCode TDySetDMConstructorF90(TDy tdy, + void (*dm_func)(MPI_Fint*, DM*, PetscErrorCode*)) { + PetscFunctionBegin; + create_dm_f90_ = dm_func; PetscFunctionReturn(0); } @@ -286,8 +309,11 @@ PetscErrorCode TDyMalloc(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + if (!tdy->dm) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"tdy->dm must be set prior to TDyMalloc()"); + SETERRQ(comm,PETSC_ERR_USER,"tdy->dm must be set prior to TDyMalloc()"); } PetscInt dim; @@ -366,23 +392,35 @@ PetscErrorCode TDyCreateGrid(TDy tdy) { PetscScalar *coords; PetscErrorCode ierr; PetscFunctionBegin; - MPI_Comm comm = PETSC_COMM_WORLD; + + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + if ((tdy->setup_flags & TDyOptionsSet) == 0) { SETERRQ(comm,PETSC_ERR_USER,"Options must be set prior to TDyCreateGrid()"); } if (!tdy->dm) { DM dm; - if (tdy->options.generate_mesh) { - // Here we lean on PETSc's DM* options. - ierr = DMCreate(comm, &dm); CHKERRQ(ierr); - ierr = DMSetType(dm, DMPLEX); CHKERRQ(ierr); - ierr = DMSetFromOptions(dm); CHKERRQ(ierr); - ierr = TDyDistributeDM(&dm); CHKERRQ(ierr); - } else { // if (tdy->options.read_mesh) + if (tdy->options.read_mesh) { ierr = DMPlexCreateFromFile(comm, tdy->options.mesh_file, PETSC_TRUE, &dm); CHKERRQ(ierr); + } else { + if (tdy->ops->create_dm) { + // We've been instructed to create a DM ourselves. + ierr = tdy->ops->create_dm(comm, tdy->create_dm_context, &dm); + } else if (create_dm_f90_) { + // Create a DM all-Fortran-like. + MPI_Fint comm_f = MPI_Comm_c2f(comm); + create_dm_f90_(&comm_f, &dm, &ierr); + } else { + ierr = DMPlexCreate(comm, &dm); CHKERRQ(ierr); + } + // Here we lean on PETSc's DM* options for overrides. + ierr = DMSetFromOptions(dm); CHKERRQ(ierr); } + // Distribute the mesh, however we got it. + ierr = TDyDistributeDM(&dm); CHKERRQ(ierr); tdy->dm = dm; } @@ -456,7 +494,7 @@ PetscErrorCode TDyDestroy(TDy *_tdy) { if (!tdy) PetscFunctionReturn(0); - ierr = TDyResetDiscretizationMethod(tdy); CHKERRQ(ierr); + ierr = TDyResetDiscretization(tdy); CHKERRQ(ierr); ierr = PetscFree(tdy->V); CHKERRQ(ierr); ierr = PetscFree(tdy->X); CHKERRQ(ierr); ierr = PetscFree(tdy->N); CHKERRQ(ierr); @@ -495,7 +533,7 @@ PetscErrorCode TDyGetDimension(TDy tdy,PetscInt *dim) { } /// Retrieves the DM used by the dycore. This must be called after -/// TDySetDM or TDySetFromOptions. +/// TDySetFromOptions. /// @param dm A pointer that stores the DM in use by the dycore PetscErrorCode TDyGetDM(TDy tdy,DM *dm) { PetscFunctionBegin; @@ -545,12 +583,15 @@ PetscErrorCode TDyGetCentroidArray(TDy tdy,PetscReal **X) { PetscFunctionReturn(0); } -PetscErrorCode TDyResetDiscretizationMethod(TDy tdy) { +PetscErrorCode TDyResetDiscretization(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidPointer(tdy,1); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + PetscInt dim; ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); @@ -572,7 +613,7 @@ PetscErrorCode TDyResetDiscretizationMethod(TDy tdy) { case 3: break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Unsupported dim in TDyResetDiscretizationMethod"); + SETERRQ(comm,PETSC_ERR_USER,"Unsupported dim in TDyResetDiscretization"); break; } // if (tdy->subc_Gmatrix) { ierr = TDyDeallocate_RealArray_4D(&tdy->subc_Gmatrix, tdy->mesh->num_cells, @@ -609,15 +650,16 @@ PetscErrorCode TDyView(TDy tdy,PetscViewer viewer) { /// Sets options for the dycore based on command line arguments supplied by a /// user. TDySetFromOptions must be called before TDySetup, /// since the latter uses options specified by the former. -/// @param tdy The dycore instance +/// @param [inout] tdy The dycore instance PetscErrorCode TDySetFromOptions(TDy tdy) { PetscErrorCode ierr; PetscFunctionBegin; - MPI_Comm comm = PETSC_COMM_WORLD; + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); if ((tdy->setup_flags & TDySetupFinished) != 0) { - SETERRQ(comm,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetup()"); + SETERRQ(comm,PETSC_ERR_USER,"You must call TDySetFromOptions before TDySetup()"); } // Collect options from command line arguments. @@ -649,10 +691,11 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Model options + TDyMode mode = options->mode; ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Model options",""); CHKERRQ(ierr); ierr = PetscOptionsEnum("-tdy_mode","Flow mode", "TDySetMode",TDyModes,(PetscEnum)options->mode, - (PetscEnum *)&options->mode, NULL); CHKERRQ(ierr); + (PetscEnum *)&mode, NULL); CHKERRQ(ierr); ierr = PetscOptionsReal("-tdy_gravity", "Magnitude of gravity vector", NULL, options->gravity_constant, &options->gravity_constant, NULL); CHKERRQ(ierr); @@ -688,10 +731,11 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Numerics options + TDyDiscretization discretization = options->discretization; ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Numerics options",""); CHKERRQ(ierr); - ierr = PetscOptionsEnum("-tdy_method","Discretization method", - "TDySetDiscretizationMethod",TDyMethods, - (PetscEnum)options->method,(PetscEnum *)&options->method, + ierr = PetscOptionsEnum("-tdy_discretization","Discretization", + "TDySetDiscretization",TDyDiscretizations, + (PetscEnum)options->discretization,(PetscEnum *)&discretization, NULL); CHKERRQ(ierr); TDyQuadratureType qtype = FULL; ierr = PetscOptionsEnum("-tdy_quadrature","Quadrature type for finite element methods","TDySetQuadratureType",TDyQuadratureTypes,(PetscEnum)qtype,(PetscEnum *)&options->qtype,NULL); CHKERRQ(ierr); @@ -713,34 +757,32 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { } // Mesh-related options - ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TDyCore: Mesh options",""); CHKERRQ(ierr); - ierr = PetscOptionsBool("-tdy_generate_mesh","Generate a mesh using provided PETSc DM options","",options->generate_mesh,&(options->generate_mesh),NULL); CHKERRQ(ierr); + ierr = PetscOptionsBegin(comm,NULL,"TDyCore: Mesh options",""); CHKERRQ(ierr); ierr = PetscOptionsGetString(NULL,NULL,"-tdy_read_mesh", options->mesh_file,sizeof(options->mesh_file),&options->read_mesh); CHKERRQ(ierr); - if (options->generate_mesh && options->read_mesh) { - SETERRQ(comm,PETSC_ERR_USER, - "Only one of -tdy_generate_mesh and -tdy_read_mesh can be specified"); - } - if ((tdy->dm != NULL) && (options->generate_mesh || options->read_mesh)) { - SETERRQ(comm,PETSC_ERR_USER, - "TDySetDM was called before TDySetFromOptions: can't generate or " - "read a mesh"); - } - if ((tdy->dm == NULL) && !(options->generate_mesh || options->read_mesh)) { - SETERRQ(comm,PETSC_ERR_USER, - "No mesh is available for TDycore: please use TDySetDM, " - "-tdy_generate_mesh, or -tdy_read_mesh to specify a mesh"); - } ierr = PetscOptionsBool("-tdy_output_mesh","Enable output of mesh attributes","",options->output_mesh,&(options->output_mesh),NULL); CHKERRQ(ierr); ierr = PetscOptionsGetString(NULL,NULL,"-tdy_output_geo_attributes", options->geom_attributes_file,sizeof(options->geom_attributes_file),&options->output_geom_attributes); CHKERRQ(ierr); ierr = PetscOptionsGetString(NULL,NULL,"-tdy_read_geo_attributes", options->geom_attributes_file,sizeof(options->geom_attributes_file),&options->read_geom_attributes); CHKERRQ(ierr); if (options->output_geom_attributes && options->read_geom_attributes){ - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Only one of -tdy_output_geom_attributes and -tdy_read_geom_attributes can be specified"); + SETERRQ(comm,PETSC_ERR_USER,"Only one of -tdy_output_geom_attributes and -tdy_read_geom_attributes can be specified"); } ierr = PetscOptionsEnd(); CHKERRQ(ierr); // Other options ierr = PetscOptionsBool("-tdy_regression_test","Enable output of a regression file","",options->regression_testing,&(options->regression_testing),NULL); CHKERRQ(ierr); + // Override the mode and/or discretization if needed. + if (options->mode != mode) { + TDySetMode(tdy, mode); + } + if (options->discretization != discretization) { + TDySetDiscretization(tdy, discretization); + } + + // Mode/discretization-specific options. + if (tdy->ops->set_from_options) { + tdy->ops->set_from_options(tdy); + } + // Wrap up and indicate that options are set. ierr = PetscOptionsEnd(); CHKERRQ(ierr); tdy->setup_flags |= TDyOptionsSet; @@ -760,61 +802,45 @@ PetscErrorCode TDySetFromOptions(TDy tdy) { PetscFunctionReturn(0); } -PetscErrorCode TDySetupDiscretizationScheme(TDy tdy) { - MPI_Comm comm; +/// Performs setup, including allocation and configuration of any bookkeeping +/// data structures and the configuration of the dycore's DM object, which can +/// subsequently be passed to a solver (e.g. SNES or TS). This function must be +/// called after TDySetOptions. +/// @param [inout] tdy the dycore instance +PetscErrorCode TDySetup(TDy tdy) { PetscErrorCode ierr; - PetscValidPointer(tdy,1); PetscFunctionBegin; - ierr = PetscObjectGetComm((PetscObject)(tdy->dm),&comm); CHKERRQ(ierr); - switch (tdy->options.method) { - case TPF: - ierr = TDyTPFInitialize(tdy); CHKERRQ(ierr); - break; - case MPFA_O: - ierr = TDyMPFAOInitialize(tdy); CHKERRQ(ierr); - break; - case MPFA_O_DAE: - ierr = TDyMPFAOInitialize(tdy); CHKERRQ(ierr); - break; - case MPFA_O_TRANSIENTVAR: - ierr = TDyMPFAOInitialize(tdy); CHKERRQ(ierr); - break; - case BDM: - ierr = TDyBDMInitialize(tdy); CHKERRQ(ierr); - break; - case WY: - ierr = TDyWYInitialize(tdy); CHKERRQ(ierr); - break; - } + TDY_START_FUNCTION_TIMER() - // Record metadata for scaling studies. - TDySetTimingMetadata(tdy); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); - PetscFunctionReturn(0); -} + if ((tdy->setup_flags & TDyDiscretizationSet) == 0) { + SETERRQ(comm,PETSC_ERR_USER, + "You must call TDySetDiscretization before TDySetup()"); + } -PetscErrorCode TDySetup(TDy tdy) { - /* must follow TDySetFromOptions() is it relies upon options set by - TDySetFromOptions */ - PetscErrorCode ierr; - PetscFunctionBegin; if ((tdy->setup_flags & TDyOptionsSet) == 0) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"TDySetFromOptions must be called prior to TDySetup()"); + SETERRQ(comm,PETSC_ERR_USER,"You must call TDySetFromOptions before TDySetup()"); } - TDY_START_FUNCTION_TIMER() TDyEnterProfilingStage("TDycore Setup"); - // TODO: Stick a call to tdy->ops->config_dm here to configure the DM for our - // TODO: specific dycore setup. - ierr = TDySetupDiscretizationScheme(tdy); CHKERRQ(ierr); + + // Perform implementation-specific setup. + // FIXME: Pass tdy->context as first argument when we break data out of TDy. + ierr = tdy->ops->setup(tdy, tdy->dm); CHKERRQ(ierr); + + // Record metadata for scaling studies. + TDySetTimingMetadata(tdy); + if (tdy->options.regression_testing) { /* must come after Sections are set up in - TDySetupDiscretizationScheme->XXXInitialize */ + TDySetupDiscretization->XXXInitialize */ ierr = TDyRegressionInitialize(tdy); CHKERRQ(ierr); } if (tdy->options.output_mesh) { - if (tdy->options.method != MPFA_O) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER, - "-tdy_output_mesh only supported for MPFA-O method"); + if (tdy->options.discretization != MPFA_O) { + SETERRQ(comm,PETSC_ERR_USER, + "-tdy_output_mesh only supported for MPFA-O discretization"); } } TDyExitProfilingStage("TDycore Setup"); @@ -823,16 +849,70 @@ PetscErrorCode TDySetup(TDy tdy) { PetscFunctionReturn(0); } -PetscErrorCode TDySetDiscretizationMethod(TDy tdy,TDyMethod method) { +/// Sets the dycore's "mode", specifying the governing equations it solves. +/// @param [inout] tdy the dycore +/// @param [in] mode the selected mode +PetscErrorCode TDySetMode(TDy tdy, TDyMode mode) { + PetscValidPointer(tdy,1); PetscFunctionBegin; - tdy->options.method = method; + tdy->options.mode = mode; + tdy->setup_flags |= TDyModeSet; + + // If we are resetting the mode and have already set the discretization, + // we call TDySetDiscretization again. + if (tdy->setup_flags & TDyDiscretizationSet) { + PetscInt ierr = TDySetDiscretization(tdy, tdy->options.discretization); CHKERRQ(ierr); + } + PetscFunctionReturn(0); } -PetscErrorCode TDySetMode(TDy tdy,TDyMode mode) { - PetscValidPointer(tdy,1); +/// Sets the discretization used by the dycore. This must be called after +/// TDySetMode. +/// @param [inout] tdy the dycore +/// @param [in] discretization the selected discretization +PetscErrorCode TDySetDiscretization(TDy tdy, TDyDiscretization discretization) { PetscFunctionBegin; - tdy->options.mode = mode; + + MPI_Comm comm; + PetscErrorCode ierr; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + + if ((tdy->setup_flags & TDyModeSet) == 0) { + SETERRQ(comm,PETSC_ERR_USER, + "You must call TDySetMode before TDySetDiscretization"); + } + + // Set function pointers for operations. + if (tdy->options.mode == RICHARDS) { + if (discretization == MPFA_O) { + tdy->ops->set_from_options = TDyMPFAO_SetFromOptions; + tdy->ops->setup = TDyRichards_MPFA_O_Setup; + } else if (discretization == MPFA_O_DAE) { + tdy->ops->set_from_options = TDyMPFAO_SetFromOptions; + tdy->ops->setup = TDyRichards_MPFA_O_DAE_Setup; + } else if (discretization == MPFA_O_TRANSIENTVAR) { + tdy->ops->set_from_options = TDyMPFAO_SetFromOptions; + tdy->ops->setup = TDyRichards_MPFA_O_TRANSIENTVAR_Setup; + } else if (discretization == BDM) { + tdy->ops->set_from_options = NULL; + tdy->ops->setup = TDyBDM_Setup; + } else if (discretization == WY) { + tdy->ops->set_from_options = NULL; + tdy->ops->setup = TDyWY_Setup; + } else { + SETERRQ(comm,PETSC_ERR_USER, "Invalid discretization given!"); + } + } else if (tdy->options.mode == TH) { + if (discretization == MPFA_O) { + tdy->ops->setup = TDyTH_MPFA_O_Setup; + } else { + SETERRQ(comm,PETSC_ERR_USER, + "The TH mode does not support the selected discretization!"); + } + } + tdy->options.discretization = discretization; + tdy->setup_flags |= TDyDiscretizationSet; PetscFunctionReturn(0); } @@ -887,10 +967,7 @@ PetscErrorCode TDySetIFunction(TS ts,TDy tdy) { ierr = PetscSectionGetNumFields(sec, &num_fields); ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); - switch (tdy->options.method) { - case TPF: - SETERRQ(comm,PETSC_ERR_SUP,"IFunction not implemented for TPF"); - break; + switch (tdy->options.discretization) { case MPFA_O: switch (tdy->options.mode) { case RICHARDS: @@ -927,10 +1004,7 @@ PetscErrorCode TDySetIJacobian(TS ts,TDy tdy) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)ts,&comm); CHKERRQ(ierr); - switch (tdy->options.method) { - case TPF: - SETERRQ(comm,PETSC_ERR_SUP,"IJacobian not implemented for TPF"); - break; + switch (tdy->options.discretization) { case MPFA_O: ierr = TDyCreateJacobian(tdy); CHKERRQ(ierr); switch (tdy->options.mode) { @@ -975,10 +1049,7 @@ PetscErrorCode TDySetSNESFunction(SNES snes,TDy tdy) { PetscInt dim; ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); - switch (tdy->options.method) { - case TPF: - SETERRQ(comm,PETSC_ERR_SUP,"SNESFunction not implemented for TPF"); - break; + switch (tdy->options.discretization) { case MPFA_O: ierr = SNESSetFunction(snes,tdy->residual,TDyMPFAOSNESFunction,tdy); CHKERRQ(ierr); break; @@ -1012,10 +1083,7 @@ PetscErrorCode TDySetSNESJacobian(SNES snes,TDy tdy) { PetscInt dim; ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); - switch (tdy->options.method) { - case TPF: - SETERRQ(comm,PETSC_ERR_SUP,"SNESJacobian not implemented for TPF"); - break; + switch (tdy->options.discretization) { case MPFA_O: ierr = SNESSetJacobian(snes,tdy->J,tdy->J,TDyMPFAOSNESJacobian,tdy); CHKERRQ(ierr); break; @@ -1042,10 +1110,7 @@ PetscErrorCode TDyComputeSystem(TDy tdy,Mat K,Vec F) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)(tdy->dm),&comm); CHKERRQ(ierr); - switch (tdy->options.method) { - case TPF: - ierr = TDyTPFComputeSystem(tdy,K,F); CHKERRQ(ierr); - break; + switch (tdy->options.discretization) { case MPFA_O: ierr = TDyMPFAOComputeSystem(tdy,K,F); CHKERRQ(ierr); break; @@ -1081,6 +1146,9 @@ PetscErrorCode TDyUpdateState(TDy tdy,PetscReal *U) { ierr = PetscMalloc((cEnd-cStart)*sizeof(PetscReal),&P);CHKERRQ(ierr); ierr = PetscMalloc((cEnd-cStart)*sizeof(PetscReal),&temp);CHKERRQ(ierr); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + if (tdy->options.mode == TH) { for (PetscInt c=0;cvg_m[c],alpha,cc->sr[i],tdy->Pref-P[i],&(cc->S[i]),&cc->dS_dP[i],&(cc->d2S_dP2[i])); break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unknown saturation function"); + SETERRQ(comm,PETSC_ERR_SUP,"Unknown saturation function"); break; } @@ -1123,7 +1191,7 @@ PetscErrorCode TDyUpdateState(TDy tdy,PetscReal *U) { RelativePermeability_Mualem(cc->mualem_m[c],cc->mualem_poly_low[c],cc->mualem_poly_coeffs[c],Se,&Kr,&dKr_dSe); break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unknown relative permeability function"); + SETERRQ(comm,PETSC_ERR_SUP,"Unknown relative permeability function"); break; } cc->Kr[i] = Kr; @@ -1143,9 +1211,9 @@ PetscErrorCode TDyUpdateState(TDy tdy,PetscReal *U) { } } - if ( (tdy->options.method == MPFA_O || - tdy->options.method == MPFA_O_DAE || - tdy->options.method == MPFA_O_TRANSIENTVAR)) { + if ( (tdy->options.discretization == MPFA_O || + tdy->options.discretization == MPFA_O_DAE || + tdy->options.discretization == MPFA_O_TRANSIENTVAR)) { PetscReal *p_vec_ptr, gz; TDyMesh *mesh = tdy->mesh; TDyCell *cells = &mesh->cells; @@ -1444,11 +1512,7 @@ PetscErrorCode TDyComputeErrorNorms(TDy tdy,Vec U,PetscReal *normp, PetscFunctionBegin; TDY_START_FUNCTION_TIMER() ierr = PetscObjectGetComm((PetscObject)(tdy->dm),&comm); CHKERRQ(ierr); - switch (tdy->options.method) { - case TPF: - if(normp != NULL) { *normp = TDyTPFPressureNorm(tdy,U); } - if(normv != NULL) { *normv = TDyTPFVelocityNorm(tdy,U); } - break; + switch (tdy->options.discretization) { case MPFA_O: if(normv) { ierr = TDyMPFAORecoverVelocity(tdy,U); CHKERRQ(ierr); @@ -1535,10 +1599,7 @@ PetscErrorCode TDyPreSolveSNESSolver(TDy tdy) { ierr = PetscObjectGetComm((PetscObject)tdy->dm,&comm); CHKERRQ(ierr); ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); - switch (tdy->options.method) { - case TPF: - SETERRQ(comm,PETSC_ERR_SUP,"TDyPreSolveSNESSolver not implemented for TPF"); - break; + switch (tdy->options.discretization) { case MPFA_O: ierr = TDyMPFAOSNESPreSolve(tdy); CHKERRQ(ierr); break; diff --git a/src/tdydriver.c b/src/tdydriver.c index 74f713b1..b5d3ecd5 100644 --- a/src/tdydriver.c +++ b/src/tdydriver.c @@ -16,21 +16,27 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { SNES snes; SNESLineSearch linesearch; + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + + DM dm; + ierr = TDyGetDM(tdy,&dm); CHKERRQ(ierr); + PetscInt dim; ierr = DMGetDimension(tdy->dm,&dim); CHKERRQ(ierr); + if (dim != 3) { - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Driver currently only supports 3D"); + SETERRQ(comm,PETSC_ERR_USER,"Driver currently only supports 3D"); } - switch(tdy->options.method) { - case TPF: + switch(tdy->options.discretization) { case MPFA_O: break; case MPFA_O_DAE: case MPFA_O_TRANSIENTVAR: case BDM: case WY: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Driver not supported for specified method."); + SETERRQ(comm,PETSC_ERR_USER,"Driver not supported for specified discretization."); break; } @@ -82,7 +88,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { case TH: break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Unrecognized flow mode."); + SETERRQ(comm,PETSC_ERR_USER,"Unrecognized flow mode."); } // check for unsupported time integration methods switch(tdy->ti->time_integration_method) { @@ -91,20 +97,20 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { case RICHARDS: break; case TH: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"SNES not supported for TH mode."); + SETERRQ(comm,PETSC_ERR_USER,"SNES not supported for TH mode."); break; } break; case TDyTS: break; default: - SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Unrecognized time integration method."); + SETERRQ(comm,PETSC_ERR_USER,"Unrecognized time integration method."); } // create time integrator switch(tdy->ti->time_integration_method) { case TDySNES: - ierr = SNESCreate(PETSC_COMM_WORLD,&snes); + ierr = SNESCreate(comm,&snes); CHKERRQ(ierr); ierr = TDySetSNESFunction(snes,tdy); CHKERRQ(ierr); ierr = TDySetSNESJacobian(snes,tdy); CHKERRQ(ierr); @@ -113,7 +119,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { tdy->ti->snes = snes; break; case TDyTS: - ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr); + ierr = TSCreate(comm,&ts); CHKERRQ(ierr); // ierr = TSSetType(ts,TSBEULER); CHKERRQ(ierr); // ierr = TSSetType(ts,TSPSEUDO); CHKERRQ(ierr); // ierr = TSPseudoSetTimeStep(ts,TSPseudoTimeStepDefault,NULL); CHKERRQ(ierr); @@ -148,6 +154,7 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { ierr = TSSetPostStep(ts,TDyRichardsTSPostStep); CHKERRQ(ierr); break; } + // FIXME: This is a different path from the one used by TDycore proper. ierr = TDyRichardsInitialize(tdy); CHKERRQ(ierr); break; case TH: @@ -161,10 +168,11 @@ PetscErrorCode TDyDriverInitializeTDy(TDy tdy) { case TDyTS: ierr = TSSetPostStep(ts,TDyTHTSPostStep); CHKERRQ(ierr); } + // FIXME: This is a different path from the one used by TDycore proper. ierr = TDyTHInitialize(tdy); CHKERRQ(ierr); break; } - PetscPrintf(PETSC_COMM_WORLD,"tdy->ti->time_integration_method = %d\n", + PetscPrintf(comm,"tdy->ti->time_integration_method = %d\n", tdy->ti->time_integration_method); size_t len; ierr = PetscStrlen(tdy->io->ic_filename, &len); CHKERRQ(ierr); diff --git a/src/tdyth.c b/src/tdyth.c index 09ed65ec..512d6357 100644 --- a/src/tdyth.c +++ b/src/tdyth.c @@ -9,11 +9,14 @@ PetscErrorCode TDyTHInitialize(TDy tdy) { PetscInt local_size; PetscFunctionBegin; - PetscPrintf(PETSC_COMM_WORLD,"Running TH mode.\n"); + MPI_Comm comm; + ierr = PetscObjectGetComm((PetscObject)tdy, &comm); CHKERRQ(ierr); + + PetscPrintf(comm,"Running TH mode.\n"); if (tdy->options.init_with_random_field) { PetscRandom rand; - ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand); CHKERRQ(ierr); + ierr = PetscRandomCreate(comm,&rand); CHKERRQ(ierr); ierr = VecGetLocalSize(tdy->solution,&local_size); CHKERRQ(ierr); ierr = DMCreateGlobalVector(tdy->dm,&temp_vec); CHKERRQ(ierr); // pressure diff --git a/src/tdytimers.c b/src/tdytimers.c index b50c5c8f..5bcfd8c3 100644 --- a/src/tdytimers.c +++ b/src/tdytimers.c @@ -11,8 +11,8 @@ khash_t(TDY_PROFILING_STAGE_MAP)* TDY_PROFILING_STAGES = NULL; // Timing metadata used by tdyperfplot and other tools. typedef struct { - TDyMethod method; TDyMode mode; + TDyDiscretization discretization; int num_cells; int num_proc; } TimingMetadata; @@ -91,8 +91,8 @@ PetscErrorCode TDySetTimingMetadata(TDy tdy) { iter = kh_put(TDY_PROFILING_MD_MAP, TDY_PROFILING_METADATA, tdy_addr, &retval); kh_val(TDY_PROFILING_METADATA, iter) = md; } - md->method = tdy->options.method; md->mode = tdy->options.mode; + md->discretization = tdy->options.discretization; if (tdy->mesh != NULL) { md->num_cells = TDyMeshGetNumberOfLocalCells(tdy->mesh); } else { @@ -136,31 +136,12 @@ PetscErrorCode TDyWriteTimingProfile(const char* filename) { int num_cells; MPI_Reduce(&(md->num_cells), &num_cells, 1, MPI_INT, MPI_SUM, 0, PETSC_COMM_WORLD); - const char* method_name; - if (md->method == TPF) { - method_name = "TPF"; - } else if (md->method == MPFA_O) { - method_name = "MPFA_O"; - } else if (md->method == MPFA_O_DAE) { - method_name = "MPFA_O_DAE"; - } else if (md->method == MPFA_O_TRANSIENTVAR) { - method_name = "MPFA_O_TRANSIENTVAR"; - } else if (md->method == BDM) { - method_name = "BDM"; - } else { // (md->method == BDM) - method_name = "WY"; - } - const char* mode_name; - if (md->mode == RICHARDS) { - mode_name = "RICHARDS"; - } else { // (md->mode == TH) - mode_name = "TH"; - } + const char* mode_name = TDyModes[md->mode]; + const char* disc_name = TDyDiscretizations[md->discretization]; FILE* f = fopen("tdycore_profile.csv", "a"); fprintf(f, "METADATA\n"); - fprintf(f, "Method,Mode,NumProc,NumCells\n"); - fprintf(f, "%s,%s,%d,%d", method_name, mode_name, - md->num_proc, num_cells); + fprintf(f, "Mode,Discretization,NumProc,NumCells\n"); + fprintf(f, "%s,%s,%d,%d", mode_name, disc_name, md->num_proc, num_cells); fclose(f); } else { // rank > 0 MPI_Reduce(&md->num_cells, NULL, 1, MPI_INT, MPI_SUM, diff --git a/src/tdytpf.c b/src/tdytpf.c deleted file mode 100644 index da0509bf..00000000 --- a/src/tdytpf.c +++ /dev/null @@ -1,281 +0,0 @@ -#include -#include - -PETSC_STATIC_INLINE void Waxpy(PetscInt dim, PetscScalar a, - const PetscScalar *x, const PetscScalar *y, PetscScalar *w) {PetscInt d; for (d = 0; d < dim; ++d) w[d] = a*x[d] + y[d];} -PETSC_STATIC_INLINE PetscScalar Dot(PetscInt dim, const PetscScalar *x, - const PetscScalar *y) {PetscScalar sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*y[d]; return sum;} -PETSC_STATIC_INLINE PetscReal Norm(PetscInt dim, const PetscScalar *x) {return PetscSqrtReal(PetscAbsScalar(Dot(dim,x,x)));} - -PetscErrorCode TDyTPFInitialize(TDy tdy) { - PetscErrorCode ierr; - MPI_Comm comm; - PetscInt dim,c,cStart,cEnd,pStart,pEnd; - PetscSection sec; - DM dm = tdy->dm; - PetscFunctionBegin; - TDY_START_FUNCTION_TIMER() - - ierr = PetscObjectGetComm((PetscObject)dm,&comm); CHKERRQ(ierr); - ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - ierr = DMPlexGetChart(dm,&pStart,&pEnd); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd); CHKERRQ(ierr); - - /* Setup the section, 1 dof per cell */ - ierr = PetscSectionCreate(comm,&sec); CHKERRQ(ierr); - ierr = PetscSectionSetNumFields(sec,1); CHKERRQ(ierr); - ierr = PetscSectionSetFieldName(sec,0,"Pressure"); CHKERRQ(ierr); - ierr = PetscSectionSetFieldComponents(sec,0,1); CHKERRQ(ierr); - ierr = DMPlexGetChart(dm,&pStart,&pEnd); CHKERRQ(ierr); - ierr = PetscSectionSetChart(sec,pStart,pEnd); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,0,&pStart,&pEnd); CHKERRQ(ierr); - for(c=cStart; cdm; - PetscFunctionBegin; - TDY_START_FUNCTION_TIMER() - if(!tdy->options.tpf_allow_all_meshes) { - ierr = TDyTPFCheckMeshSuitability(tdy); CHKERRQ(ierr); - } - ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - dim2 = dim*dim; - ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd); CHKERRQ(ierr); - MaterialProp *matprop = tdy->matprop; - - // Here's our domain boundary label. - DMLabel bnd_label; - ierr = DMGetLabel(tdy->dm, "boundary", &bnd_label); CHKERRQ(ierr); - - // Loop over all faces and compute fluxes. - for (PetscInt f = fStart; f < fEnd; ++f) { - - const PetscInt *supp; - ierr = DMPlexGetSupport(dm,f,&supp); CHKERRQ(ierr); - - // wrt - // v = -Ki (pd-p0)/dist - // 00 Ki/dist - // 0 -Ki pd / dist - - if (tdy->ops->compute_boundary_pressure) { - // Enforce boundary pressures if needed. Boundary faces have a label - // value of 1. - PetscInt bnd_label_val; - ierr = DMLabelGetValue(bnd_label, f, &bnd_label_val); - if (bnd_label_val == 1) { - - ierr = (*tdy->ops->compute_boundary_pressure)(tdy, &(tdy->X[f*dim]),&p, tdy->boundary_pressure_ctx);CHKERRQ(ierr); - if (fabs(p+999.)<1.e-40) continue; - ierr = DMPlexGetPointGlobal(dm,supp[0],&row,&junk); CHKERRQ(ierr); - Waxpy(dim,-1,&(tdy->X[supp[0]*dim]),&(tdy->X[f*dim]),pnt2pnt); - dist = Norm(dim,pnt2pnt); - Ki = matprop->K[(supp[0]-cStart)*dim2]; - ierr = MatSetValue(K,row,row,Ki/dist*tdy->V[f]/tdy->V[supp[0]],ADD_VALUES); - CHKERRQ(ierr); - ierr = VecSetValue(F,row,p*Ki/dist*tdy->V[f]/tdy->V[supp[0]],ADD_VALUES); - CHKERRQ(ierr); - continue; - } - } - - Waxpy(dim,-1,&(tdy->X[supp[0]*dim]),&(tdy->X[supp[1]*dim]),pnt2pnt); - dist = Norm(dim,pnt2pnt); - Ki = 0.5*(matprop->K[(supp[0]-cStart)*dim2]+matprop->K[(supp[1]-cStart)*dim2]); - - //wrt 0 - // v = -Ki (p1-p0)/dist - // 00 Ki/dist - // 01 -Ki/dist - - //wrt 1 - // v = -Ki (p0-p1)/dist - // 11 Ki/dist - // 10 -Ki/dist - - ierr = DMPlexGetPointGlobal(dm,supp[0],&row,&junk); CHKERRQ(ierr); - ierr = DMPlexGetPointGlobal(dm,supp[1],&col,&junk); CHKERRQ(ierr); - ierr = MatSetValue(K,row,row, Ki/dist*tdy->V[f]/tdy->V[supp[0]],ADD_VALUES); - CHKERRQ(ierr); - ierr = MatSetValue(K,row,col,-Ki/dist*tdy->V[f]/tdy->V[supp[0]],ADD_VALUES); - CHKERRQ(ierr); - ierr = MatSetValue(K,col,col, Ki/dist*tdy->V[f]/tdy->V[supp[1]],ADD_VALUES); - CHKERRQ(ierr); - ierr = MatSetValue(K,col,row,-Ki/dist*tdy->V[f]/tdy->V[supp[1]],ADD_VALUES); - CHKERRQ(ierr); - } - - if(tdy->ops->computeforcing) { - for(c=cStart; cops->computeforcing)(tdy,&(tdy->X[c*dim]),&force,tdy->forcingctx);CHKERRQ(ierr);; - ierr = VecSetValue(F,row,force,ADD_VALUES); CHKERRQ(ierr); - } - } - - // max = 3.63988 - // min = 3.15550 - - ierr = VecAssemblyBegin(F); CHKERRQ(ierr); - ierr = VecAssemblyEnd (F); CHKERRQ(ierr); - ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); - ierr = MatAssemblyEnd (K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); - - TDY_STOP_FUNCTION_TIMER() - PetscFunctionReturn(0); -} - -/* - norm = sqrt( sum_i^n V_i * ( p(X_i) - P_i )^2 ) - - where n is the number of cells. - */ -PetscReal TDyTPFPressureNorm(TDy tdy,Vec U) { - PetscFunctionBegin; - TDY_START_FUNCTION_TIMER() - PetscErrorCode ierr; - PetscSection sec; - PetscInt c,cStart,cEnd,offset,dim,gref,junk; - PetscReal p,*u,norm,norm_sum; - DM dm = tdy->dm; - if(!(tdy->ops->compute_boundary_pressure)) { - SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_USER, - "You must set the pressure function with TDySetBoundaryPressureFn"); - } - norm = 0; - ierr = VecGetArray(U,&u); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd); CHKERRQ(ierr); - ierr = DMGetSection(dm,&sec); CHKERRQ(ierr); - ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - for(c=cStart; cops->compute_boundary_pressure)(tdy,&(tdy->X[c*dim]),&p,tdy->boundary_pressure_ctx);CHKERRQ(ierr); - norm += tdy->V[c]*PetscSqr(u[offset]-p); - } - ierr = MPI_Allreduce(&norm,&norm_sum,1,MPIU_REAL,MPI_SUM, - PetscObjectComm((PetscObject)U)); CHKERRQ(ierr); - norm_sum = PetscSqrtReal(norm_sum); - ierr = VecRestoreArray(U,&u); CHKERRQ(ierr); - TDY_STOP_FUNCTION_TIMER() - PetscFunctionReturn(norm_sum); -} - -PetscReal TDyTPFVelocityNorm(TDy tdy,Vec U) { - PetscErrorCode ierr; - PetscInt dim,dim2,fStart,fEnd,cStart,cEnd,row,junk; - PetscReal pnt2pnt[3],dist,Ki,p,vel[3],va,ve,*u,sign,face_error,norm,norm_sum; - DM dm = tdy->dm; - PetscFunctionBegin; - TDY_START_FUNCTION_TIMER() - ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - dim2 = dim*dim; - ierr = VecGetArray(U,&u); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd); CHKERRQ(ierr); - face_error = 0; - norm = 0; - MaterialProp *matprop = tdy->matprop; - - // Here's our domain boundary label. - DMLabel bnd_label; - ierr = DMGetLabel(tdy->dm, "boundary", &bnd_label); CHKERRQ(ierr); - - for(PetscInt c=cStart; cops->compute_boundary_pressure) { - ierr = DMPlexGetPointGlobal(dm,supp[0],&row,&junk); CHKERRQ(ierr); - Waxpy(dim,-1,&(tdy->X[supp[0]*dim]),&(tdy->X[f*dim]),pnt2pnt); - dist = Norm(dim,pnt2pnt); - Ki = matprop->K[(supp[0]-cStart)*dim2]; - ierr = (*tdy->ops->compute_boundary_pressure)(tdy, &(tdy->X[f*dim]),&p, tdy->boundary_pressure_ctx);CHKERRQ(ierr); - sign = PetscSign(TDyADotBMinusC(&(tdy->N[dim*f]),&(tdy->X[dim*f]), - &(tdy->X[dim*supp[0]]),dim)); - va = sign* -Ki*(p-u[(supp[0]-cStart)])/dist; - ierr = (*tdy->ops->compute_boundary_velocity)(tdy,&(tdy->X[f*dim]),&(vel[0]),tdy->boundary_velocity_ctx);CHKERRQ(ierr); - ve = TDyADotB(vel,&(tdy->N[dim*f]),dim); - face_error = PetscSqr((va-ve)/tdy->V[f]); - } else { - Waxpy(dim,-1,&(tdy->X[supp[0]*dim]),&(tdy->X[supp[1]*dim]),pnt2pnt); - dist = Norm(dim,pnt2pnt); - Ki = 0.5*(matprop->K[(supp[0]-cStart)*dim2]+matprop->K[(supp[1]-cStart)*dim2]); - sign = PetscSign(TDyADotBMinusC(&(tdy->N[dim*f]),&(tdy->X[dim*f]), - &(tdy->X[dim*c]),dim)); - va = sign* -Ki*(u[(supp[1]-cStart)]-u[(supp[0]-cStart)])/dist *tdy->V[f]; - if(c==supp[1]) va *= -1; - ierr = (*tdy->ops->compute_boundary_velocity)(tdy,&(tdy->X[f*dim]),&(vel[0]),tdy->boundary_velocity_ctx);CHKERRQ(ierr); - ve = TDyADotB(vel,&(tdy->N[dim*f]),dim)*tdy->V[f]; - face_error = PetscSqr((va-ve)/tdy->V[f]); - } - norm += face_error*tdy->V[c]; - } - } - ierr = MPI_Allreduce(&norm,&norm_sum,1,MPIU_REAL,MPI_SUM, - PetscObjectComm((PetscObject)dm)); CHKERRQ(ierr); - norm_sum = PetscSqrtReal(norm_sum); - ierr = VecRestoreArray(U,&u); CHKERRQ(ierr); - TDY_STOP_FUNCTION_TIMER() - PetscFunctionReturn(norm_sum); -} - -PetscErrorCode TDyTPFCheckMeshSuitability(TDy tdy) { - PetscFunctionBegin; - TDY_START_FUNCTION_TIMER() - PetscErrorCode ierr; - PetscInt dim,f,fStart,fEnd; - PetscReal diff,dist,pnt2pnt[3]; - DM dm = tdy->dm; - ierr = DMGetDimension(dm,&dim); CHKERRQ(ierr); - ierr = DMPlexGetHeightStratum(dm,1,&fStart,&fEnd); CHKERRQ(ierr); - for(f=fStart; fX[supp[1]*dim]),&(tdy->X[supp[0]*dim]),pnt2pnt); - dist = Norm(dim,pnt2pnt); - diff = PetscAbsReal(TDyADotB(&(tdy->N[f*dim]),pnt2pnt,dim)/dist); - if(PetscAbsReal(diff-1) > 10*PETSC_MACHINE_EPSILON) { - SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP, - "Mesh is unsuitable for a consistent two point flux method. To force rerun with -tdy_tpf_allow_unsuitable_mesh"); - } - } - TDY_STOP_FUNCTION_TIMER() - PetscFunctionReturn(0); -} diff --git a/src/tdywy.c b/src/tdywy.c index 7b80db00..aa9d3802 100644 --- a/src/tdywy.c +++ b/src/tdywy.c @@ -217,7 +217,7 @@ PetscErrorCode TDyWYLocalElementCompute(TDy tdy) { PetscFunctionReturn(0); } -PetscErrorCode TDyWYInitialize(TDy tdy) { +PetscErrorCode TDyWY_Setup(TDy tdy, DM dm) { PetscFunctionBegin; TDY_START_FUNCTION_TIMER() PetscErrorCode ierr; @@ -225,7 +225,6 @@ PetscErrorCode TDyWYInitialize(TDy tdy) { PetscInt i,dim,ncv,nfv,c,cStart,cEnd,f,fStart,fEnd,vStart,vEnd,p,pStart,pEnd; PetscInt closureSize, *closure; PetscSection sec; - DM dm = tdy->dm; MaterialProp *matprop = tdy->matprop; ierr = PetscObjectGetComm((PetscObject)dm,&comm); CHKERRQ(ierr);