[cig-commits] r11890 - in long/3D/Gale/trunk: . src/StgFEM/SLE/SystemSetup/src
walter at geodynamics.org
walter at geodynamics.org
Thu May 1 16:54:26 PDT 2008
Author: walter
Date: 2008-05-01 16:54:26 -0700 (Thu, 01 May 2008)
New Revision: 11890
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c
Log:
r2135 at earth: boo | 2008-05-01 16:55:05 -0700
Decouple the BC and nonlinear iterations
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2130
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2135
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c 2008-05-01 22:59:14 UTC (rev 11889)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c 2008-05-01 23:54:26 UTC (rev 11890)
@@ -654,116 +654,120 @@
Iteration_Index maxIterations = self->nonLinearMaxIterations;
Bool converged;
Stream* errorStream = Journal_Register( Error_Type, self->type );
- double wallTime;
+ double wallTime;
Iteration_Index minIterations = self->nonLinearMinIterations;
- FeVariable* stress = NULL;
+ FeVariable* stress = NULL;
+ int bc_iterations=0;
Journal_Printf( self->info, "In %s\n", __func__ );
Stream_IndentBranch( StgFEM_Debug );
wallTime = MPI_Wtime();
-
- /* First Solve */
- self->nonLinearIteration_I = 0;
- Journal_Printf(self->info,"\nNon linear solver - iteration %d\n", self->nonLinearIteration_I);
-
- /* Add this to create a new set of bcs dependent on the stress used in friction bc's. */
+
+ /* Set up the current BC's (only really useful for friction BC's). */
stress = (FeVariable*)FieldVariable_Register_GetByName( ((FiniteElementContext *)_context)->fieldVariable_Register, "StressField" );
-
if(stress!=NULL)
{
ParticleFeVariable_Update( stress );
_SystemLinearEquations_Build(sle,_context);
}
- self->linearExecute( self, _context );
- self->hasExecuted = True;
+ /* Iterate until the boundary conditions are no longer
+ changing. This involves a full nonlinear solve each time
+ to get a real solution. */
+ do {
+ self->linearExecute( self, _context );
+ self->hasExecuted = True;
+ currentVector =
+ SystemLinearEquations_GetSolutionVectorAt( self, 0 )->vector;
+
+ /* Non linear iterations */
+ converged=0;
+ for ( self->nonLinearIteration_I = 1 ;
+ self->nonLinearIteration_I < maxIterations ;
+ self->nonLinearIteration_I++ ) {
+ Journal_Printf(self->info,"\nNon linear solver - iteration %d\n",
+ self->nonLinearIteration_I);
+
+ if(previousVector!=NULL)
+ FreeObject(previousVector);
+ Vector_Duplicate( currentVector, (void**)&previousVector );
+ Vector_SetLocalSize(previousVector,
+ Vector_GetLocalSize(currentVector));
+ Vector_CopyEntries( currentVector, previousVector );
- /* TODO - Give option which solution vector to test */
- currentVector = SystemLinearEquations_GetSolutionVectorAt( self, 0 )->vector;
-
- for ( self->nonLinearIteration_I = 1 ; self->nonLinearIteration_I < maxIterations ; self->nonLinearIteration_I++ ) {
+ self->linearExecute( self, _context );
+ self->hasExecuted = True;
+
+ /* TODO - Give option which solution vector to test */
+ currentVector =
+ SystemLinearEquations_GetSolutionVectorAt( self, 0 )->vector;
- if(previousVector!=NULL)
- _Vector_Delete(previousVector);
- Vector_Duplicate( currentVector, (void**)&previousVector );
- Vector_SetLocalSize(previousVector,
- Vector_GetLocalSize(currentVector));
- Vector_CopyEntries( currentVector, previousVector );
-
- /* Add this to create a new set of bcs dependent on the
- stress used in friction bc's. */
+ /* Calculate Residual */
+ Vector_AddScaled( previousVector, -1.0, currentVector );
+ residual = Vector_L2Norm( previousVector )
+ / Vector_L2Norm( currentVector );
+
+ Journal_Printf( self->info, "In func %s: Iteration %u of %u - Residual %.5g - Tolerance = %.5g\n",
+ __func__, self->nonLinearIteration_I,
+ maxIterations, residual, tolerance );
+ if ( self->makeConvergenceFile ) {
+ Journal_Printf( self->convergenceStream,
+ "%d\t\t%d\t\t%.5g\t\t%.5g\n",
+ self->context->timeStep,
+ self->nonLinearIteration_I,
+ residual, tolerance );
+ }
+
+ /* Check if residual is below tolerance */
+ converged = (residual < tolerance);
+
+ Journal_Printf(self->info,"Non linear solver - Residual %.8e; Tolerance %.4e%s%s - %g (secs)\n\n", residual, tolerance,
+ (converged) ? " - Converged" : " - Not converged",
+ (self->nonLinearIteration_I < maxIterations) ? ""
+ : " - Reached iteration limit",
+ MPI_Wtime() - wallTime );
+
+ if ( (converged) && (self->nonLinearIteration_I>=minIterations) )
+ break;
+ }
+
+ /* Print Info */
+ if ( converged ) {
+ Journal_Printf( self->info, "In func %s: Converged after %u iterations.\n",
+ __func__, self->nonLinearIteration_I );
+ }
+ else {
+ Journal_Printf( errorStream, "In func %s: Failed to converge after %u iterations.\n",
+ __func__, self->nonLinearIteration_I);
+ if ( self->killNonConvergent ) {
+ abort();
+ }
+ }
if(stress!=NULL)
{
ParticleFeVariable_Update( stress );
_SystemLinearEquations_Build(sle,_context);
}
+ currentVector =
+ SystemLinearEquations_GetSolutionVectorAt( self, 0 )->vector;
-
- Journal_Printf(self->info,"Non linear solver - iteration %d\n", self->nonLinearIteration_I);
+ Journal_Printf
+ (self->info,
+ "BC iteration %d - The number of elements: %d %d %d %d\n",
+ bc_iterations,
+ Vector_GetLocalSize(currentVector),
+ Vector_GetLocalSize(previousVector),
+ Vector_GetGlobalSize(currentVector),
+ Vector_GetGlobalSize(previousVector));
- self->linearExecute( self, _context );
+ ++bc_iterations;
+ } while(Vector_GetLocalSize(currentVector)
+ !=Vector_GetLocalSize(previousVector)
+ || Vector_GetGlobalSize(currentVector)
+ !=Vector_GetGlobalSize(previousVector));
- currentVector = SystemLinearEquations_GetSolutionVectorAt( self, 0 )->vector;
- /* Check whether the number of boundary elements has
- changed. I don't really like this, since it seems
- you could get it converging on some processors and
- not on others if a boundary element is lost on one
- processor and gained on another. */
-
- printf("resid %lf %lf\n",Vector_L2Norm(previousVector),Vector_L2Norm(currentVector));
- if(Vector_GetLocalSize(currentVector)
- ==Vector_GetLocalSize(previousVector)
- && Vector_GetGlobalSize(currentVector)
- ==Vector_GetGlobalSize(previousVector))
- {
- /* Calculate Residual */
- Vector_AddScaled( previousVector, -1.0, currentVector );
- residual = Vector_L2Norm( previousVector )
- / Vector_L2Norm( currentVector );
-
- Journal_Printf( self->info, "In func %s: Iteration %u of %u - Residual %.5g - Tolerance = %.5g\n",
- __func__, self->nonLinearIteration_I, maxIterations, residual, tolerance );
- if ( self->makeConvergenceFile ) {
- Journal_Printf( self->convergenceStream, "%d\t\t%d\t\t%.5g\t\t%.5g\n",
- self->context->timeStep, self->nonLinearIteration_I, residual, tolerance );
- }
-
- /* Check if residual is below tolerance */
- converged = (residual < tolerance);
-
- Journal_Printf(self->info,"Non linear solver - Residual %.8e; Tolerance %.4e%s%s - %g (secs)\n\n", residual, tolerance,
- (converged) ? " - Converged" : " - Not converged",
- (self->nonLinearIteration_I < maxIterations) ? "" : " - Reached iteration limit",
- MPI_Wtime() - wallTime );
-
- if ( (converged) && (self->nonLinearIteration_I>=minIterations) )
- break;
- }
- else
- {
- Journal_Printf(self->info,"Non linear solver - The number of boundary elements has changed %d %d %d %d: Not converged\n",
- Vector_GetLocalSize(currentVector),
- Vector_GetLocalSize(previousVector),
- Vector_GetGlobalSize(currentVector),
- Vector_GetGlobalSize(previousVector));
- }
- }
-
- /* Print Info */
- if ( converged ) {
- Journal_Printf( self->info, "In func %s: Converged after %u iterations.\n",
- __func__, self->nonLinearIteration_I );
- }
- else {
- Journal_Printf( errorStream, "In func %s: Failed to converge after %u iterations.\n",
- __func__, self->nonLinearIteration_I);
- if ( self->killNonConvergent ) {
- abort();
- }
- }
-
Stream_UnIndentBranch( StgFEM_Debug );
FreeObject( previousVector );
More information about the cig-commits
mailing list