[petsc-users] Back to struct in Fortran to represent field with dof > 1
Christophe Ortiz
christophe.ortiz at ciemat.es
Fri Sep 27 05:40:41 CDT 2013
Hi Barry,
I made some tests with your code that allows to use struct in Fortran90. It
compiles and does not complain with the struct with 2 components but...it
exhibits some strange behaviour. Here are the behaviours I observed. To
make it simple, I used n=10.
- After using call VecSet(x,one,ierr), I checked the values of each
component of the vector using
call VecGetArrayMyStruct(x,xarray,ierr)
do i = 0 , n
write(*,*) i , xarray(i)%a, xarray(i)%b
end do
call VecRestoreArrayMyStruct(x,xarray,ierr)
I checked the values from i= 0 to n since I was not sure what was the
convention (0- or 1-based). I obtained:
0 0.000000000000000E+000 1.00000000000000
1 1.00000000000000 1.00000000000000
2 1.00000000000000 1.00000000000000
3 1.00000000000000 1.00000000000000
4 1.00000000000000 1.00000000000000
5 1.00000000000000 1.00000000000000
6 1.00000000000000 1.00000000000000
7 1.00000000000000 1.00000000000000
8 1.00000000000000 1.00000000000000
9 1.00000000000000 1.00000000000000
10 1.00000000000000 1.996650376645798E-314
As you can see, for component a, there is a 0 for i=0, then the values are
correct. For component b, correct values go from i=1 to 9.There is garbage
in second component for i=10 . I checked that if you define a third
component c, its values go from i=-1 to n-2. It's like if values of second
component started at the end of first component, and same thing for third
component.
Then, instead of using VecSet(), I tried to assign values directly to each
component of the vector with
call VecGetArrayMyStruct(x,xarray,ierr)
do i = 0 , n-1
xarray(i)%a = 2.0
xarray(i)%b = 3.0
end do
call VecRestoreArrayMyStruct(x,xarray,ierr)
Here I checked that indices actually start at 0 and not at 1. With a loop
from i=1 to n I got a Memory corruption message.
and I checked that the values are correctly assigned with
call VecGetArrayMyStruct(x,xarray,ierr)
do i = 0 , n-1
write(*,*) i, xarray(i)%a, xarray(i)%b
end do
call VecRestoreArrayMyStruct(x,xarray,ierr)
Then, I obtained:
0 2.00000000000000 2.00000000000000
1 2.00000000000000 2.00000000000000
2 2.00000000000000 2.00000000000000
3 2.00000000000000 2.00000000000000
4 2.00000000000000 2.00000000000000
5 2.00000000000000 2.00000000000000
6 2.00000000000000 2.00000000000000
7 2.00000000000000 2.00000000000000
8 2.00000000000000 2.00000000000000
9 2.00000000000000 3.00000000000000
Here the problem seems more severe. Values are not correctly assigned to
the second component. There should be 3.0 in the second column. It seems
values of the first component are copied to the second one. Only for i=10
the value of xarray(i)%b is correct (3.0).
Any idea where it could come from ?
I guess it needs some fixes here and there but I think your idea is good
and that it could work.
Christophe
On Thu, Sep 26, 2013 at 5:45 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> I just put it with the Fortran source code and compile it with the rest
> of the application code; here is the makefile I used
>
>
>
>
> On Sep 26, 2013, at 10:14 AM, Christophe Ortiz <christophe.ortiz at ciemat.es>
> wrote:
>
> > By the way, what should I do with the small .c code ? Where should I put
> it ?
> >
> > Christophe
> >
> >
> >
> > On Thu, Sep 26, 2013 at 4:09 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > Christophe,
> >
> > Despite my LinkIn Endorsement for expertise in Fortran :-) I cannot
> pretend to be an expert in FortranXX but I have cooked up an example
> demonstrating accessing the Vec entries as if they are in an array of
> derived types. I've attached the example code; there needs to be a small C
> stub that defines the functions for your specific derived type name.
> > Note that it will only work I think if your N is a compile time constant.
> >
> >
> >
> >
> > It worked with
> > ~/Downloads$ gfortran --version
> > GNU Fortran (GCC) 4.8.1 20130404 (prerelease)
> >
> >
> > I do not understand exactly why it works since it uses
> F90Array1dCreate(fa,PETSC_SCALAR,1,len,ptr PETSC_F90_2PTR_PARAM(ptrd));
> which has a single PETSC_SCALAR as a building block but … I hope it works
> for you. If it doesn't, let us know the compiler you are using and we may
> be able to get it working for that compiler.
> >
> > Barry
> >
> >
> >
> >
> > On Sep 26, 2013, at 4:41 AM, Christophe Ortiz <
> christophe.ortiz at ciemat.es> wrote:
> >
> > > Hi all,
> > >
> > > Me again !
> > >
> > > I have read in previous posts that it is hard in Fortran to declare
> something similar to a typedef struct in C to manage a multicomponent
> problem.
> > >
> > > Is it still the case ? Has the problem been solved ?
> > >
> > > I am asking because my plan is to implement a multicomponent problem
> (in 1D), with many components that will be organized in arrays of two
> dimensions. In C I could define
> > >
> > > typedef struct{
> > > PetscScalar U;
> > > PetscScalar V;
> > > PetscScalar A[N][N];
> > > } Field;
> > >
> > > and then I could calculate the residual function with
> > >
> > > F[i].U = ...
> > > F[i].V = ...
> > > F[i].A[k][n] = ...
> > >
> > > which is quite convenient.
> > >
> > > If in Fortran it is not possible to use a struct as in C, I am afraid
> I'll have to deal with
> > >
> > > F(jdof,i) = ...
> > >
> > > where I will have only jdof to address U, V and A[ ][ ], which can be
> difficult and not very convenient I guess. Before I move all my code to C,
> does anyone have an alternative idea to manage this multi(many)component
> problem in Fortran ?
> > >
> > > Many thanks in advance for your help and suggestion !
> > >
> > > Christophe
> >
> >
> >
> >
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130927/bb668f1d/attachment.html>
More information about the petsc-users
mailing list