181 |
} else { |
} else { |
182 |
dt2=dt; |
dt2=dt; |
183 |
} |
} |
184 |
while(t<dt && Paso_noError()) { |
while( (t<dt*(1.-sqrt(EPSILON))) && Paso_noError()) { |
185 |
printf("substep step %d at t=%e (step size= %e)\n",i_substeps+1,t+dt2,dt2); |
printf("substep step %d at t=%e (step size= %e)\n",i_substeps+1,t+dt2,dt2); |
186 |
Paso_FCT_setUp(fctp,dt2,sourceN,sourceP,b_n,uTilde_n,uTilde_n_coupler,QN_n,QN_n_coupler,QP_n,QP_n_coupler, |
Paso_FCT_setUp(fctp,dt2,sourceN,sourceP,b_n,uTilde_n,uTilde_n_coupler,QN_n,QN_n_coupler,QP_n,QP_n_coupler, |
187 |
options,&pp); |
options,&pp); |
198 |
* now we solve the linear system to get the correction dt: |
* now we solve the linear system to get the correction dt: |
199 |
* |
* |
200 |
*/ |
*/ |
201 |
|
/* for (i=0;i<n;++i) printf("%d %f \n", i,z_m[i]); */ |
202 |
if (fctp->theta > 0) { |
if (fctp->theta > 0) { |
203 |
omega=1./(dt2*fctp->theta); |
omega=1./(dt2*fctp->theta); |
204 |
Paso_Solver_solvePreconditioner(fctp->iteration_matrix,du_m,z_m); |
Paso_Solver_solvePreconditioner(fctp->iteration_matrix,du_m,z_m); |
279 |
double dt_max, dt_max_loc; |
double dt_max, dt_max_loc; |
280 |
register double l_ii,m; |
register double l_ii,m; |
281 |
n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix); |
n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix); |
282 |
if (! fctp->valid_matrices) { |
if ( ! fctp->valid_matrices) { |
283 |
fctp->dt_max=LARGE_POSITIVE_FLOAT; |
/* extract the row sum of the advective part */ |
284 |
/* extract the row sum of the advective part */ |
Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix); |
|
Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix); |
|
285 |
|
|
286 |
/* set low order transport operator */ |
/* set low order transport operator */ |
287 |
Paso_FCTransportProblem_setLowOrderOperator(fctp); |
Paso_FCTransportProblem_setLowOrderOperator(fctp); |
288 |
|
|
289 |
|
if (Paso_noError()) { |
290 |
/* |
/* |
291 |
* calculate time step size: |
* calculate time step size: |
292 |
*/ |
*/ |
293 |
dt_max=LARGE_POSITIVE_FLOAT; |
dt_max=LARGE_POSITIVE_FLOAT; |
294 |
if (fctp->theta < 1.) { |
#pragma omp parallel private(dt_max_loc) |
295 |
#pragma omp parallel private(dt_max_loc) |
{ |
|
{ |
|
296 |
dt_max_loc=LARGE_POSITIVE_FLOAT; |
dt_max_loc=LARGE_POSITIVE_FLOAT; |
297 |
#pragma omp for schedule(static) private(i,l_ii,m) |
#pragma omp for schedule(static) private(i,l_ii,m) |
298 |
for (i=0;i<n;++i) { |
for (i=0;i<n;++i) { |
306 |
{ |
{ |
307 |
dt_max=MIN(dt_max,dt_max_loc); |
dt_max=MIN(dt_max,dt_max_loc); |
308 |
} |
} |
|
} |
|
|
#ifdef PASO_MPI |
|
|
dt_max_loc = dt_max; |
|
|
MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm); |
|
|
#endif |
|
|
if (dt_max<LARGE_POSITIVE_FLOAT) dt_max*=1./(1.-fctp->theta); |
|
309 |
} |
} |
310 |
|
#ifdef PASO_MPI |
311 |
|
dt_max_loc = dt_max; |
312 |
|
MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm); |
313 |
|
#endif |
314 |
|
printf("dt_max = %e %e\n",dt_max, fctp->dt_factor); |
315 |
|
if (dt_max<LARGE_POSITIVE_FLOAT) dt_max*=fctp->dt_factor; |
316 |
if (dt_max <= 0.) { |
if (dt_max <= 0.) { |
317 |
Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive."); |
Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive."); |
318 |
} else { |
} else { |
319 |
if (dt_max<LARGE_POSITIVE_FLOAT) |
if (dt_max<LARGE_POSITIVE_FLOAT) |
320 |
printf("maximum time step size is %e (theta = %e).\n",dt_max,fctp->theta); |
printf("maximum time step size is %e (theta = %e).\n",dt_max,fctp->theta); |
321 |
} |
} |
322 |
fctp->dt_max=dt_max; |
fctp->dt_max=dt_max; |
323 |
fctp->valid_matrices=Paso_noError(); |
fctp->valid_matrices=Paso_noError(); |
324 |
|
} |
325 |
} |
} |
326 |
return fctp->dt_max; |
return fctp->dt_max; |
327 |
} |
} |