Re: [petsc-users] problem using 64-bit-indices

2015-11-17 Thread Randall Mackie
Thanks Satish,

This fixed all error messages.

Randy


> On Nov 17, 2015, at 9:12 AM, Satish Balay  wrote:
> 
> ierr is of type 'PetscErrorCode' - not 'PetscInt'
> 
> Perhaps the dmdavecgetarrayf90 issue is also due to the difference in types 
> used for parameters
> 
> Satish
> 
> On Tue, 17 Nov 2015, Randall Mackie wrote:
> 
>> I ran into a problem yesterday where a call to DMDACreate3d gave an error 
>> message about the size being too big and that I should use 
>> —with-64-bit-indices.
>> 
>> So I recompiled PETSc (latest version, 3.6.2) with that option, but when I 
>> recompiled my code, I found the following errors:
>> 
>> 
>>  call VecGetArrayF90(vSeqZero,ptr_a,ierr)
>>   1
>> Error: Type mismatch in argument 'ierr' at (1); passed INTEGER(8) to 
>> INTEGER(4)
>> 
>>  call VecRestoreArrayF90(vSeqZero,ptr_a,ierr)
>>   1
>> Error: Type mismatch in argument 'ierr' at (1); passed INTEGER(8) to 
>> INTEGER(4)
>> 
>>  call DMDAVecGetArrayF90(da,J,ptr_j,ierr)
>> 1
>> Error: There is no specific subroutine for the generic 'dmdavecgetarrayf90' 
>> at (1)
>> 
>> 
>> So strangely, even though all my variables are defined like PetscInt, 
>> PetscReal, etc, VecGetArrayF90 is still expecting a 4 byte integer, which is 
>> easy enough for me to do, but there seems to be a problem with 
>> DMDAVecGetArrayF90. In the meantime, I’ll switch to VecGetArrayF90, and 
>> retry my code, but I thought I’d pass along these error messages.
>> 
>> 
>> Randy M.



Re: [petsc-users] problem using 64-bit-indices

2015-11-17 Thread Satish Balay
ierr is of type 'PetscErrorCode' - not 'PetscInt'

Perhaps the dmdavecgetarrayf90 issue is also due to the difference in types 
used for parameters

Satish

On Tue, 17 Nov 2015, Randall Mackie wrote:

> I ran into a problem yesterday where a call to DMDACreate3d gave an error 
> message about the size being too big and that I should use 
> —with-64-bit-indices.
> 
> So I recompiled PETSc (latest version, 3.6.2) with that option, but when I 
> recompiled my code, I found the following errors:
> 
> 
>   call VecGetArrayF90(vSeqZero,ptr_a,ierr)
>1
> Error: Type mismatch in argument 'ierr' at (1); passed INTEGER(8) to 
> INTEGER(4)
> 
>   call VecRestoreArrayF90(vSeqZero,ptr_a,ierr)
>1
> Error: Type mismatch in argument 'ierr' at (1); passed INTEGER(8) to 
> INTEGER(4)
> 
>   call DMDAVecGetArrayF90(da,J,ptr_j,ierr)
>  1
> Error: There is no specific subroutine for the generic 'dmdavecgetarrayf90' 
> at (1)
> 
> 
> So strangely, even though all my variables are defined like PetscInt, 
> PetscReal, etc, VecGetArrayF90 is still expecting a 4 byte integer, which is 
> easy enough for me to do, but there seems to be a problem with 
> DMDAVecGetArrayF90. In the meantime, I’ll switch to VecGetArrayF90, and retry 
> my code, but I thought I’d pass along these error messages.
> 
> 
> Randy M.


[petsc-users] problem using 64-bit-indices

2015-11-17 Thread Randall Mackie
I ran into a problem yesterday where a call to DMDACreate3d gave an error 
message about the size being too big and that I should use —with-64-bit-indices.

So I recompiled PETSc (latest version, 3.6.2) with that option, but when I 
recompiled my code, I found the following errors:


  call VecGetArrayF90(vSeqZero,ptr_a,ierr)
   1
Error: Type mismatch in argument 'ierr' at (1); passed INTEGER(8) to INTEGER(4)

  call VecRestoreArrayF90(vSeqZero,ptr_a,ierr)
   1
Error: Type mismatch in argument 'ierr' at (1); passed INTEGER(8) to INTEGER(4)

  call DMDAVecGetArrayF90(da,J,ptr_j,ierr)
 1
Error: There is no specific subroutine for the generic 'dmdavecgetarrayf90' at 
(1)


So strangely, even though all my variables are defined like PetscInt, 
PetscReal, etc, VecGetArrayF90 is still expecting a 4 byte integer, which is 
easy enough for me to do, but there seems to be a problem with 
DMDAVecGetArrayF90. In the meantime, I’ll switch to VecGetArrayF90, and retry 
my code, but I thought I’d pass along these error messages.


Randy M.

Re: [petsc-users] Problem Iterating DMPlex

2015-11-17 Thread Matthew Knepley
On Tue, Nov 17, 2015 at 10:11 AM, Morten Nobel-Jørgensen 
wrote:

> After distributing a DMPlex it seems like my cells are appearing twice (or
> rather multiple cells maps onto the same vertices).
>

Your code is creating the serial mesh on each process. You only want
nonzero sizes on one proc.

You can see the mesh using -dm_view, and everything using -dm_view
::ascii_info_detail

  Thanks,

 Matt


> I’m assuming the way I’m iterating the DMPlex is wrong. Essentially I
> iterate the DMPlex the following way after distribution (see code snippet
> below – or attached file).
>
> A related problem; Since distribution of a DMPlex reorders the point
> indices, how to do I map between distributed point indices and the original
> point indices.
>
> And a final question: After distributing a DMPlex, some of the vertices
> are shared and exists in multiple instances. When adding dofs to these, how
> I I know if dof is owned by the current instance or it is a ghost dof?
>
> I hope someone can point me in the right direction :)
>
> Kind regards,
> Morten
>
>
> Code snippet
>
> PetscInt from,to,dof,off;
> DMPlexGetHeightStratum(dm, 0,&from, &to);
> for (int cellIndex=from;cellIndex const PetscInt *edges;
> PetscInt   numEdges;
> DMPlexGetConeSize(dm, cellIndex, &numEdges);
> DMPlexGetCone(dm, cellIndex, &edges);
> for (int e = 0;e int edgeIndex = edges[e];
> const PetscInt *vertices;
> PetscInt   numVertices;
> DMPlexGetConeSize(dm, edgeIndex, &numVertices);
> DMPlexGetCone(dm, edgeIndex, &vertices);
> for (int v = 0;v int vertexIndex = vertices[v];
>
> For a non distibuted mesh the top of the hasse diagram looks like this:
> 0 --> 2
> 0 --> 3
> 0 --> 4
> 1 --> 4
> 1 --> 5
> 1 --> 6
>
> But when distributing (on two cores) it looks like this, where both cells
> maps to the same edges (true for both cores):
> 0 --> 11
> 0 --> 12
> 0 --> 13
> 1 --> 11
> 1 --> 12
> 1 --> 13
>
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener


[petsc-users] Problem Iterating DMPlex

2015-11-17 Thread Morten Nobel-Jørgensen
After distributing a DMPlex it seems like my cells are appearing twice (or 
rather multiple cells maps onto the same vertices).

I’m assuming the way I’m iterating the DMPlex is wrong. Essentially I iterate 
the DMPlex the following way after distribution (see code snippet below – or 
attached file).

A related problem; Since distribution of a DMPlex reorders the point indices, 
how to do I map between distributed point indices and the original point 
indices.

And a final question: After distributing a DMPlex, some of the vertices are 
shared and exists in multiple instances. When adding dofs to these, how I I 
know if dof is owned by the current instance or it is a ghost dof?

I hope someone can point me in the right direction :)

Kind regards,
Morten


Code snippet

PetscInt from,to,dof,off;
DMPlexGetHeightStratum(dm, 0,&from, &to);
for (int cellIndex=from;cellIndex 2
0 --> 3
0 --> 4
1 --> 4
1 --> 5
1 --> 6

But when distributing (on two cores) it looks like this, where both cells maps 
to the same edges (true for both cores):
0 --> 11
0 --> 12
0 --> 13
1 --> 11
1 --> 12
1 --> 13




twotriangles.cc
Description: twotriangles.cc


[petsc-users] Understanding PETSC: boundary layer flow with SNES

2015-11-17 Thread Pavel Schor

Dear PETSC users,
I am a newbie to PETSC and I try to understand it. My first attempt is 
to solve 2D laminar boundary layer equations:


u*du/dx +v*du/dy  = Ue*(dUe/dx). -nu*d2u/dy2
dudx + dvdy = 0

B.C.:
y=0:u=0,v=0
y=H:u(x)=Ue(x)
Where H is height of the domain, Ue(x) is edge velocity given by a 
potential flow solver. I assume Ue(x)=10.0



I have several questions regarding the workflow in PETSC, I took SNES 
example29 (2D cavity flow) and I tried to modify it:


1) There are vectors x and f. I suppose vector x to represent the values 
of velocity and vector f to represent velocity residuals. Please is this 
correct?


Based on example29 and my assumption, I calculated the velocity:
  dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
  dx = 1.0/dhx;   dy = 1./dhy;
u  = x[j][i].u;
v  = x[j][i].v;
dudx= (x[j][i+1].u - x[j][i-1].u)/2.0/dx;
dudy= (x[j+1][i].u - x[j-1][i].u)/2.0/dy;
dvdy= (x[j+1][i].v - x[j-1][i].v)/2.0/dy;
d2udy2= (-2.0*u + x[j-1][i].u + x[j+1][i].u)/dy/dy;

/* U velocity */
f[j][i].u  = u*dudx +v*dudy -Ue*dUedx -nu*d2udy2;

/* V velocity */
f[j][i].v  = dudx + dvdy;

The code does not work. With initial conditions x[j][i].u=Ue; and 
x[0][i].u=0.0; the result is the same as the initial conditions.


2) In SNES example29, there are boundary conditions specified on vector 
f? For example:


  /* Test whether we are on the top edge of the global array */
  if (yinte == info->my) {
j = info->my - 1;
yinte = yinte - 1;
/* top edge */
for (i=info->xs; ixs+info->xm; i++) {
f[j][i].u = x[j][i].u - lid;
f[j][i].v = x[j][i].v;

I don't understand the last two lines. I just deleted the conditions for 
left and right edges and replaced f[j][i].u = x[j][i].u - lid; with 
f[j][i].u = x[j][i].u - Ue;


3) Please could someone explain the normalization on following lines? 
(Taken from example 29)


  /*
 Define mesh intervals ratios for uniform grid.

 Note: FD formulae below are normalized by multiplying through by
 local volume element (i.e. hx*hy) to obtain coefficients O(1) 
in two dimensions.


  */
  dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal) (info->my-1);
  hx = 1.0/dhx;  hy = 1.0/dhy;
  hxdhy = hx*dhy; hydhx = hy*dhx;


Thanks in advance & Kind regards
Pavel Schor
PhD. student, Institute of aerospace engineering, Brno University of 
technology