static PetscErrorCode TSAdaptChoose_Basic(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept,PetscReal *wlte)
{
TSAdapt_Basic *basic = (TSAdapt_Basic*)adapt->data;
PetscInt order = PETSC_DECIDE;
PetscReal enorm = -1;
PetscReal safety = basic->safety;
PetscReal hfac_lte,h_lte;
PetscErrorCode ierr;
PetscFunctionBegin;
*next_sc = 0; /* Reuse the same order scheme */
if (ts->ops->evaluatewlte) {
ierr = TSEvaluateWLTE(ts,adapt->wnormtype,&order,&enorm);CHKERRQ(ierr);
if (enorm >= 0 && order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed error order %D must be positive",order);
} else if (ts->ops->evaluatestep) {
if (adapt->candidates.n < 1) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"No candidate has been registered");
if (!adapt->candidates.inuse_set) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"The current in-use scheme is not among the %D candidates",adapt->candidates.n);
if (!basic->Y) {ierr = VecDuplicate(ts->vec_sol,&basic->Y);CHKERRQ(ierr);}
order = adapt->candidates.order[0];
ierr = TSEvaluateStep(ts,order-1,basic->Y,NULL);CHKERRQ(ierr);
ierr = TSErrorWeightedNorm(ts,ts->vec_sol,basic->Y,adapt->wnormtype,&enorm);CHKERRQ(ierr);
}
if (enorm < 0) {
*accept = PETSC_TRUE;
*next_h = h; /* Reuse the old step */
*wlte = -1; /* Weighted local truncation error was not evaluated */
PetscFunctionReturn(0);
}
/* Determine whether the step is accepted of rejected */
if (enorm > 1) {
if (!*accept) safety *= basic->reject_safety; /* The last attempt also failed, shorten more aggressively */
if (h < (1 + PETSC_SQRT_MACHINE_EPSILON)*adapt->dt_min) {
ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting because step size %g is at minimum\n",(double)enorm,(double)h);CHKERRQ(ierr);
*accept = PETSC_TRUE;
} else if (basic->always_accept) {
ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g because always_accept is set\n",(double)enorm,(double)h);CHKERRQ(ierr);
*accept = PETSC_TRUE;
} else {
ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, rejecting step of size %g\n",(double)enorm,(double)h);CHKERRQ(ierr);
*accept = PETSC_FALSE;
}
} else {
ierr = PetscInfo2(adapt,"Estimated scaled local truncation error %g, accepting step of size %g\n",(double)enorm,(double)h);CHKERRQ(ierr);
*accept = PETSC_TRUE;
}
/* The optimal new step based purely on local truncation error for this step. */
if (enorm > 0)
hfac_lte = safety * PetscPowReal(enorm,((PetscReal)-1)/order);
else
hfac_lte = safety * PETSC_INFINITY;
h_lte = h * PetscClipInterval(hfac_lte,basic->clip[0],basic->clip[1]);
*next_h = PetscClipInterval(h_lte,adapt->dt_min,adapt->dt_max);
*wlte = enorm;
PetscFunctionReturn(0);
}
PetscErrorCode PCMGMCycle_Private(PC pc,PC_MG_Levels **mglevelsin,PCRichardsonConvergedReason *reason)
{
PC_MG *mg = (PC_MG*)pc->data;
PC_MG_Levels *mgc,*mglevels = *mglevelsin;
PetscErrorCode ierr;
PetscInt cycles = (mglevels->level == 1) ? 1 : (PetscInt) mglevels->cycles;
PC subpc;
PCFailedReason pcreason;
PetscFunctionBegin;
if (mglevels->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
ierr = KSPSolve(mglevels->smoothd,mglevels->b,mglevels->x);CHKERRQ(ierr); /* pre-smooth */
ierr = KSPGetPC(mglevels->smoothd,&subpc);CHKERRQ(ierr);
ierr = PCGetSetUpFailedReason(subpc,&pcreason);CHKERRQ(ierr);
if (pcreason) {
pc->failedreason = PC_SUBPC_ERROR;
}
if (mglevels->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
if (mglevels->level) { /* not the coarsest grid */
if (mglevels->eventresidual) {ierr = PetscLogEventBegin(mglevels->eventresidual,0,0,0,0);CHKERRQ(ierr);}
ierr = (*mglevels->residual)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);CHKERRQ(ierr);
if (mglevels->eventresidual) {ierr = PetscLogEventEnd(mglevels->eventresidual,0,0,0,0);CHKERRQ(ierr);}
/* if on finest level and have convergence criteria set */
if (mglevels->level == mglevels->levels-1 && mg->ttol && reason) {
PetscReal rnorm;
ierr = VecNorm(mglevels->r,NORM_2,&rnorm);CHKERRQ(ierr);
if (rnorm <= mg->ttol) {
if (rnorm < mg->abstol) {
*reason = PCRICHARDSON_CONVERGED_ATOL;
ierr = PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n",(double)rnorm,(double)mg->abstol);CHKERRQ(ierr);
} else {
*reason = PCRICHARDSON_CONVERGED_RTOL;
ierr = PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than relative tolerance times initial residual norm %g\n",(double)rnorm,(double)mg->ttol);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
}
mgc = *(mglevelsin - 1);
if (mglevels->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = MatRestrict(mglevels->restrct,mglevels->r,mgc->b);CHKERRQ(ierr);
if (mglevels->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = VecSet(mgc->x,0.0);CHKERRQ(ierr);
while (cycles--) {
ierr = PCMGMCycle_Private(pc,mglevelsin-1,reason);CHKERRQ(ierr);
}
if (mglevels->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);CHKERRQ(ierr);
if (mglevels->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
if (mglevels->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
ierr = KSPSolve(mglevels->smoothu,mglevels->b,mglevels->x);CHKERRQ(ierr); /* post smooth */
if (mglevels->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
}
PetscFunctionReturn(0);
}
开发者ID:ziolai，项目名称:petsc，代码行数:56，代码来源:mg.c
示例3: PetscDrawSave
/*@C
PetscDrawSetSave - Saves images produced in a PetscDraw into a file
Collective on PetscDraw
Input Parameter:
+ draw - the graphics context
. filename - name of the file, if .ext then uses name of draw object plus .ext using .ext to determine the image type
- movieext - if not NULL, produces a movie of all the images
Options Database Command:
+ -draw_save <filename> - filename could be name.ext or .ext (where .ext determines the type of graphics file to save, for example .png)
. -draw_save_movie <.ext> - saves a movie to filename.ext
. -draw_save_final_image [optional filename] - saves the final image displayed in a window
- -draw_save_single_file - saves each new image in the same file, normally each new image is saved in a new file with filename/filename_%d.ext
Level: intermediate
Concepts: X windows^graphics
Notes: You should call this BEFORE creating your image and calling PetscDrawSave().
The supported image types are .png, .gif, .jpg, and .ppm (PETSc chooses the default in that order).
Support for .png images requires configure --with-libpng.
Support for .gif images requires configure --with-giflib.
Support for .jpg images requires configure --with-libjpeg.
Support for .ppm images is built-in. The PPM format has no compression (640x480 pixels ~ 900 KiB).
The ffmpeg utility must be in your path to make the movie.
.seealso: PetscDrawSetFromOptions(), PetscDrawCreate(), PetscDrawDestroy(), PetscDrawSetSaveFinalImage()
@*/
PetscErrorCode PetscDrawSetSave(PetscDraw draw,const char filename[],const char movieext[])
{
const char *savename = NULL;
const char *imageext = NULL;
char buf[PETSC_MAX_PATH_LEN];
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(draw,PETSC_DRAW_CLASSID,1);
if (filename) PetscValidCharPointer(filename,2);
if (movieext) PetscValidCharPointer(movieext,2);
/* determine save filename and image extension */
if (filename && filename[0]) {
ierr = PetscStrchr(filename,'.',(char **)&imageext);CHKERRQ(ierr);
if (!imageext) savename = filename;
else if (imageext != filename) {
size_t l1 = 0,l2 = 0;
ierr = PetscStrlen(filename,&l1);CHKERRQ(ierr);
ierr = PetscStrlen(imageext,&l2);CHKERRQ(ierr);
ierr = PetscStrncpy(buf,filename,l1-l2+1);CHKERRQ(ierr);
savename = buf;
}
}
if (!savename) {ierr = PetscObjectGetName((PetscObject)draw,&savename);CHKERRQ(ierr);}
ierr = PetscDrawImageCheckFormat(&imageext);CHKERRQ(ierr);
if (movieext) {ierr = PetscDrawMovieCheckFormat(&movieext);CHKERRQ(ierr);}
if (movieext) draw->savesinglefile = PETSC_FALSE; /* otherwise we cannot generage movies */
if (draw->savesinglefile) {
ierr = PetscInfo2(NULL,"Will save image to file %s%s\n",savename,imageext);CHKERRQ(ierr);
} else {
ierr = PetscInfo3(NULL,"Will save images to file %s/%s_%%d%s\n",savename,savename,imageext);CHKERRQ(ierr);
}
if (movieext) {
ierr = PetscInfo2(NULL,"Will save movie to file %s%s\n",savename,movieext);CHKERRQ(ierr);
}
draw->savefilecount = 0;
ierr = PetscFree(draw->savefilename);CHKERRQ(ierr);
ierr = PetscFree(draw->saveimageext);CHKERRQ(ierr);
ierr = PetscFree(draw->savemovieext);CHKERRQ(ierr);
ierr = PetscStrallocpy(savename,&draw->savefilename);CHKERRQ(ierr);
ierr = PetscStrallocpy(imageext,&draw->saveimageext);CHKERRQ(ierr);
ierr = PetscStrallocpy(movieext,&draw->savemovieext);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:plguhur，项目名称:petsc，代码行数:78，代码来源:dsave.c
示例4: Collective
/*@C
PetscHMPISpawn - Initialize additional processes to be used as "worker" processes. This is not generally
called by users. One should use -hmpi_spawn_size <n> to indicate that you wish to have n-1 new MPI
processes spawned for each current process.
Not Collective (could make collective on MPI_COMM_WORLD, generate one huge comm and then split it up)
Input Parameter:
. nodesize - size of each compute node that will share processors
Options Database:
. -hmpi_spawn_size nodesize
Notes: This is only supported on systems with an MPI 2 implementation that includes the MPI_Comm_Spawn() routine.
$ Comparison of two approaches for HMPI usage (MPI started with N processes)
$
$ -hmpi_spawn_size <n> requires MPI 2, results in n*N total processes with N directly used by application code
$ and n-1 worker processes (used by PETSc) for each application node.
$ You MUST launch MPI so that only ONE MPI process is created for each hardware node.
$
$ -hmpi_merge_size <n> results in N total processes, N/n used by the application code and the rest worker processes
$ (used by PETSc)
$ You MUST launch MPI so that n MPI processes are created for each hardware node.
$
$ petscmpiexec -n 2 ./ex1 -hmpi_spawn_size 3 gives 2 application nodes (and 4 PETSc worker nodes)
$ petscmpiexec -n 6 ./ex1 -hmpi_merge_size 3 gives the SAME 2 application nodes and 4 PETSc worker nodes
$ This is what would use if each of the computers hardware nodes had 3 CPUs.
$
$ These are intended to be used in conjunction with USER HMPI code. The user will have 1 process per
$ computer (hardware) node (where the computer node has p cpus), the user's code will use threads to fully
$ utilize all the CPUs on the node. The PETSc code will have p processes to fully use the compute node for
$ PETSc calculations. The user THREADS and PETSc PROCESSES will NEVER run at the same time so the p CPUs
$ are always working on p task, never more than p.
$
$ See PCHMPI for a PETSc preconditioner that can use this functionality
$
For both PetscHMPISpawn() and PetscHMPIMerge() PETSC_COMM_WORLD consists of one process per "node", PETSC_COMM_LOCAL_WORLD
consists of all the processes in a "node."
In both cases the user's code is running ONLY on PETSC_COMM_WORLD (that was newly generated by running this command).
Level: developer
Concepts: HMPI
.seealso: PetscFinalize(), PetscInitializeFortran(), PetscGetArgs(), PetscHMPIFinalize(), PetscInitialize(), PetscHMPIMerge(), PetscHMPIRun()
@*/
PetscErrorCode PetscHMPISpawn(PetscMPIInt nodesize)
{
PetscErrorCode ierr;
PetscMPIInt size;
MPI_Comm parent,children;
PetscFunctionBegin;
ierr = MPI_Comm_get_parent(&parent);CHKERRQ(ierr);
if (parent == MPI_COMM_NULL) { /* the original processes started by user */
char programname[PETSC_MAX_PATH_LEN];
char **argv;
ierr = PetscGetProgramName(programname,PETSC_MAX_PATH_LEN);CHKERRQ(ierr);
ierr = PetscGetArguments(&argv);CHKERRQ(ierr);
ierr = MPI_Comm_spawn(programname,argv,nodesize-1,MPI_INFO_NULL,0,PETSC_COMM_SELF,&children,MPI_ERRCODES_IGNORE);CHKERRQ(ierr);
ierr = PetscFreeArguments(argv);CHKERRQ(ierr);
ierr = MPI_Intercomm_merge(children,0,&PETSC_COMM_LOCAL_WORLD);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
ierr = PetscInfo2(0,"PETSc HMPI successfully spawned: number of nodes = %d node size = %d\n",size,nodesize);CHKERRQ(ierr);
saved_PETSC_COMM_WORLD = PETSC_COMM_WORLD;
} else { /* worker nodes that get spawned */
ierr = MPI_Intercomm_merge(parent,1,&PETSC_COMM_LOCAL_WORLD);CHKERRQ(ierr);
ierr = PetscHMPIHandle(PETSC_COMM_LOCAL_WORLD);CHKERRQ(ierr);
PetscHMPIWorker = PETSC_TRUE; /* so that PetscHMPIFinalize() will not attempt a broadcast from this process */
PetscEnd(); /* cannot continue into user code */
}
PetscFunctionReturn(0);
}
/*@
MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
zero, decreases h until this is satisfied.
Logically Collective on Vec
Input Parameters:
+ U - base vector that is added to
. a - vector that is added
. h - scaling factor on a
- dummy - context variable (unused)
Options Database Keys:
. -mat_mffd_check_positivity
Level: advanced
Notes: This is rarely used directly, rather it is passed as an argument to
MatMFFDSetCheckh()
.seealso: MatMFFDSetCheckh()
@*/
PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
{
PetscReal val, minval;
PetscScalar *u_vec, *a_vec;
PetscErrorCode ierr;
PetscInt i,n;
MPI_Comm comm;
PetscFunctionBegin;
ierr = PetscObjectGetComm((PetscObject)U,&comm);CHKERRQ(ierr);
ierr = VecGetArray(U,&u_vec);CHKERRQ(ierr);
ierr = VecGetArray(a,&a_vec);CHKERRQ(ierr);
ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr);
minval = PetscAbsScalar(*h*1.01);
for (i=0; i<n; i++) {
if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
val = PetscAbsScalar(u_vec[i]/a_vec[i]);
if (val < minval) minval = val;
}
}
ierr = VecRestoreArray(U,&u_vec);CHKERRQ(ierr);
ierr = VecRestoreArray(a,&a_vec);CHKERRQ(ierr);
ierr = MPI_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);CHKERRQ(ierr);
if (val <= PetscAbsScalar(*h)) {
ierr = PetscInfo2(U,"Scaling back h from %G to %G\n",PetscRealPart(*h),.99*val);CHKERRQ(ierr);
if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
else *h = -0.99*val;
}
PetscFunctionReturn(0);
}
开发者ID:ZJLi2013，项目名称:petsc，代码行数:52，代码来源:mffd.c
示例8: iterate
/* @ TaoApply_LineSearch - This routine takes step length of 1.0.
Input Parameters:
+ tao - TAO_SOLVER context
. X - current iterate (on output X contains new iterate, X + step*S)
. S - search direction
. f - objective function evaluated at X
. G - gradient evaluated at X
. W - work vector
. gdx - inner product of gradient and the direction of the first linear manifold being searched
- step - initial estimate of step length
Output parameters:
+ f - objective function evaluated at new iterate, X + step*S
. G - gradient evaluated at new iterate, X + step*S
. X - new iterate
- step - final step length
Info is set to 0.
@ */
static int TaoApply_UnitStep(TAO_SOLVER tao,TaoVec* X,TaoVec* G,TaoVec* S,TaoVec* W,double *f, double *f_full,
double *step,TaoInt *info2,void*ctx)
{
int info;
double fnew;
TaoVec *XL,*XU;
TaoFunctionBegin;
tao->new_search=TAO_TRUE;
info = W->CopyFrom(X); CHKERRQ(info);
info = W->Axpy(*step,S);CHKERRQ(info);
info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
if (XL && XU){
info = W->Median(XL,W,XU);CHKERRQ(info);
}
info = TaoComputeMeritFunctionGradient(tao,W,&fnew,G); CHKERRQ(info);
info = X->CopyFrom(W); CHKERRQ(info);
info = PetscInfo1(tao,"Tao Apply Unit Step: %4.4e\n",*step);
CHKERRQ(info);
if (*f<fnew){
info = PetscInfo2(tao,"Tao Apply Unit Step, FINCREASE: F old:= %12.10e, F new: %12.10e\n",*f,fnew); CHKERRQ(info);
}
*f=fnew;
*f_full = fnew;
*info2 = 0;
TaoFunctionReturn(0);
}
PetscErrorCode PetscLs(MPI_Comm comm,const char libname[],char found[],size_t tlen,PetscBool *flg)
{
PetscErrorCode ierr;
size_t len;
char *f,program[PETSC_MAX_PATH_LEN];
FILE *fp;
PetscFunctionBegin;
ierr = PetscStrcpy(program,"ls ");CHKERRQ(ierr);
ierr = PetscStrcat(program,libname);CHKERRQ(ierr);
#if defined(PETSC_HAVE_POPEN)
ierr = PetscPOpen(comm,PETSC_NULL,program,"r",&fp);CHKERRQ(ierr);
#else
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP_SYS,"Cannot run external programs on this machine");
#endif
f = fgets(found,tlen,fp);
if (f) *flg = PETSC_TRUE; else *flg = PETSC_FALSE;
while (f) {
ierr = PetscStrlen(found,&len);CHKERRQ(ierr);
f = fgets(found+len,tlen-len,fp);
}
if (*flg) {ierr = PetscInfo2(0,"ls on %s gives \n%s\n",libname,found);CHKERRQ(ierr);}
#if defined(PETSC_HAVE_POPEN)
ierr = PetscPClose(comm,fp);CHKERRQ(ierr);
#else
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP_SYS,"Cannot run external programs on this machine");
#endif
PetscFunctionReturn(0);
}
开发者ID:Kun-Qu，项目名称:petsc，代码行数:29，代码来源:ftest.c
示例10: variable
/*@
MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
zero, decreases h until this is satisfied.
Logically Collective on Vec
Input Parameters:
+ U - base vector that is added to
. a - vector that is added
. h - scaling factor on a
- dummy - context variable (unused)
Options Database Keys:
. -mat_mffd_check_positivity
Level: advanced
Notes:
This is rarely used directly, rather it is passed as an argument to
MatMFFDSetCheckh()
.seealso: MatMFFDSetCheckh()
@*/
PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
{
PetscReal val, minval;
PetscScalar *u_vec, *a_vec;
PetscErrorCode ierr;
PetscInt i,n;
MPI_Comm comm;
PetscFunctionBegin;
PetscValidHeaderSpecific(U,VEC_CLASSID,2);
PetscValidHeaderSpecific(a,VEC_CLASSID,3);
PetscValidPointer(h,4);
ierr = PetscObjectGetComm((PetscObject)U,&comm);CHKERRQ(ierr);
ierr = VecGetArray(U,&u_vec);CHKERRQ(ierr);
ierr = VecGetArray(a,&a_vec);CHKERRQ(ierr);
ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr);
minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
for (i=0; i<n; i++) {
if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
val = PetscAbsScalar(u_vec[i]/a_vec[i]);
if (val < minval) minval = val;
}
}
ierr = VecRestoreArray(U,&u_vec);CHKERRQ(ierr);
ierr = VecRestoreArray(a,&a_vec);CHKERRQ(ierr);
ierr = MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);CHKERRQ(ierr);
if (val <= PetscAbsScalar(*h)) {
ierr = PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));CHKERRQ(ierr);
if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
else *h = -0.99*val;
}
PetscFunctionReturn(0);
}
客服电话
APP下载
官方微信
请发表评论