diff --git a/example/test_problem/ELBDM/Soliton/Input__TestProb b/example/test_problem/ELBDM/Soliton/Input__TestProb index 889cb6ff5f..0cc0967910 100644 --- a/example/test_problem/ELBDM/Soliton/Input__TestProb +++ b/example/test_problem/ELBDM/Soliton/Input__TestProb @@ -5,4 +5,5 @@ Soliton_CoreRadiusAll 1.0e1 # core radius for ALL solitons (in code uni Soliton_EmptyRegion 0.0 # soliton-free region [0.0] # -> solitons only reside in [Soliton_EmptyRegion : BOX_SIZE-Soliton_EmptyRegion] # -> useful only when Soliton_RSeed >= 0 +Soliton_InitMode 1 # soliton initialization mode (1=table of the density profile, 2=analytical function of the density profile) [1] Soliton_DensProf_Filename SolitonDensityProfile_Lambda0.0 # filename of the input soliton density profile diff --git a/example/test_problem/ELBDM/Soliton/README b/example/test_problem/ELBDM/Soliton/README index 74960f341b..1ad35d9af0 100644 --- a/example/test_problem/ELBDM/Soliton/README +++ b/example/test_problem/ELBDM/Soliton/README @@ -18,3 +18,4 @@ Note: Soliton_N > 1: soliton merger --> Set Soliton_RSeed >= 0 Set Soliton_EmptyRegion > 0.0 (e.g., 1.0e2) +2. Input density profile tables should use $| \psi |^2$ units. The code automatically converts this to mass density. \ No newline at end of file diff --git a/src/TestProblem/ELBDM/Soliton/Init_TestProb_ELBDM_Soliton.cpp b/src/TestProblem/ELBDM/Soliton/Init_TestProb_ELBDM_Soliton.cpp index 2f384ab4be..afeaf914d7 100644 --- a/src/TestProblem/ELBDM/Soliton/Init_TestProb_ELBDM_Soliton.cpp +++ b/src/TestProblem/ELBDM/Soliton/Init_TestProb_ELBDM_Soliton.cpp @@ -11,6 +11,7 @@ static double Soliton_CoreRadiusAll; // core radius for all // (<=0.0 --> hard coding each soliton) static double Soliton_EmptyRegion; // soliton-free region from the boundary // (useful only when Soliton_RSeed>=0) +static int Soliton_InitMode; // soliton initialization mode (1=table of the density profile, 2=analytical function of the density profile) static char Soliton_DensProf_Filename[MAX_STRING]; // filename of the reference soliton density profile static int Soliton_DensProf_NBin; // number of radial bins of the soliton density profile @@ -20,7 +21,7 @@ static double (*Soliton_Center)[3] = NULL; // center coordinates o static double *Soliton_ScaleL = NULL; // L/D: length/density scale factors of each soliton // (defined as the ratio between the core radii/peak // density of the target and reference soliton profiles) -static double *Soliton_ScaleD = NULL; +static double *Soliton_ScaleD = NULL; // ( for Soliton_InitMode==1 only) // ======================================================================================= static void BC( real Array[], const int ArraySize[], real fluid[], const int NVar_Flu, @@ -124,6 +125,7 @@ void LoadInputTestProb( const LoadParaMode_t load_mode, ReadPara_t *ReadPara, HD LOAD_PARA( load_mode, "Soliton_RSeed", &Soliton_RSeed, 0, NoMin_int, NoMax_int ); LOAD_PARA( load_mode, "Soliton_CoreRadiusAll", &Soliton_CoreRadiusAll, NoDef_double, NoMin_double, NoMax_double ); LOAD_PARA( load_mode, "Soliton_EmptyRegion", &Soliton_EmptyRegion, 0.0, NoMin_double, NoMax_double ); + LOAD_PARA( load_mode, "Soliton_InitMode", &Soliton_InitMode, 1, 1, 2 ); LOAD_PARA( load_mode, "Soliton_DensProf_Filename", Soliton_DensProf_Filename, NoDef_str, Useless_str, Useless_str ); } // FUNCITON : LoadInputTestProb @@ -171,13 +173,21 @@ void SetParameter() if ( Soliton_CoreRadiusAll == NoDef_double ) Aux_Error( ERROR_INFO, "Runtime parameter \"Soliton_CoreRadiusAll\" is not set !!\n" ); + if ( !OPT__UNIT && Soliton_InitMode == 2 ) + Aux_Error( ERROR_INFO, "OPT__UNIT must be enabled for Soliton_InitMode == 2 !!\n" ); // (2) set the problem-specific derived parameters // (2-1) allocate memory Soliton_CoreRadius = new double [Soliton_N]; Soliton_Center = new double [Soliton_N][3]; - Soliton_ScaleL = new double [Soliton_N]; - Soliton_ScaleD = new double [Soliton_N]; + if ( Soliton_InitMode == 1 ) + { + Soliton_ScaleL = new double [Soliton_N]; + Soliton_ScaleD = new double [Soliton_N]; + } + else if ( Soliton_InitMode == 2 ) {} + else + Aux_Error( ERROR_INFO, "unsupported Soliton_InitMode (%d) !!\n", Soliton_InitMode ); // (2-2) soliton core radii if ( Soliton_CoreRadiusAll > 0.0 ) @@ -231,41 +241,47 @@ void SetParameter() // (3) load the reference soliton density profile and evaluate the scale factors if ( OPT__INIT != INIT_BY_RESTART ) { -// load the reference profile - const bool RowMajor_No = false; // load data into the column-major order - const bool AllocMem_Yes = true; // allocate memory for Soliton_DensProf - const int NCol = 2; // total number of columns to load - const int Col[NCol] = {0, 1}; // target columns: (radius, density) + if ( Soliton_InitMode == 1 ) + { + // load the reference profile + const bool RowMajor_No = false; // load data into the column-major order + const bool AllocMem_Yes = true; // allocate memory for Soliton_DensProf + const int NCol = 2; // total number of columns to load + const int Col[NCol] = {0, 1}; // target columns: (radius, density) - Soliton_DensProf_NBin = Aux_LoadTable( Soliton_DensProf, Soliton_DensProf_Filename, NCol, Col, RowMajor_No, AllocMem_Yes ); + Soliton_DensProf_NBin = Aux_LoadTable( Soliton_DensProf, Soliton_DensProf_Filename, NCol, Col, RowMajor_No, AllocMem_Yes ); -// get the core radius of the reference profile - const double *RadiusRef = Soliton_DensProf + 0*Soliton_DensProf_NBin; - const double *DensRef = Soliton_DensProf + 1*Soliton_DensProf_NBin; - const double DensCore = 0.5*DensRef[0]; // define core radius as the half-density radius + // get the core radius of the reference profile + const double *RadiusRef = Soliton_DensProf + 0*Soliton_DensProf_NBin; + const double *DensRef = Soliton_DensProf + 1*Soliton_DensProf_NBin; + const double DensCore = 0.5*DensRef[0]; // define core radius as the half-density radius - double CoreRadiusRef = NULL_REAL; + double CoreRadiusRef = NULL_REAL; - for (int b=1; b= DensCore && DensRef[b+1] <= DensCore ) + for (int b=1; b= DensCore && DensRef[b+1] <= DensCore ) + { + CoreRadiusRef = RadiusRef[b] + (DensCore - DensRef[b])*(RadiusRef[b+1]-RadiusRef[b])/(DensRef[b+1]-DensRef[b]); + break; + } } - } - if ( CoreRadiusRef == NULL_REAL ) - Aux_Error( ERROR_INFO, "cannot determine the reference core radius !!\n" ); + if ( CoreRadiusRef == NULL_REAL ) + Aux_Error( ERROR_INFO, "cannot determine the reference core radius !!\n" ); -// evaluate the scale factors of each soliton - for (int t=0; t use the table of the density profile)\n" ); + Aux_Message( stdout, " %7s %13s %13s %13s %13s %13s %13s\n", + "ID", "CoreRadius", "ScaleL", "ScaleD", "Center_X", "Center_Y", "Center_Z" ); + for (int t=0; t use the analytical function of the density profile)\n" ); + Aux_Message( stdout, " %7s %13s %13s %13s %13s\n", "ID", "CoreRadius", "Center_X", "Center_Y", "Center_Z" ); + for (int t=0; t reference radius) - r_ref = r_tar / Soliton_ScaleL[t]; + // loop over all solitons to get the total density + for (int t=0; t reference radius) + r_ref = r_tar / Soliton_ScaleL[t]; - if ( dens_ref == NULL_REAL ) - { - if ( r_ref < Table_Radius[0] ) - dens_ref = Table_Density[0]; + // linear interpolation + dens_ref = Mis_InterpolateFromTable( Soliton_DensProf_NBin, Table_Radius, Table_Density, r_ref ); - else if ( r_ref >= Table_Radius[Soliton_DensProf_NBin-1] ) - dens_ref = Table_Density[Soliton_DensProf_NBin-1]; + if ( dens_ref == NULL_REAL ) + { + if ( r_ref < Table_Radius[0] ) + dens_ref = Table_Density[0]; - else - Aux_Error( ERROR_INFO, "interpolation failed at radius %13.7e (min/max radius = %13.7e/%13.7e) !!\n", - r_ref, Table_Radius[0], Table_Radius[Soliton_DensProf_NBin-1] ); - } + else if ( r_ref >= Table_Radius[Soliton_DensProf_NBin-1] ) + dens_ref = Table_Density[Soliton_DensProf_NBin-1]; -// rescale density (reference density --> target density) and add to the fluid array - fluid[DENS] += dens_ref*Soliton_ScaleD[t]; - } // for (int t=0; t target density) and add to the fluid array + fluid[DENS] += dens_ref*Soliton_ScaleD[t]; + } // for (int t=0; t