diff --git a/environ/param.gpu.water b/environ/param.gpu.water
index 2982676bc96c9aad99bd5fde9e8562fb90c68a81..0fdd8726f0cf57906762dfba1c2dc2185c7b83c1 100644
--- a/environ/param.gpu.water
+++ b/environ/param.gpu.water
@@ -21,6 +21,8 @@ cm_solver_max_iters 20 ! max solver iterations
cm_solver_restart 100 ! inner iterations of GMRES before restarting
cm_solver_q_err 1e-6 ! relative residual norm threshold used in solver
cm_domain_sparsity 1.0 ! scalar for scaling cut-off distance, used to sparsify charge matrix (between 0.0 and 1.0)
+cm_init_guess_extrap1 3 ! order of spline extrapolation for initial guess (s)
+cm_init_guess_extrap2 2 ! order of spline extrapolation for initial guess (t)
cm_solver_pre_comp_type 1 ! method used to compute preconditioner, if applicable
cm_solver_pre_comp_refactor 1000 ! number of steps before recomputing preconditioner
cm_solver_pre_comp_droptol 0.0 ! threshold tolerance for dropping values in preconditioner computation, if applicable
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 5d2faa75397beec670a3dedc38c758d02b4d5a2e..b61db4b85bc42ee5767f4d74aef0c2f1f5f7edf8 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -39,7 +39,7 @@ static void Extrapolate_Charges_QEq( const reax_system * const system,
int i;
real s_tmp, t_tmp;
- /* extrapolation for s & t */
+ /* spline extrapolation for s & t */
//TODO: good candidate for vectorization, avoid moving data with head pointer and circular buffer
#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
@@ -47,29 +47,61 @@ static void Extrapolate_Charges_QEq( const reax_system * const system,
#endif
for ( i = 0; i < system->N_cm; ++i )
{
- // no extrapolation
- //s_tmp = workspace->s[0][i];
- //t_tmp = workspace->t[0][i];
-
- // linear
- //s_tmp = 2 * workspace->s[0][i] - workspace->s[1][i];
- //t_tmp = 2 * workspace->t[0][i] - workspace->t[1][i];
-
- // quadratic
- //s_tmp = workspace->s[2][i] + 3 * (workspace->s[0][i]-workspace->s[1][i]);
- t_tmp = workspace->t[2][i] + 3 * (workspace->t[0][i] - workspace->t[1][i]);
-
- // cubic
- s_tmp = 4 * (workspace->s[0][i] + workspace->s[2][i]) -
- (6 * workspace->s[1][i] + workspace->s[3][i] );
- //t_tmp = 4 * (workspace->t[0][i] + workspace->t[2][i]) -
- // (6 * workspace->t[1][i] + workspace->t[3][i] );
-
- // 4th order
- //s_tmp = 5 * (workspace->s[0][i] - workspace->s[3][i]) +
- // 10 * (-workspace->s[1][i] + workspace->s[2][i] ) + workspace->s[4][i];
- //t_tmp = 5 * (workspace->t[0][i] - workspace->t[3][i]) +
- // 10 * (-workspace->t[1][i] + workspace->t[2][i] ) + workspace->t[4][i];
+ /* no extrapolation */
+ if ( control->cm_init_guess_extrap1 == 0 )
+ {
+ s_tmp = workspace->s[0][i];
+ }
+ /* linear */
+ else if ( control->cm_init_guess_extrap1 == 1 )
+ {
+ s_tmp = 2 * workspace->s[0][i] - workspace->s[1][i];
+ }
+ /* quadratic */
+ else if ( control->cm_init_guess_extrap1 == 2 )
+ {
+ s_tmp = workspace->s[2][i] + 3 * (workspace->s[0][i]-workspace->s[1][i]);
+ }
+ /* cubic */
+ else if ( control->cm_init_guess_extrap1 == 3 )
+ {
+ s_tmp = 4 * (workspace->s[0][i] + workspace->s[2][i]) -
+ (6 * workspace->s[1][i] + workspace->s[3][i] );
+ }
+ /* 4th order */
+ else if ( control->cm_init_guess_extrap1 == 4 )
+ {
+ s_tmp = 5 * (workspace->s[0][i] - workspace->s[3][i]) +
+ 10 * (-workspace->s[1][i] + workspace->s[2][i] ) + workspace->s[4][i];
+ }
+
+ /* no extrapolation */
+ if ( control->cm_init_guess_extrap2 == 0 )
+ {
+ t_tmp = workspace->t[0][i];
+ }
+ /* linear */
+ else if ( control->cm_init_guess_extrap2 == 1 )
+ {
+ t_tmp = 2.0 * workspace->t[0][i] - workspace->t[1][i];
+ }
+ /* quadratic */
+ else if ( control->cm_init_guess_extrap2 == 2 )
+ {
+ t_tmp = workspace->t[2][i] + 3 * (workspace->t[0][i] - workspace->t[1][i]);
+ }
+ /* cubic */
+ else if ( control->cm_init_guess_extrap2 == 3 )
+ {
+ t_tmp = 4.0 * (workspace->t[0][i] + workspace->t[2][i]) -
+ (6.0 * workspace->t[1][i] + workspace->t[3][i] );
+ }
+ /* 4th order */
+ else if ( control->cm_init_guess_extrap2 == 4 )
+ {
+ t_tmp = 5.0 * (workspace->t[0][i] - workspace->t[3][i]) +
+ 10.0 * (-workspace->t[1][i] + workspace->t[2][i] ) + workspace->t[4][i];
+ }
workspace->s[4][i] = workspace->s[3][i];
workspace->s[3][i] = workspace->s[2][i];
@@ -93,7 +125,7 @@ static void Extrapolate_Charges_EE( const reax_system * const system,
int i;
real s_tmp;
- /* extrapolation for s */
+ /* spline extrapolation for s */
//TODO: good candidate for vectorization, avoid moving data with head pointer and circular buffer
#ifdef _OPENMP
#pragma omp parallel for schedule(static) \
@@ -101,22 +133,33 @@ static void Extrapolate_Charges_EE( const reax_system * const system,
#endif
for ( i = 0; i < system->N_cm; ++i )
{
- // no extrapolation
- //s_tmp = workspace->s[0][i];
-
- // linear
- //s_tmp = 2 * workspace->s[0][i] - workspace->s[1][i];
-
- // quadratic
- //s_tmp = workspace->s[2][i] + 3 * (workspace->s[0][i]-workspace->s[1][i]);
-
- // cubic
- s_tmp = 4 * (workspace->s[0][i] + workspace->s[2][i]) -
- (6 * workspace->s[1][i] + workspace->s[3][i] );
-
- // 4th order
- //s_tmp = 5 * (workspace->s[0][i] - workspace->s[3][i]) +
- // 10 * (-workspace->s[1][i] + workspace->s[2][i] ) + workspace->s[4][i];
+ /* no extrapolation */
+ if ( control->cm_init_guess_extrap1 == 0 )
+ {
+ s_tmp = workspace->s[0][i];
+ }
+ /* linear */
+ else if ( control->cm_init_guess_extrap1 == 1 )
+ {
+ s_tmp = 2.0 * workspace->s[0][i] - workspace->s[1][i];
+ }
+ /* quadratic */
+ else if ( control->cm_init_guess_extrap1 == 2 )
+ {
+ s_tmp = workspace->s[2][i] + 3.0 * (workspace->s[0][i]-workspace->s[1][i]);
+ }
+ /* cubic */
+ else if ( control->cm_init_guess_extrap1 == 3 )
+ {
+ s_tmp = 4.0 * (workspace->s[0][i] + workspace->s[2][i]) -
+ (6.0 * workspace->s[1][i] + workspace->s[3][i] );
+ }
+ /* 4th order */
+ else if ( control->cm_init_guess_extrap1 == 4 )
+ {
+ s_tmp = 5.0 * (workspace->s[0][i] - workspace->s[3][i]) +
+ 10.0 * (-workspace->s[1][i] + workspace->s[2][i] ) + workspace->s[4][i];
+ }
workspace->s[4][i] = workspace->s[3][i];
workspace->s[3][i] = workspace->s[2][i];
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index bec04123391c1a165a1252d8d7ccb66a8deb37c0..5325a929dfa08855e406d1dc6c3dbf5a0b453326 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -76,6 +76,8 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->cm_solver_q_err = 0.000001;
control->cm_domain_sparsify_enabled = FALSE;
control->cm_domain_sparsity = 1.0;
+ control->cm_init_guess_extrap1 = 3;
+ control->cm_init_guess_extrap2 = 2;
control->cm_solver_pre_comp_type = ICHOLT_PC;
control->cm_solver_pre_comp_sweeps = 3;
control->cm_solver_pre_comp_sai_thres = 0.1;
@@ -276,6 +278,16 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->cm_domain_sparsify_enabled = TRUE;
}
}
+ else if ( strcmp(tmp[0], "cm_init_guess_extrap1") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_init_guess_extrap1 = ival;
+ }
+ else if ( strcmp(tmp[0], "cm_init_guess_extrap2") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_init_guess_extrap2 = ival;
+ }
else if ( strcmp(tmp[0], "cm_solver_pre_comp_type") == 0 )
{
ival = atoi( tmp[1] );
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index dfadda825b596470163f7b396abb01688954af1d..0abe07f208d1b617b050a6bb94af410116f17fea 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -609,6 +609,8 @@ typedef struct
real cm_solver_q_err;
real cm_domain_sparsity;
unsigned int cm_domain_sparsify_enabled;
+ unsigned int cm_init_guess_extrap1;
+ unsigned int cm_init_guess_extrap2;
unsigned int cm_solver_pre_comp_type;
unsigned int cm_solver_pre_comp_refactor;
real cm_solver_pre_comp_droptol;
diff --git a/tools/run_sim.py b/tools/run_sim.py
index c9d4e4687c011222ffaed8f935272818a77428eb..ebaaba35fc0107a3be082a5030552d638937fb3c 100644
--- a/tools/run_sim.py
+++ b/tools/run_sim.py
@@ -47,6 +47,10 @@ class TestCase():
r'(?P<key>cm_solver_q_err\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
'cm_domain_sparsity': lambda l, x: sub(
r'(?P<key>cm_domain_sparsity\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+ 'cm_init_guess_extrap1': lambda l, x: sub(
+ r'(?P<key>cm_init_guess_extrap1\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+ 'cm_init_guess_extrap2': lambda l, x: sub(
+ r'(?P<key>cm_init_guess_extrap2\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
'cm_solver_pre_comp_type': lambda l, x: sub(
r'(?P<key>cm_solver_pre_comp_type\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
'cm_solver_pre_comp_droptol': lambda l, x: sub(
@@ -252,6 +256,8 @@ if __name__ == '__main__':
'cm_solver_restart': ['100'],
'cm_solver_q_err': ['1e-6'],
'cm_domain_sparsity': ['1.0'],
+ 'cm_init_guess_extrap1': ['3'],
+ 'cm_init_guess_extrap2': ['2'],
'cm_solver_pre_comp_type': ['2'],
'cm_solver_pre_comp_refactor': ['100'],
'cm_solver_pre_comp_droptol': ['0.0'],