[cig-commits] r12075 - in long/3D/Gale/trunk: . src/Gale/plugins/SurfaceProcess
walter at geodynamics.org
walter at geodynamics.org
Tue Jun 3 02:47:16 PDT 2008
Author: walter
Date: 2008-06-03 02:47:16 -0700 (Tue, 03 Jun 2008)
New Revision: 12075
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Gale/plugins/SurfaceProcess/SurfaceProcess.c
Log:
r2208 at earth: boo | 2008-06-03 02:42:08 -0700
Make SurfaceProcess work on velocities rather than heights, making it more robust
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2207
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2208
Modified: long/3D/Gale/trunk/src/Gale/plugins/SurfaceProcess/SurfaceProcess.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/plugins/SurfaceProcess/SurfaceProcess.c 2008-06-03 09:47:12 UTC (rev 12074)
+++ long/3D/Gale/trunk/src/Gale/plugins/SurfaceProcess/SurfaceProcess.c 2008-06-03 09:47:16 UTC (rev 12075)
@@ -51,101 +51,97 @@
void Gale_SurfaceProcess_Execute( TimeIntegratee* crdAdvector, GaleContext* galeCtx ) {
Gale_SurfaceProcess_Context* spCtx;
- Mesh* surface;
- double* heights;
double dt;
-
+
assert( galeCtx );
/* Get the surface process context. */
spCtx = ExtensionManager_Get( galeCtx->extensionMgr, galeCtx, Gale_SurfaceProcess_ContextHandle );
/* Extract information from contexts. */
- surface = spCtx->surface;
- heights = spCtx->heights;
dt = galeCtx->dt;
+
/* Transfer information to the surface. */
- Gale_SurfaceProcess_Send( spCtx );
+/* Gale_SurfaceProcess_Send( spCtx ); */
/*
- ** SURFACE PROCESS CODE GOES HERE, SHOULD MODIFY THE HEIGHTS ONLY.
+ ** SURFACE PROCESS CODE GOES HERE, SHOULD MODIFY THE VELOCITIES ONLY.
*/
#if 1
{
Dictionary* dictionary = galeCtx->dictionary;
-
- double K = Dictionary_GetDouble(dictionary,
- "diffusionCoefficient");
- if( surface->topo->nDims == 1 ) {
- int n_i, nx;
- double *newVal;
+ FeVariable *velocity;
+ double K = Dictionary_GetDouble(dictionary,"diffusionCoefficient");
+
+ Stream* errorStr = Journal_Register( Error_Type, Gale_SurfaceProcess_Type );
+ Grid* grid;
+ Mesh* mesh;
+ unsigned nDims;
+
+ mesh=spCtx->mesh;
+
+ nDims = Mesh_GetDimSize( mesh );
+ grid = *(Grid**)ExtensionManager_Get( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info,
+ "vertexGrid" ) );
- nx=Mesh_GetLocalSize(surface,MT_VERTEX);
- newVal=malloc(sizeof(double)*nx);
+ velocity = (FeVariable*)FieldVariable_Register_GetByName
+ ( galeCtx->fieldVariable_Register, "VelocityField" );
+ Journal_Firewall( ( velocity!=NULL ), errorStr,
+ "Error: In KineticFriction, the name provided for the velocity field \"%s\" does not exist\n",
+ "VelocityField");
+ if( nDims == 2 ) {
+ Node_LocalIndex n_i;
+ unsigned nNodes;
+ IJK ijk;
- for( n_i = 0; n_i < nx; n_i++ ) {
- if( n_i > 0 && n_i < nx - 1 ) {
- double dx = surface->verts[n_i][0]
- - surface->verts[n_i - 1][0];
- newVal[n_i] = (heights[n_i - 1] - 2.0 * heights[n_i]
- + heights[n_i + 1]) / (dx * dx);
- }
- }
- for( n_i = 0; n_i < nx; n_i++ ) {
- if( n_i > 0 && n_i < nx - 1 ) {
- heights[n_i] = heights[n_i] + newVal[n_i]*dt*K;
- }
- }
- free(newVal);
- }
- else if( surface->topo->nDims == 2 ) {
- int n_i, n, nx, ny;
- double *newVal;
- CartesianGenerator *gen;
+ nNodes = FeMesh_GetNodeLocalSize( mesh);
+ for( n_i = 0; n_i < nNodes; n_i++ ) {
+ RegularMeshUtils_Node_1DTo3D
+ ( mesh, Mesh_DomainToGlobal( mesh, MT_VERTEX, n_i ), ijk );
+
+ if(ijk[1]==grid->sizes[1]-1
+ && ijk[0]!=0 && ijk[0]!=grid->sizes[0]-1)
+ {
+ IJK ijk_minus, ijk_plus;
+ Node_LocalIndex n_plus, n_minus;
+ double y, y_plus, y_minus, dx, delta_v, v[3];
- n=Mesh_GetLocalSize(surface,MT_VERTEX);
- gen=((CartesianGenerator *)(surface->generator));
- nx=gen->range[0]+1;
- ny=gen->range[1]+1;
+ Vec_Set2D(ijk_minus,ijk);
+ ijk_minus[0]-=1;
+ n_minus=RegularMeshUtils_Node_3DTo1D(mesh,ijk_minus);
- newVal=malloc(sizeof(double)*n);
-
- for( n_i = 0; n_i < n; n_i++ ) {
- newVal[n_i]=0;
- /* Make sure we are not on the boundary */
- if( n_i%nx!=0 && (n_i+1)%nx!=0 /* x boundary */
- && n_i>nx && n_i<n-nx) { /* y boundary */
- double dx = surface->verts[n_i][0]
- - surface->verts[n_i - 1][0];
- double dy = surface->verts[n_i][1]
- - surface->verts[n_i - nx][1];
-
- newVal[n_i] =
- (heights[n_i - 1] - 2*heights[n_i] + heights[n_i + 1])
- / (dx * dx)
- + (heights[n_i-nx] - 2*heights[n_i] + heights[n_i+nx])
- / (dy * dy);
- }
+ Vec_Set2D(ijk_plus,ijk);
+ ijk_plus[0]+=1;
+ n_plus=RegularMeshUtils_Node_3DTo1D(mesh,ijk_plus);
+
+ y_plus=mesh->verts[n_plus][1];
+ y_minus=mesh->verts[n_minus][1];
+ y=mesh->verts[n_i][1];
+
+ dx=mesh->verts[n_i][0]-mesh->verts[n_minus][0];
+ delta_v=K*(y_plus+y_minus-2*y)/(dx*dx);
+
+ FeVariable_GetValueAtNode(velocity,n_i,v);
+ v[1]+=delta_v;
+ FeVariable_SetValueAtNode(velocity,n_i,v);
+ }
}
- for( n_i = 0; n_i < nx; n_i++ ) {
- heights[n_i] = heights[n_i] + newVal[n_i]*dt*K;
- }
- free(newVal);
-
}
else {
abort();
}
+ FeVariable_SyncShadowValues(velocity);
}
#endif
/*
** END SURFACE PROCESS CODE.
*/
- /* Transfer information back to the mesh. */
- Gale_SurfaceProcess_Recv( spCtx );
+ /* Sync new velocities. */
+/* Gale_SurfaceProcess_Recv( spCtx ); */
}
Index _Gale_SurfaceProcess_Register( PluginsManager* pluginsMgr ) {
@@ -389,7 +385,7 @@
/* Append to the list of time integratee finish routines. It
should come last, because EulerDeform will reset the
values. */
- TimeIntegrator_PrependFinishEP( galeCtx->timeIntegrator,
+ TimeIntegrator_AppendSetupEP( galeCtx->timeIntegrator,
"Gale_SurfaceProcess_Execute",
Gale_SurfaceProcess_Execute,
"SurfaceProcess",
More information about the cig-commits
mailing list